Problem with numerical solution to Sch eq in Mathematica

[naruto365]

I am trying to solve the potential barrier problem numerically in Mathematica
I am giving the following command but it is showing error which i could not figure out

a = 2;
s = NDSolve[{y1''[x] + k1*y1[x] == 0, y2''[x] + k2*y2[x] == 0,
y3''[x] + k1*y3[x] == 0}, {y1[0] == y2[0], y2[a] == y3[a],
y1'[0] == y2'[0], y2'[a] == y3'[a]}, {y1, y2, y3}, {x, -a, 2*a}]


Or is there any other way to do this numerically in mathematica

 

[Simon_Tyler]

There's a few things wrong with what you're doing.

1) The first argument of NDSolve should contain both the differential equations and the boundary/initial conditions.

2) You don't have enough boundary conditions

3) k1 and k2 have not been given numerical values.

Fixing these, something like the following works:

a = 2;
k1 = k2 = 1;
s = NDSolve[{
   y1''[x] + k1*y1[x] == 0,  y2''[x] + k2*y2[x] == 0, y3''[x] + k1*y3[x] == 0,
   y1[0] == y2[0], y2[a] == y3[a],
   y1'[0] == y2'[0], y2'[a] == y3'[a],
   y1[a] == 1, y1'[a] == 1},
  {y1, y2, y3}, {x, -a, 2*a}]

Of course, being a constant coefficient linear DE, you can also get an analytic solution.

Clear[a, k1, k2, k3]
DSolve[{
  y1''[x] + k1*y1[x] == 0,
  y2''[x] + k2*y2[x] == 0,
  y3''[x] + k3*y3[x] == 0,
  y1[0] == y2[0], y2[a] == y3[a],
  y1'[0] == y2'[0], y2'[a] == y3'[a]},
 {y1, y2, y3}, x]

The output will contain 2 integration constants C[_] that come from the fact you only have 4 conditions in 3 2nd order DEs. To solve by hand, you can vectorize the problem - or reduce to first order DEs and just use the standard matrix exponential approach.

 

[naruto365]

Thanks for the reply

1--
could you please tell me why r u putting the condition
y1[a] == 1, y1'[a] == 1
when solving numerically.

2--
again even if u do it numerically the plots of y1 ,y2, y3 look same.

 

[Simon_Tyler]

1) I told you why. You need 6 conditions to fix the 6 integration constants that come from having a system of three 2nd order differential equations. Your original code only had 4, so there is no way to numerically integrate the system. This shows up as the two leftover integration constants C[1] and C[2] in the analytic solution.

2) Your system has three non-coupled harmonic oscillators all of the same frequency:
k1=k2=1;
y1''[x] + k1*y1[x] == 0, y2''[x] + k2*y2[x] == 0, y3''[x] + k1*y3[x] == 0
Your initial/boundary conditions force y1 to be the same as y2 and y2 to be the same as y3 -- ie they are all the same.
y1[0] == y2[0], y2[a] == y3[a], y1'[0] == y2'[0], y2'[a] == y3'[a]
This is not effected by my choice of numerical initial condition.
So of course the results will be the same for y1, y2 and y3.
(Although, I didn't notice this when I wrote my previous answer - I didn't even think to check!)

If you want different results you could
- change the frequency (or wavenumber, since you use k and x) values.
- change the matching "boundary" conditions so that they no longer match.
- or add different couplings between the oscillators.

 

[Simon_Tyler]

I forgot to mention that the reiterate, as it stands, is simple enough that you can get an analytic solution. See the second code block in my first post. This might be a preferable option for you...

 

[naruto365]

Thanks again for the replyI am now learning to use mathematica for my computation and simulation(i used to do these with C++).I want to write a code in mathematica using Numerov method to solve this problem.But i need to know how to write code in mathematica. so is there any quick references for the same or i have to drill a lot.

 

[naruto365]

I have a question

instead of using
y1[a] == 1, y1'[a] == 1
as other two conditions if i will use
y1[-5]==1,y1'[5]==0
as my condition then will that work?

i am taking this condition as my wave func will be a plane wave for a<=0 and i can take a condition like above.

please tell me where did i go wrong(if any)

 

[Simon_Tyler]

It's a strange choice to make, since it's not in the range that you were numerically integrating over, but Mathematica seems to be able to handle it.
Of course, from the point of view of the analytic solution, there is no problem.
However, most numerical integrators would not be happy with such a condition - most are only happy with either initial conditions or boundary conditions.
See http://en.wikipedia.org/wiki/Boundary_value_problem

