# 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. Commented Oct 5, 2012 at 10:46
• Would you care to comment about the field in which you've found this problem? (seems pretty difficult!) Commented Oct 5, 2012 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... Commented Oct 5, 2012 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. Commented Oct 5, 2012 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. Commented Oct 5, 2012 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
Commented Nov 30, 2012 at 16:59
• @plato Please take a look at this one mathematica.stackexchange.com/q/58149/193 Commented Aug 25, 2014 at 17:14
• @belisarius Done! Thx for notifying though..:) Commented Aug 25, 2014 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? Commented Jan 6, 2019 at 13:49
• @PlatoManiac, thanks for this excellent answer. Have you considered solving for f(y)?
– Dean
Commented Aug 14, 2022 at 19:57

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[0] == 1, f, x]
(*
->{{f -> Function[{x}, E^(x/2 - 1/4 Sin[2 x])]}}
*)


Testing it

Integrate[g[u] (f /. fs[[1]])[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... Commented Oct 5, 2012 at 13:30
• @J.M. SolveForUniverseHamiltonian[ currentUniverse ]` Commented Oct 5, 2012 at 13:32