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.

I have a recurrence relation that grows exponentially if expressed in closed-form. I need to integrate this function. If I plot it, I think Mathematica evaluates the recurrence numerically at every level, so the recurrence doesn't grow to be a huge expression. NIntegrate, however, seems to be trying to find a closed-form equation to integrate, even when I turn off symbolicprocessing. I've tried all the options I can think of, including extracting my relation to a function, adding Simplify to the recurrence, and I'm just completely stumped at how I can force Mathematica to stop using symbolic manipulation for an expression that has $2^{100}$ terms. Anyone have any ideas?

HOWavefunctions = Function[{q, n}, 
    RecurrenceTable[{y[i + 1] == 
       Sqrt[2/(i + 1)]*q*y[i] - Sqrt[i/(i + 1)]*
           y[i - 1], y[0] == 1/(Pi^0.25*E^(0.5*q^2)), 
     y[1] == Sqrt[2]*q*y[0]}, y, {i, 0, n}]]
n = 100;
Plot[HOWavefunctions[q, n][[n + 1]]^2, {q, 14.05, 15}]
NIntegrate[HOWavefunctions[q, n][[n + 1]]^2, {q, 14.05, 15}]

(If you turn $n$ down to $15$ or so, both work just fine.)

I guess, ultimately, what I'm looking for is a way to tell RecurrenceTable to evaluate 'lazily', that is, only evaluate a term once all variables have numeric values, and never ever expand the relationship symbolically. I tried inserting a Simplify in the recurrence, but Mathematica apparently ignored me.

share|improve this question

2 Answers 2

HOWavefunctions[q_?NumericQ, n_] := 
 RecurrenceTable[{y[i + 1] == Sqrt[2/(i + 1)]*q*y[i] - Sqrt[i/(i + 1)]*y[i - 1], 
                y[0] == 1/(Pi^0.25*E^(0.5*q^2)), y[1] == Sqrt[2]*q*y[0]}, y, {i, 0, n}]
n = 100;

Plot[HOWavefunctions[q, n][[n + 1]]^2, {q, 14.05, 15}]
NIntegrate[HOWavefunctions[q, n][[n + 1]]^2, {q, 14.05, 15}, 
           Method -> {Automatic, "SymbolicProcessing" -> None}]

Mathematica graphics

(*
 0.0280796
*)
share|improve this answer

I present this as an alternative to evaluation of the function (exploiting the insights of belisarius).

 fun[q_?NumericQ, n_] :=
 Last@Fold[{{0, 1}, {-Sqrt[#2/(#2 + 1)], 
       Sqrt[2/(#2 + 1)] q}}.#1 &, {1/(Pi^0.25*E^(0.5*q^2)), 
    Sqrt[2]*q*1/(Pi^0.25*E^(0.5*q^2))}, Range[1, n]]

I will use h for the HOWavefunctions. Note the indexing is slightly different (viz, h[q,n]=fun[q,n-1]).

h[q_?NumericQ, n_] := 
 RecurrenceTable[{y[i + 1] == 
    Sqrt[2/(i + 1)]*q*y[i] - Sqrt[i/(i + 1)]*y[i - 1], 
   y[0] == 1/(Pi^0.25*E^(0.5*q^2)), y[1] == Sqrt[2]*q*y[0]}, 
  y, {i, 0, n}]
n = 100;

Timing:

Timing[Do[h[14.05, 100], {1000}]]
Timing[Do[fun[14.05, 100], {1000}]]

give: 0.66, 0.063 respectively

Visualizing function to be integrated:

g[w_?NumericQ, r_] := fun[w, r]^2;
Plot[g[s,99],{s,14.05,15}]

enter image description here

Integrating:

NIntegrate[g[s, 99], {s, 14.05, 15}, 
 Method -> {Automatic, "SymbolicProcessing" -> None}]

yields:

(* 0.0280796 *)

You could also (depending on your precision requirment and mesh choice) generate a table interpolate and integrate interpolating function.

share|improve this answer

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.