The other issue, is that you still have not differentiated y1, y2 and y3, so they will still all have the same solution...

 

[naruto365]

i modified my code a little bit like the following

a = 1;
k1 = 2;
k2 = 1;
s = NDSolve[{y1''[x] + k1*y1[x] == 0, y2''[x] + k2*y2[x] == 0,
y3''[x] + k1*y3[x] == 0, y1[0] == y2[0], y2[a] == y3[a],
y1'[0] == y2'[0], y2'[a] == y3'[a], y1[-3] == 1,
y1'[-3] == 0}, {y1, y2, y3}, {x, -10 a, 10 a}]

Plot[{y3[x] /. Out[correspoinding number]}, {x, 1, 5}]
Similarly for other two with different ranges.

but the graphs are not as expected i.e y1 and y3 should have plane wave solution and y2 should have exponentially decaying solution.
Is there anythig wrong with code.

 

[Simon_Tyler]

I don't think that Numerov's method is built into Mathematica. So you have two options.
1) Just program a numerical solver using loops like you would in C++. Using Compile[] with the CompilationTarget -> "C" option will help speed it up.
2) Write a NDSolve Method plugin: http://reference.wolfram.com/mathematica/tutorial/NDSolvePlugIns.html

Of course, Numerov's method (naively applied) only works for initial value problems. You seem to want to use various intermediate values to fix your DE's solution. This is quite a bit trickier, as you can't simply step through the solution like you can when using initial values.

As for your modified code,

