# Integral equation numerical solution with NDSolve

I'm trying to solve something like: f[x] == Integrate[f[x]*g[x]] where g[x] is known and f[x] is the function to solve (numerically is fine).

I've seen some examples of integral equations but wonder if there is a direct way to solve it with Mathematica.

• NDSolve[] is meant for differential equations, and there isn't a built-in function (yet) for solving integral equations. OTOH, solutions to an inhomogeneous Fredholm equation of the second kind, like in your example, can be solved with the Liouville-Neumann series; such an expansion ought to be doable with the built-in functions of Mathematica. – J. M. is in limbo Oct 5 '12 at 10:46
• Would you care to comment about the field in which you've found this problem? (seems pretty difficult!) – Dr. belisarius Oct 5 '12 at 11:49
• Oh dear, your actual problem looks nothing like f[x] == Integrate[f[x] g[x], x] (thus, your "a bit more complex" is one hell of an understatement). This is a nonlinear problem, whose solution is much more difficult... – J. M. is in limbo Oct 5 '12 at 12:04
• PCL, thanks for registering an account. I merged it with the unregistered one so you should now be able to edit your own posts without the edit needing approval. – Mr.Wizard Oct 5 '12 at 12:51
• There are a few examples of numerically solving integral eqns in Mathematica.StackExchange.com and the older StackOverflow version. A search should locate them. That said, I do not know if any will help for your particular problem. – Daniel Lichtblau Oct 5 '12 at 15:42

If you want to solve the Fredholm equation of the second kind which is an integral equation of the form

$$f(x) - \lambda\int_{a}^{b} K(x,y)f(y)dy = g(x), \quad \forall x \in [a, b]$$ you can use the following code. We follow the simplistic source. One can make it much better by considering further implementation aspects as well as injecting more state of the art algorithmic antibiotic!!

## Code:

Options[FredholmKind2] = {Method -> Automatic};
FredholmKind2[{a_, b_, lambda_, k_, g_}, n_?IntegerQ,OptionsPattern[]] :=
Block[{step, SI, GI, KMatrix, W, DMatrix, f, deltaX, delta},
step = (b - a)/n;
SI = Range[a, b, step];
GI = g /@ SI;
KMatrix = Outer[k, SI, SI];
W = {step/2}~Join~ConstantArray[step, n - 1]~Join~{step/2};
DMatrix = DiagonalMatrix[W];
f = If[OptionValue[Method] === NIntegrate,
deltaX[x_?NumericQ] := W . (k[x, #] & /@ SI) -
NIntegrate[k[x, y], {y, a, b},AccuracyGoal -> 4];
(*If the integral is expensive ParallelMap is an option here *)
delta = deltaX /@ SI;
Interpolation[
Transpose@{SI,
LinearSolve[IdentityMatrix[n + 1] + lambda*(DiagonalMatrix[delta] -
KMatrix . DMatrix), GI]}]
,
Interpolation[
Transpose@{SI,
LinearSolve[IdentityMatrix[n + 1] - lambda*(KMatrix . DMatrix),GI]}]
];
f
]


## Testing:

Now lets test it for the following equation which has an exact solution Sin[x]!! $$f(x) - \frac{4}{\pi}\int_{0}^{\frac{\pi}{2}} \cos(x-y)f(y)dy = -\frac{2}{\pi}\cos(x), \quad \forall x \in \left[0, \frac{\pi}{2}\right]$$ There are two separate methods available called Automatic and NIntegrate. The second method is computationally more expensive but produces better result. The above function returns the solution f as an InterpolatingFunction which you can later use an ordinary function in MMA.

n=90;(*number of discretization*)
a = 0.;
b = 0.5*Pi;
lambda = 4./Pi;
Kpart[x_, y_] := Cos[x - y];
Gpart[x_] := -2. Cos[x]/Pi;
f1 = FredholmKind2[{a,b,lambda,Kpart, Gpart},n,Method -> Automatic];
f2 = FredholmKind2[{a,b,lambda,Kpart, Gpart},n,Method -> NIntegrate];
Needs["PlotLegends"];
Plot[Evaluate@Sqrt@((Sin[x] - #)^2 & /@ {f1[x], f2[x]}), {x, a, b},Frame -> True,
Axes -> False,PlotLegend -> {"Automatic", "NIntegrate"},LegendPosition -> {1.1, -0.4}] In the above plot one can see how accurate the numerical solutions perform w.r.t the exact solution.

## To Do!

Now if you want to solve the Fredholm integral equation of the first kind which looks like

$$\int_{a}^{b} K(x,y) f(y) dy = g(x), \quad \forall x \in [a,b],$$

I suggest that you take look at this at page 213. After the above example it will be easy for you to implement the solution for the problem of first kind.

BR

• Awake dead thread! @PlatoManiac, should you expect things to go badly for homogeneous equations (that is, with g(x) = 0)? – Ian Nov 30 '12 at 16:59
• @plato Please take a look at this one mathematica.stackexchange.com/q/58149/193 – Dr. belisarius Aug 25 '14 at 17:14
• @belisarius Done! Thx for notifying though..:) – PlatoManiac Aug 25 '14 at 20:50
• @PlatoManiac, I tried applying this solution for Kpart[x_, y_] := Piecewise[{{1/Abs[x - y]^(8/3), Abs[x - y] >= u}, {0, Abs[x - y] < u}}]; and Gpart[x_] := Piecewise[{{-1/2 + 1/2 (u/(1 - x))^(5/3), x < u}, {-1/2 (u/x)^(5/3) + 1/2, 1 - x < u}, {-1/2 (u/x)^(5/3) + 1/2 (u/(1 - x))^(5/3), u <= x <= 1 - u}}]; for u <<1 and for a=0 and b=1. However, the solution doesn't satisfy the equation for x<u or x>1-u. I guess the problem is in the use of the Trapezoidal rule? Do you perhaps have a workaround for these kernel and source terms? – Asaf Miron Jan 6 '19 at 13:49

I was not going to post this because in a previous edit (now gone) you showed a very difficult problem.

Anyway, the current question could be answered as follows:

You want to find f[x] satisfying

f[x] == Integrate[f[x] g[x], x]


for a known g[x].

Differentiating:

f'[x] == g[x] f[x]


So for example

g[x_] := Sin@x^2;
fs = DSolve[f'[x] == g[x] f[x] && f == 1, f, x]
(*
->{{f -> Function[{x}, E^(x/2 - 1/4 Sin[2 x])]}}
*)


Testing it

Integrate[g[u] (f /. fs[])[u], u]
(*
-> E^(u/2 - 1/4 Sin[2 u])
*)

• Sometimes, one is lucky, and integral equations (as well as integro-differential equations) can be reformulated as DEs. Still, one hopes that one does not have to resort to such trickery in future versions of Mathematica... – J. M. is in limbo Oct 5 '12 at 13:30
• @J.M. SolveForUniverseHamiltonian[ currentUniverse ]` – Dr. belisarius Oct 5 '12 at 13:32