Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

How to obtain high quantum harmonic eigenfunctions? Everything works until, let say 65th eigenfunction, but for higher eigenfunctions values are too high. My code

NoH = Table[1/Sqrt[2^n n!] Pi^(-1/4), {n, 0, ns - 1}];
phi[x_] = Table[NoH[[i]] HermiteH[i - 1, x] Exp[-x^2/2], {i, ns}];

NoH are normalization factors, phi[x] are eigenfunctions.

Best regards.

share|improve this question
    
I don't see what the problem is. Could you add detail? Also, why does it matter than they're eigenfunction to the QHO? this is irrelevant. –  acl Mar 14 '13 at 17:14
4  
Chances are that you get into trouble using machine precision, since for higher eigenfunctions massive cancellations of large quantities lead to accumulation of errors. Try to increase the precision, and, generally, work at fixed precision - then Mathematica will track the precision for you and you will see how reliable your results are. –  Leonid Shifrin Mar 14 '13 at 17:17
1  
Also, it would be best to include complete self-contained examples, so that we don't have to guess what your real problem is. –  Leonid Shifrin Mar 14 '13 at 17:19
1  
But that's what I mean: those oscillations are constructed by a high-degree polynomial plus the exponent - so the right oscillatory behavior is a result of non-trivial cancellation of individual terms in the polynomials (and exponent), which are quite large for large x. –  Leonid Shifrin Mar 14 '13 at 17:22
1  
At some point, you may actually be better off, in terms of computations, by using the WKB - approximated functions, for large values of the quantum number (where WKB works quite well). –  Leonid Shifrin Mar 14 '13 at 17:24
show 6 more comments

1 Answer

It's very easy to get around the numerical problem:

Instead of asking Mathematica to generate a list of 100 or so function values at every value of the independent variable x, just define the functions individually for each index:

phi1[n_][x_] := 1/Sqrt[2^n n!] Pi^(-1/4) HermiteH[n - 1, x] Exp[-x^2/2]

Plot[phi1[99][x], {x, 0, 20}, PlotRange -> All]

plot

Plot[phi1[200][x], {x, 0, 20}, PlotRange -> All]

plot200

In this way, the numerical accuracy is tracked individually for each index, not for the whole table.

Edit: what about tables of functions?

From the fact that you defined a table, one can guess that you wanted to plot a set of these functions in a single plot. In that case, you again have the numerical stability problem, so it's best to separate the individual plots and combine them with Show:

Show[Table[
  Plot[phi1[n][x], {x, 0, 21}, PlotRange -> All, 
   PlotStyle -> ColorData[1][n]], {n, 100, 150, 10}]]

Show

Here, I used ColorData[1] as the default color scheme that you also get when plotting a table of functions. The PlotRange is determined by the first of the plots in Show (but since I chose All it will make everything fit into the frame), but you can change that explicitly if needed.

I had some custom style options set when making these plots, so yours will have no frame and a different font by default...

share|improve this answer
    
Something I don't understand: even with the function defined this way, the problem arises if the x in the expression isn't local to the Plot. e.g. expr = phi1[99][x]; Plot[expr, {x,0,20}] In fact the OP's code works if phi is defined using SetDelayed. –  Simon Woods Mar 14 '13 at 21:17
    
@SimonWoods It looks like the numerical errors appear whenever the high-degree polynomials are passed to Plot in expanded form. When Plot sees that they are HermiteH polynomials, it does the right thing. But since we can't look inside the Plot code, that's all I can hypothesize. –  Jens Mar 14 '13 at 21:52
    
Interesting that when we define phi as a function of two variables it works, and when we define it as a table it doesn't. Thank you. –  weisskreuz Mar 14 '13 at 22:12
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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