# Numerical solution of an ordinary integro-differential equation

I am trying to solve the following integro-differential equation: $$a(\varepsilon) f'(\varepsilon) + b(\varepsilon) f(\varepsilon) + \int_{\varepsilon}^{\varepsilon + \varepsilon_1} R_1(\epsilon) f(\epsilon) d\epsilon - \int_{\varepsilon - \varepsilon_2}^{\varepsilon} R_2(\epsilon) f(\epsilon) d\epsilon = 0,$$ where $$\varepsilon_1 = 10$$, $$\varepsilon_2 = 20$$.

Note that the reduction of this equation into a second-order ODE is possible, but is much harder to solve. See, e.g., this question.

Taking

eps1 = 10;
eps2 = 20;
epsmin = .025;
epsmax = 150;
a[eps_] := 1 + eps
b[eps_] := eps^1.4
R1[eps_] := Piecewise[{{eps, eps1 <= eps <= epsmax}}, 0]
R2[eps_] := Piecewise[{{eps, eps2 <= eps <= epsmax}}, 0]


We know that f[eps] > 0 for epsmin <= eps <= epsmax, and that f[eps] = 0 for eps <= epsmin and eps >= epsmax. For eps >= epsmax, all derivatives of f[eps] also become 0. Also,

Integrate[Sqrt[eps] f[eps], {eps, epsmin, epsmax}] == 1;


I have tried to use NDSolve with NIntegrate as follows:

solBW = NDSolve[{a[eps]f'[eps] + b[eps]f[eps] + NIntegrate[R1[y]f[y],{y,eps,eps + eps1}] - NIntegrate[R2[y]f[y],{y, eps - eps2, eps}] == 0, f[epsmax] == 0}, f, {eps, epsmax, epsmin}];


but it returns

NIntegrate: y = eps is not a valid limit of integration.


Any hint on how to solve this equation numerically is appreciated.

• Differentiate your equation and solve for f'. Subsequently integrate f'. May 24, 2022 at 7:47
• The derivative of the integral terms will produce f[eps + eps1] and f[eps - eps2] terms. I can't proceed from that point. May 24, 2022 at 7:58
• Why not to use collocation method? May 24, 2022 at 8:00
• The term in your equation 1/a[eps] Integrate[R2[eps1] f[eps1], {eps1, eps - eps2, eps}] is in discordance with eps1 = 10; for Mathematica syntax. Upgrade your math culture. May 24, 2022 at 8:27
• @AsaturKhurshudyan It looks like it can be solved with FDM - see my answer. May 24, 2022 at 23:55

It looks like we can solve this problem with using FDM and LinearSolve[] as follows

Clear["Global*"]

y1 = 10;
y2 = 20;
ymin = .025;
ymax = 150.025;
a[y_] := 1 + y;
b[y_] := y^1.4
R1[y_] := Piecewise[{{y, y1 <= y <= ymax}}, 0];
R2[y_] := Piecewise[{{y, y2 <= y <= ymax}}, 0];

h = 1/5; xcol = Range[ymin, ymax, h]; nn = Length[xcol]; u =
Array[f, {nn}];

eq = Join[
Table[a[xcol[[j]]] 1/2 (f[j + 1] - f[j - 1])/h +
b[xcol[[j]]] f[j] +
h Sum[If[k <= Length[xcol], R1[xcol[[k]]] f[k], 0], {k, j,
j + Round[y1/h]}] -
h Sum[If[k >= 1, R2[xcol[[k]]] f[k], 0], {k, j - Round[y2/h],
j}] == 0, {j, 2, nn - 1}], {h Sqrt[xcol] . u == 1,
f[nn] == 0}];
{vec, mat} = CoefficientArrays[eq, u];

sol = LinearSolve[mat, -vec];


Visualization

{ListLinePlot[Table[{xcol[[i]], sol[[i]]}, {i, nn}],
PlotRange -> All],
ListLinePlot[Table[{xcol[[i]], sol[[i]]}, {i, Round[20/h]}],
PlotRange -> All],
ListLinePlot[Table[{xcol[[i]], sol[[i]]}, {i, Round[20/h] + 1, nn}],
PlotRange -> All]}


Now if we turn step to h=1/10, then we have same picture, therefore numerical method is stable

• Sorry for being picky, but why the tail does not decrease to 10^-40 or lower? Both integral terms decrease pretty fast and so should f. May 25, 2022 at 1:18
• May be you need to increase WorkingPrecision up to 40 and decrease h to 1/30, also use rational inputs like ymin = 25/1000. May 25, 2022 at 1:39
• I used \$PreRead = (#/.s_String/;StringMatchQ[s,NumberString]&&Precision@ToExpression@s == MachinePrecision :>s<>" 150."&) ymin =Rationalize[ .025,0] (same for ymax) and h = 1/30`. No difference. May 25, 2022 at 7:01
• @AsaturKhurshudyan We can use dense grid with FEM, but dense grid is not the only method to avoid oscillations. Also oscillations may by exited due to model instability and not as numerical solution instability. May 26, 2022 at 3:08
• @AsaturKhurshudyan This question seems to be answered. Could you upload code with 10 terms in a new topic? May 26, 2022 at 4:50