4
$\begingroup$

I am trying to solve a Schrödinger equation for a particle hitting a step potential using NDSolve in Mathematica. Here is my code:

mu = 6.;
m = mu;
R  = 5.;
Vs2 = 4./(2*m*R^2);
Vs = -10./(2*m*R^2) + Vs2;
Energy = 0.001
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}]

There should be decaying wave when particle hits the potential step, but NDSolve gives an increasing result. I am sure there is some trick to fix this, so I am waiting for you help.

| improve this question | |
$\endgroup$
  • 2
    $\begingroup$ Before posting this question, you should first have followed the advice in the comments to your identical question on StackOverflow: stackoverflow.com/questions/9855619/… $\endgroup$ – Jens Mar 25 '12 at 18:38
  • $\begingroup$ What is your m ? Also take a look at this (it's a bit different case but you can see how everything is set up): demonstrations.wolfram.com/ScatteringOverPotentialStep $\endgroup$ – Vitaliy Kaurov Mar 25 '12 at 18:39
  • $\begingroup$ By the way, this has exact solutions, so I'm guessing this could be a homework problem and hence not appropriate for this forum. $\endgroup$ – Jens Mar 25 '12 at 18:43
  • $\begingroup$ you forgot to define m here (this is to make it simpler for someone to cut and paste this code directly) $\endgroup$ – acl Mar 25 '12 at 18:43
  • $\begingroup$ and, the root of you problem is mathematical, not an error in programming. did you try to do it analytically, for instance? $\endgroup$ – acl Mar 25 '12 at 18:44
12
$\begingroup$

noeckel’s answer on StackOverflow is spot on. This is not a Mathematica issue, this is a mathematical issue. Namely, Mathematica is giving you the correct solution to the system of differential equation and boundary conditions given. The conditions given (and in particular the derivative imposed at the origin) are incompatible with the expected decay. Bear in mind that, at $r \geq 5$, your wavefunction will have two components of the form $\exp(\alpha r)$ and $\exp(-\alpha r)$. For each set of boundary conditions, you get a different linear combination of these two, and the only conditions that make sense are those for which the diverging term is zero.

| improve this answer | |
$\endgroup$
  • 2
    $\begingroup$ +1. If memory serves me well, usual conditions for the radial wave function are the absence of increasing exponent at infinity and finiteness at zero (or at least, convergence of the integral of Abs[f[r]]^2 over the volume around zero. The divergence of this integral would correspond to the classical fall-on-the-center scenario, for potentials ~ 1/r^a with a > 2 IIRC). This rules out fixing the derivative at zero, as you said, and fixing the w.f. value itself does not make sense either, since it is fixed by w.f. normalization (that is, unless I forgot everything, it's been a while :)) $\endgroup$ – Leonid Shifrin Mar 25 '12 at 19:24
  • $\begingroup$ sorry guys, m = mu. Also, I guess this is not exactly mathematical issue because mathematica does not realize bounded problems (f''(r)-k^2*f(r)==0 type problems), so you have to tune one of the parameters (energy, potential etc.) in order to get decaying solution. I am able to do that up to r=50-60, but I would like to get it for farther distance. $\endgroup$ – serelha Mar 25 '12 at 23:11
  • $\begingroup$ By the way, I have an analytic solution for this type problem, and my purpose is to compare my analytic result to numerical calculation. This is not homework. $\endgroup$ – serelha Mar 25 '12 at 23:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.