# Recursive Function Numerical Integration

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. Jan 28, 2017 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, 2017 at 15:31
• Memoization helps a bit with the speed at the expense of memory. Jan 29, 2017 at 13:43

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, 2017 at 19:29