Forcing RecurrenceTable + NIntegrate to evaluate lazily

Question

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.

Answer 1

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
*)

Answer 2

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.