1
$\begingroup$

I am trying to solve the Schrödinger equation for the hydrogen atom with Mathematica 10.0. I use the following code:

NDSolve[{D[R[r],{r,2}]-2 Z^2/(2n^2) R[r]+2 Z/r R[r]-l(l+1)/r^2 R[r]==0,R[0]==Sqrt[2],R'[0]==-Sqrt[2]},R,{r,0,1}]

but I always get an error 1/0 encountered; if I solve between {r,10^-10,1}, I get a solution but it is wrong. The initial conditions come from my old quantum mechanics book which says that the radial solution for Z=1, n=1 is R[r]=Sqrt[2]*Exp[-r]

$\endgroup$
  • $\begingroup$ What's the value of l? $\endgroup$ – xzczd Apr 14 at 10:19
  • $\begingroup$ @xzczd I forgot that l=0 $\endgroup$ – mattiav27 Apr 14 at 10:25
  • $\begingroup$ Something is wrong with your equation. The given solution doesn't satisfy it: D[R[r], {r, 2}] - 2 Z^2/(2 n^2) R[r] + 2 Z/r R[r] - l (l + 1)/r^2 R[r] == 0 /. R -> Function[r, Sqrt[2]*Exp[-r]] /. n -> 1 /. l -> 0 /. Z -> 1 // Simplify. $\endgroup$ – xzczd Apr 14 at 10:31
3
$\begingroup$
but I always get an error 1/0 encountered

This is because there is singularity at $r=0$. Just look at your ODE

ode = D[R[r], {r, 2}] - 2 Z^2/(2 n^2) R[r] + 2 Z/r R[r] - l (l + 1)/r^2 R[r] == 0
ic = {R[0] == Sqrt[2], R'[0] == -Sqrt[2]}

Mathematica graphics

And you are giving initial conditions at $r=0$.

In addition, you are using numerical solver and did not give any numerical values to l and Z and n, so I am not sure how you expected NDSolve to work (side note, do not use l as variable name, it looks like 1.

l = 0; Z = 1; n=1;
ode = D[R[r], {r, 2}] - 2 Z^2/(2 n^2) R[r] + 2 Z/r R[r] - l (l + 1)/r^2 R[r] == 0
ic = {R[0] == Sqrt[2], R'[0] == -Sqrt[2]}

Mathematica graphics

 NDSolve[{ode, ic}, R, {r, 0, 1}]
 (*NDSolve::ndnum Encountered non-numerical value for a derivative at r == 0.*)

One way is to push ic a little to the right of zero and hope for the best (there are better ways to find the starting point)

l = 1; Z = 1; n = 1;
eps = 1/10^6;
ode = D[R[r], {r, 2}] - 2 Z^2/(2 n^2) R[r] + 2 Z/r R[r] - l (l + 1)/r^2 R[r] == 0
ic = {R[eps] == Sqrt[2], R'[eps] == -Sqrt[2]};
NDSolve[{ode, ic}, R, {r, eps, 1}, Method -> "StiffnessSwitching"]

Mathematica graphics

 Plot[Evaluate[R[r] /. sol], {r, 0.0001, 1}]

Mathematica graphics

Or you could solve it analytically

l = 0; Z = 1; n = 1;
ode = D[R[r], {r, 2}] - 2 Z^2/(2 n^2) R[r] + 2 Z/r R[r] - l (l + 1)/r^2 R[r] == 0;
eps = 1/10^6;
ic = {R[eps] == Sqrt[2], R'[eps] == -Sqrt[2]};
DSolve[{ode, ic}, R[r], r]

Mathematica graphics

BUT

It seems your ODE is wrong. I looked it up at this site and it looks like you are solving the radial equation

Mathematica graphics

Where

Mathematica graphics

So using the above, the solution now very close to what your book says it is supposed to be Sqrt[2]*Exp[-r]. Used n=2 here which gave better approximation. It is possible also the solution Sqrt[2]*Exp[-r] is meant to be good approximation for small $r$ only, that is why it agrees well with numerical solution in the plot for small $r$.

l = 0; Z = 1; n = 2;
En = -Z^2/n^2;
ode = D[ R[r], {r, 2}] + 2/r R'[r] + (2 (En + Z/r) - l (l + 1)/r^2) R[r] == 0
eps = 10^(-6);
ic = {R[eps] == Sqrt[2], R'[eps] == -Sqrt[2]};
sol = NDSolve[{ode, ic}, R, {r, eps, 1},  Method -> "StiffnessSwitching"]

Mathematica graphics

Plot[{Sqrt[2]*Exp[-r], Evaluate[R[r] /. sol]}, {r, 0.0001, 1}]

Mathematica graphics

May be it is best to post exact ODE to solve as it is shown in book so one does not have to guess.

| improve this answer | |
$\endgroup$

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.