a = 1; k1 = 2; k2 = 1;
s = NDSolve[{y1''[x] + k1*y1[x] == 0, y2''[x] + k2*y2[x] == 0, y3''[x] + k1*y3[x] == 0,
y1[0] == y2[0], y2[a] == y3[a], y1'[0] == y2'[0], y2'[a] == y3'[a], y1[-3] == 1, y1'[-3] == 0},
{y1, y2, y3}, {x, -10 a, 10 a}];
Plot[Evaluate[Through[{y1, y2, y3}[x]] /. s], {x, -10 a, 10 a}]

I don't know why you'd expect y2 to have a decaying solution... y2 still satisfies the basic harmonic oscillator equation.

 

[naruto365]

i actually want to calculate the reflection coefficient for any potential.so could you suggest any numerical method to do this as i want to change the parameters (like height ,length of the potential) and get the corresponding T and R.

 

[Simon_Tyler]

See the wiki article http://en.wikipedia.org/wiki/Quantum_tunneling

The easiest way is to approximate the transmission coefficient T with the WKB approximation, then for any 1D potential V(X)

$T = 1- R = \frac{e^{-2\int_{x_1}^{x_2} dx \sqrt{\frac{2m}{\hbar^2} \left( V(x) - E \right)}}}{ \left( 1 + \frac{1}{4} e^{-2\int_{x_1}^{x_2} dx \sqrt{\frac{2m}{\hbar^2} \left( V(x) - E \right)}} \right)^2}$

Or you can integrate the equations numerically and define (in 1D) $T = \frac{j_\text{trans}}{j_\text{inc}}$,
where $j_\text{inc}$ is the incident probability current and $j_\text{trans}$ is the transmitted probability current.

The incident prob current should be worked out from your initial conditions. The transmitted probability current should be numerically calculated from the solution on the opposite side of the barrier.

Note that it probably doesn't make sense to talk about transmission/reflection using anything other than an initial value problem. If you have a boundary value problem, then it's not clear which side is incident and which is transmitted. And the crazy conditions you used in your original post would be hard to make any sense of.

I suggest you try with the 1D step potential. The analytic solution is a standard textbook problem. Then compare it to a numerical solution.
Then move onto the 1D rectangular bump potential. Again, then analytic solution is standard. You can compare it to both the WKB and the numerical approaches.
Then you can move on to any harder problems you might be interested in.

 

[naruto365]

i have one problem in numerical approach that how to calculate j(inc) as how shall i know that which part is incident and which part is reflected.
I also tried the transfer matrix approach but it is giving the correct result for one rectangular potential barrier but not for two.

 

[Simon_Tyler]

You're right about j_inc. I thought using the initial conditions would be ok, but I guess I didn't think about it enough!

j_inc should be constructed out of the forward moving components of the wave on the incident side of barrier. So I guess you need to perform a Fourier decomposition on a decent sized chunk of the incident side of the barrier. You can then throw away the components with negative frequencies and transform back to get the incident wave.

This might not be the best approach and it's getting out of any area that I've had much experience with! You should read a book or lecture series on numerical quantum mechanics.

Actually, I just had a quick look at https://www.amazon.com/dp/0387979735/?tag=pfamazon01-20 and there it is suggested to use Fit[] to match up the plane waves for x<<0 and x>>0. I'm not sure whether Fit or Fourier would give better results. (Aside: I'm not a big fan of the style of Mathematica code used in this book...)

Also see the course http://www.fisica.uniud.it/~giannozz/Corsi/MQ/mq.html

If you do get some good code working for numerically calculating the scattering over a general 1D barrier, then maybe you could make a demonstration (http://demonstrations.wolfram.com) out of it.

 

[naruto365]

I am giving the following command for calculation of transmission probability through a double barrier but mathematica is not giving any solution to this infact it is taking a lot of time saying 'running'.

I want to know that whether there is something wrong with the command lines(I am sure that the integration is a valid one)

Clear[k, \[Kappa], a, b, V, \[Sigma], T, g]
\[Sigma] = 1; a = 0.1; V = 3;
\[Kappa] = Sqrt[V] - k*k;
T = Abs[(16*k*k*\[Kappa]*\[Kappa]*
      Exp[- 2 *I* k* 
        a])/((((k + \[Kappa])^2 )*
          Exp[- I \[Kappa] a] - ((k - \[Kappa])^2 )*
          Exp[ I \[Kappa] a])^2  +  
      4*((k*k - \[Kappa]*\[Kappa])^2 )*
       Exp[2 I b k]*((Sin[\[Kappa] a])^2))]^2
g = Exp[-k*k*\[Sigma]*\[Sigma]] 
Integrate[(1/(2*\[Pi]))*g*T , {k, -3*\[Sigma], 3*\[Sigma]}]

 

[naruto365]

IN THE FOLLOWING IS IT POSSIBLE TO MODIFY THE CODE SOTHAT IT WILL GIVE ME OUTPUT FOR DIFFERENT VALUE OF 'b' (LIKE ATTACHING A FOR LOOP IN THAT AS WE DO IN C++.)Clear[k, \[Kappa], a, b, V, \[Sigma], T, g]
\[Sigma] = 1; a = 0.1; b = 26; V = 3;
\[Kappa] = Sqrt[V] - k*k;
T = Abs[(16*k*k*\[Kappa]*\[Kappa]*
Exp[-2*I*k*
a])/((((k + \[Kappa])^2)*
Exp[-I \[Kappa] a] - ((k - \[Kappa])^2)*
Exp[I \[Kappa] a])^2 +
4*((k*k - \[Kappa]*\[Kappa])^2)*
Exp[2 I b k]*((Sin[\[Kappa] a])^2))]^2
g = Exp[-k*k*\[Sigma]*\[Sigma]]
NIntegrate[(1/(2*\[Pi]))*g*T, {k, -3*\[Sigma], 3*\[Sigma]}]

 

[Simon_Tyler]

Of course it's possible to modify (and there is NO NEED TO SHOUT!)

Just use a Do[] loop or a Table[] construction etc... For example:

Clear[k, \[Kappa], a, b, V, \[Sigma], T, g]
\[Sigma] = 1; a = 0.1; V = 3;
\[Kappa] = Sqrt[V] - k*k;
T = Abs[(16*k*k*\[Kappa]*\[Kappa]*
      Exp[-2*I*k*
        a])/((((k + \[Kappa])^2)*
          Exp[-I \[Kappa] a] - ((k - \[Kappa])^2)*
          Exp[I \[Kappa] a])^2 + 
      4*((k*k - \[Kappa]*\[Kappa])^2)*
       Exp[2 I b k]*((Sin[\[Kappa] a])^2))]^2;
g = Exp[-k*k*\[Sigma]*\[Sigma]];
Table[NIntegrate[(1/(2*\[Pi]))*g*T, {k, -3*\[Sigma], 3*\[Sigma]}], {b,
    2, 26,2}] // Chop

It also might be possible to do the analytic integral, or at least some approximations of it and/or series expansions if you work on it a bit.

Anyway, I don't have time to play with it at the moment (or anytime in the near future), so good luck with it!

Simon