Help with Mathematica NDSolved


by serelha
Tags: mathematica, ndsolved
serelha
serelha is offline
#1
Mar7-12, 03:04 PM
P: 5
I am trying to use NDSolve to solve a Schrodinger equation.

My question is that if I want to solve the Schrodinger equation for an infinite wall potential at r = R, I would want my wave function to die at r = R.

Here is the code that I know for a square well potential;


NDSolve[ f''[r] + 2*mu*Energy - 2*mu*VOO[r]f[r] == 0, f[0] == 1.,
f'[0] == 0.0]


However, I do not know how to define specially boundary conditions for the finite wall problem.

If you have any suggestion, I will appreciate it.


Thanks
Phys.Org News Partner Science news on Phys.org
Going nuts? Turkey looks to pistachios to heat new eco-city
Space-tested fluid flow concept advances infectious disease diagnoses
SpaceX launches supplies to space station (Update)
FunkyDwarf
FunkyDwarf is offline
#2
Mar17-12, 10:48 PM
P: 479
You can do it two ways.

You can either assume a solution in the exponentially decaying region and match that to your function in the 'well' part of the potential (i.e. this will define some boundary conditions at your finite wall)

OR

you can solve the whole thing numerically but then you need to think about a method to find the energy eigenvalues (i.e. you will probably need to vary the energy in order to get exponentially decaying solutions in the appropriate region, rather than decaying and increasing)
serelha
serelha is offline
#3
Mar21-12, 10:14 AM
P: 5
Quote Quote by FunkyDwarf View Post
You can do it two ways.

You can either assume a solution in the exponentially decaying region and match that to your function in the 'well' part of the potential (i.e. this will define some boundary conditions at your finite wall)

OR

you can solve the whole thing numerically but then you need to think about a method to find the energy eigenvalues (i.e. you will probably need to vary the energy in order to get exponentially decaying solutions in the appropriate region, rather than decaying and increasing)
Thanks for your help.

I want to solve the whole system numerically, and I have already tried the method that you suggested. In this problem since I want to calculate the phase shift, I can vary the energy as much as I want.

Here is a sample code


mu = 6.;
width = 5.;
R = width;

Vb = 4./(2*mu*R^2);
Vs = -10./(2*mu*R^2) + Vb;

p = .3;
Energy = p^2/(2*mu);


VCC[r_] =
Vs*UnitStep[width - r] +Vb*UnitStep[r - width];
eps = .01;

"Numerical Results"
system1 = {RC''[r] + (-2*mu*(VCC[r] - Energy))*RC[r] == 0,
RC[eps] == 1.0, RC'[eps] == 0.0 };



syssol1 =
NDSolve[system1, { RC[r]}, { r, eps, 50./p}, MaxSteps -> 10000000];

Plot[Evaluate[{ RC[r]} /. syssol1], {r, eps, 10.0/p},
PlotRange -> {-100, 100}]

No matter how you vary p parameter here, you always get an increasing solution. However, what I want a decaying around R, which is physical.

FunkyDwarf
FunkyDwarf is offline
#4
Mar21-12, 10:36 PM
P: 479

Help with Mathematica NDSolved


I think for starters your boundary conditions are wrong. If you are using a second order DE with no first order derivative that means at some point you made the substitution R = u(r)/r, this implies your near zero solution should be zero:

R(r) ~ r^l
u(r) = r^l+1 so near r = 0, u(r) = 0 and u'(r) = 1.

Also you will never get the growing solution equal to exactly zero unless you are exactly at the energy eigenvalue. You can test this with your code by modeling the hydrogen case and using known analytic solutions to compare with.
serelha
serelha is offline
#5
Mar21-12, 11:05 PM
P: 5
Thanks for your respond..

Regarding the boundary conditions;

my original differential equation is Bessel differential equations, so as you guessed I did substitution R[r]=RC[r]/r, however, the boundary conditions are correct for L = 0 (s-wave). The conditions that you suggested are correct for L>0 solutions (p, d, f...waves). As I said I tried to give a simple example to discuss.

Regarding the solution with p=0.4;

functions becomes constant for r=>R, but not zero as you could see in the attachment. Do you think this corresponds to the physical solution, I am skeptical about it.

Thanks
Attached Thumbnails
Untitled-1.jpg  
serelha
serelha is offline
#6
Mar24-12, 05:00 PM
P: 5
Here is the Schrodinger equation I mentioned

Manipulate[
m = 6.;
R = 5.;

Vs2 = 4./(2*m*R^2);
Vs = -10./(2*m*R^2) + Vs2;


VCC[r_] = Vs*UnitStep[R - r] + Vs2*UnitStep[r - R];

L = 0;

system = {RC''[r] +
2/r*RC'[r] + (-L*(L + 1)/r^2 - 2*mu*(VCC[r] - Energy))*RC[r] ==
0, RC[0.001] == 1.0, RC'[0.001] == 0.0 };


syssol =
NDSolve[system, {RC[r]}, { r, 0.001, 1000.}, MaxSteps -> 10000000];

Plot[Evaluate[{RC[r]} /. syssol], {r, 0.001, 200.0},
PlotRange -> {-1.1, 1.1}]
, {Energy, 0.001, 0.012, 1/100}]

The value that you point out p = 0.4 corresponds an energy value the same as the barrier, so the particle is not bounded for that value.
FunkyDwarf
FunkyDwarf is offline
#7
Mar27-12, 05:36 PM
P: 479
See Landau Lifgarbagez Vol 3, section 32:

If you have the full radial Schrod equation (where you've not done any transformations) the solution near the origin (for a potential V(r) such that r^2 V(r) vanishes near the origin) is R ~ r^l.
serelha
serelha is offline
#8
Mar27-12, 05:41 PM
P: 5
You are right, but L=0 here. So the boundary conditions are correct. If L>0, then the w.f vanishes at the origin.
FunkyDwarf
FunkyDwarf is offline
#9
Mar29-12, 12:52 AM
P: 479
Yes but your two equations are inconsistent, in the step well case the equation appears transformed (no first derivative term), in the other one it is not yet the BCs are the same, thus something is wrong.

In any event the real answer (which I missed because i'm currently working on a relativistic problem and that's messing with my head a bit) is that you are only looking at positive energies.

If you want bound state solutions use this code:

Manipulate[mu=6.;

R=5.;
Vs2=4./(2*mu*R^2);
Vs=-10./(2*mu*R^2)+Vs2;
VCC[r_]=Vs*UnitStep[R-r]+Vs2*UnitStep[r-R];
L=0;
system={RC''[r]+2/r*RC'[r]+(-L*(L+1)/r^2-2*mu*(VCC[r]-Energy))*RC[r]==0,RC[0.001]==1.0,RC'[0.001]==0.0};
syssol=NDSolve[system,{RC[r]},{r,0.001,1000.},MaxSteps->10000000];
Plot[Evaluate[{RC[r]}/.syssol],{r,0.001,200.0},PlotRange->{-1.1,1.1}],{Energy,-.1,0.12,1/10000}]

You also had mu and m as different things but I assume they are the same as m lacked a numeric definition so i fixed that.

It's a rather finicky system: your ground state is somewhere between E = -0.0021 and E = -0.002


Register to reply

Related Discussions
Mathematica 7 Math & Science Software 11
mathematica help Calculus & Beyond Homework 1
Mathematica code validator Math & Science Software 1
athematica seems to be asking for more arguments Math & Science Software 2
Mathematica help Math & Science Software 2