# Forcing RecurrenceTable + NIntegrate to evaluate lazily

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.

-

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}]


(*
0.0280796
*)

-

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}]


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.

-