Group Abstract Group Abstract

Message Boards Message Boards

0
|
6.8K Views
|
6 Replies
|
0 Total Likes
View groups...
Share
Share this post:

Schrodinger Equation by Shooting Method

Sara Zamharir
Posted 10 years ago

Hello everybody

I want to solve the radial part of Schrodinger equation by shooting method. In this method I set the initial value of energy and appropriate boundary conditions. By varying the value of energy the best wave function, which satisfies the boundary conditions, is found. My code is:

ClearAll;
ClearSystemCache;
v[x_]:=If[0<x< 10,0,100];
energy=13.3;
xMax=12;
mc2=0.511 * 10^9;
hc = 1240 * 10^3 ;
solution = NDSolve[{x psi''[x]+ 2 psi'[x]== (8 Pi2 (mc2))/hc2 x ( v[x] - energy) psi[x], psi[0.01] == 1, psi'[0.01] == 0.001}, psi, {x, 0.01, xMax}];
Plot[psi[x]/. solution, {x, 0.01, xMax}]

Now I want to make a loop to change the value of energy and again recalculate the wave function until the correct wave function is obtained. I used Do[expr,{i,imin,imax,di}] but it didn’t work. Please guide me how to make a loop to change the energy and plot wave functions for each value of energy.
POSTED BY: Sara Zamharir
6 Replies
Sort By:

The PlotLegends option for your plot function may do what you want.
POSTED BY: Frank Kampas

It would be easier to read your code if you put it into a "code sample" block, which is the first icon.

I suggest that you first test your solution[energy] function to see if it gives the correct result for the correct energy. Also note that the code I posted as an illustration can be split up to get the value found by FindRoot. Since there are multiple eigenvalues, you need to give FindRoot starting values near the correct result.

In[1]:= shootsoln[d_?NumberQ] := 
 NDSolve[{y''[x] - 5 Sinh[5 y[x]] == 0, y[0] == 0, y'[0] == d}, 
   y, {x, 0, 1}][[1]]

In[2]:= yd[d_?NumberQ] := y[1] /. shootsoln[d]

In[3]:= fr = FindRoot[yd[d] == 1, {d, -.05, .05}]

Out[3]= {d -> 0.0457504}

In[4]:= soln = y /. (shootsoln[d] /. fr)

Out[4]= InterpolatingFunction[{{0., 1.}}, <>]
POSTED BY: Frank Kampas
Sara Zamharir
Sara Zamharir
Posted 10 years ago

Dear Dr. Kampas

Thanks a million for all your hints and suggestions. I used ParametricNDSolve function and fortunately it worked. now I want to have legend on the graph to label wave function corresponding to each energy value. how should I do that? is that PlotLegend function?
POSTED BY: Sara Zamharir

Here's an example of how it could be done. I wrote this before ParametricNDSolve was available.

In[1]:= shootsoln[d_?NumberQ] := 
 NDSolve[{y''[x] - 5 Sinh[5 y[x]] == 0, y[0] == 0, y'[0] == d}, 
   y, {x, 0, 1}][[1]]

In[2]:= yd[d_?NumberQ] := y[1] /. shootsoln[d]

In[4]:= soln = y /. (shootsoln[d] /.
           FindRoot[yd[d] == 1, {d, -.05, .05}][[1]])

Out[4]= InterpolatingFunction[{{0., 1.}}, <>]
POSTED BY: Frank Kampas
Sara Zamharir
Sara Zamharir
Posted 10 years ago

Dear Dr. Kampas Thanks for your reply. My new code is:

CearAll; ClearSystemCache; v[x_] := If[0 < x < 10, 0, 100]; xMax = 12; mc2 = 0.511 * 10^9; hc = 1240 * 10^3 ; system [energy_] := { psi''[x] + 2 psi'[x] == (8 Pi^2 (mc2))/ hc^2 x ( v[x] - energy) psi[x], psi[0] == 1, psi'[0] == 0.001}; solution[energy_] = NDSolve[system [energy], psi[x], {x, 0, xMax}]; yend[energy_?NumericQD] := First[psi[x] /. solution [energy] /. x]; bc = FindRoot[yend [energy] == 0, {energy, 13, 13.5}]; Plot[ Evaluate[psi[x] /. solution[energy /. bc]], {x, 0, xMax}]

?again it doesn't give the correct answer. it should be something like the attached file. one more thing is that this approach doesn't give us the energy value for which the wave function has been plotted. What should I do?

Attachments:
POSTED BY: Sara Zamharir

I think it would be better to use FindRoot .
POSTED BY: Frank Kampas
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard