Recursive Function Numerical Integration

Question

I just tried to modify a code snippet for Neumann series to the Fredholm integral equation

$ f(x)= \sqrt{x} - \int^1_0 dy\; f(y) \sqrt{x y} $

which I read in another post: How to solve an integral equation by iteration method? with a slightly different integral, namely :

f[0]:=Sqrt[x];
f[n_]:=Sqrt[x]-Integrate[(f[n-1]/.x->y)*Sqrt[x*y],{y,0,1}];
data=Table[f[i],{i,0,10}];

The code works well but for my later application I need to solve those types of integrals numerically because the Kernel is too complicated to be solved by Integrate[]. Therefore I want to replace the integration by a numerical integration in the above code and naively tried

h[x_,0]:=Sqrt[x];
h[x_,n_]:=InterpolatingPolynomial[
Table[{x,Sqrt[x]-NIntegrate[h[y,n-1]*Sqrt[x*y],{y,0,1},AccuracyGoal->4]},{x,0,1,.1}], x];
data2=Table[h[x,i],{i,0,3}];

However even for a short series up to 3 iterations the computing time "explodes" and already exceeds the first code. How can I improve the code / computation?

I tried to replace the InterpolatingPolynomial by Interpolation[]

h[x_, 0] := Sqrt[x];
hint := Interpolation[Table[{x,Sqrt[x]-NIntegrate[h[y,n-1]*Sqrt[x*y],{y,0,1},AccuracyGoal->4]},{x,0,1,.1}]];
h[x_, n_] := hint[x];
data = Table[h[x, i], {i, 0, 3}];

but that code does not even compute the first iteration step h[x,1]

The only thing I can suggest is to post your actual function and see if anyone can find a fast method for integrating. — Jack LaVigne, Jan 28 '17 at 15:22
@JackLaVigne my kernel is a rather complicated function of temperature derivatives of bose distributions and a scattering matrix and has no easy analytical solution. I think as soon as I understand how to numerically calculate the Neumann series for an easy integral like this one here with a numerical integration I can apply it to my actual problem. — Rico, Jan 28 '17 at 15:31
Memoization helps a bit with the speed at the expense of memory. — Jack LaVigne, Jan 29 '17 at 13:43

score 2 · Accepted Answer · edited Apr 13 '17 at 12:55

At the expense of memory you can get some increase in speed by using memoization.

Here is the timing using the original defintion:

h[x_, 0] := Sqrt[x]

h[x_, n_] := 
 InterpolatingPolynomial[
  Table[{x, Sqrt[x] - 
     NIntegrate[h[y, n - 1]*Sqrt[x*y], {y, 0, 1}, 
      AccuracyGoal -> 4]}, {x, 0, 1, .1}], x]

Now

AbsoluteTiming[data = Table[h[x, i], {i, 0, 3}];]

(* {5.54783, Null} *)

Memorize each step as you compute

hr[x_, 0] = Sqrt[x]

hr[x_, n_] := 
 hr[x, n] = 
  InterpolatingPolynomial[
   Table[{x, Sqrt[x] - 
      NIntegrate[hr[y, n - 1]*Sqrt[x*y], {y, 0, 1}, 
       AccuracyGoal -> 4]}, {x, 0, 1, .1}], x]

The timing is much better

AbsoluteTiming[data = Table[hr[x, i], {i, 0, 3}];]

(* {0.293722, Null} *)

After clearing everything here is the time to compute the table with i set to 10.

AbsoluteTiming[data = Table[hr[x, i], {i, 0, 10}];]

(* {1.35537, Null} *)

Thanks a lot, that solved my problem! – Rico Jan 29 '17 at 19:29 — Rico, Jan 29 '17 at 19:29

Recursive Function Numerical Integration

1 Answers1