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]

Answer 1

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