4

When i try to plot the numerical solution of quantum harmonic oscillator , it blows up even for the correct energies. Why does this happen and how can it be fixed? The code i used is as following.

sol = 
  First[
    NDSolve[
      {-1/2 D[ψ[x], {x, 2}] + 1/2 x^2 ψ[x] == (1 + 1/2) ψ[x], 
       ψ[0] == 0, ψ'[0] == 1}, 
      ψ[x], {x, -a, a}], 
      MaxSteps -> 10000000] 

Plot[ψ[x] /. sol, {x, -a, a}, PlotRange -> 1]
m_goldberg
  • 107,779
  • 16
  • 103
  • 257
Shubham Raghuvanshi
  • 91
  • 1
  • What is the numerical value for a? – zhk Dec 22 '16 at 08:04
  • typically it is 10 but you can play around with it – Shubham Raghuvanshi Dec 22 '16 at 08:07
  • 1
    The problem is that you're trying to get the solution that's decaying very sharply (Gaussian times a Hermite polynomial); increasing WorkingPrecision will delay this, but you'll hit the effect at even higher values of a. BTW, this ODE can be solved by DSolve[]. – J. M.'s missing motivation Dec 22 '16 at 08:10
  • Yes it can be solved by DSolve but i wanted to write c programm to solve it numerically and for that is used NDSolve just to see what kind of behavior is expected from numerical solutions. Can i change the method of solving this ODE in mathematica to get better results. – Shubham Raghuvanshi Dec 22 '16 at 08:20
  • 1
    It's kind of an inherent property of such DEs: when your DE has a growing and a decaying solution, the result of a numerical integration method is actually some linear combination of those solutions, where the coefficients multiplying the growing solution are hopefully small. So, a numerical solution will be fine over a small-ish interval, but will certainly blow up past it. This is even more of a problem if you'll be doing it in C, since you don't have arbitrary precision there (usually). – J. M.'s missing motivation Dec 22 '16 at 08:26
  • Note that you've set MaxSteps-> inside of First, not NDSolve[], you might also be interested in NDSolveValue[] – Feyre Dec 22 '16 at 08:54
  • , WorkingPrecision -> 50 gives you an accurate plot for the entire domain given. But like @J.M. said, this won't help you in C. – Feyre Dec 22 '16 at 08:56
  • 1
    I'm voting to close this question as off-topic because the question has no solution other than the partial fixes mentioned in the comments. – Feyre Dec 22 '16 at 10:01

0 Answers0