Setting boundary conditions for infinite wall potential in NDSolve

[serelha]

I am trying to use NDSolve to solve a Schrödinger equation.

My question is that if I want to solve the Schrödinger 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

 

[FunkyDwarf]

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]

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]

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]

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

 

[serelha]

Here is the Schrödinger 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]

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]

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]

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