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 :


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

Table[{x,Sqrt[x]-NIntegrate[h[y,n-1]*Sqrt[x*y],{y,0,1},AccuracyGoal->4]},{x,0,1,.1}], x];

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]

  • $\begingroup$ The only thing I can suggest is to post your actual function and see if anyone can find a fast method for integrating. $\endgroup$ Jan 28, 2017 at 15:22
  • $\begingroup$ @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. $\endgroup$
    – Rico
    Jan 28, 2017 at 15:31
  • $\begingroup$ Memoization helps a bit with the speed at the expense of memory. $\endgroup$ Jan 29, 2017 at 13:43

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_] := 
  Table[{x, Sqrt[x] - 
     NIntegrate[h[y, n - 1]*Sqrt[x*y], {y, 0, 1}, 
      AccuracyGoal -> 4]}, {x, 0, 1, .1}], x]


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] = 
   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} *)
  • $\begingroup$ Thanks a lot, that solved my problem! $\endgroup$
    – Rico
    Jan 29, 2017 at 19:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.