# Solving Integral Equation -numerical solution

I am trying to find a (numerical) solution to the following integral equation:

where $$\epsilon$$ is a real valued function and $$\beta$$ and $$c_0$$ are constants.

MMA code:

-(β*Integrate[ϵ[s]/E^((-s + t)*β), {s, 0, t}]) + β*Integrate[ϵ[s]^2/E^((-s + t)*β), {s, 0, t}] +
β*Integrate[(ϵ[s]^2*ϵ[t])/E^((-s + t)*β), {s,
0, t}] + ϵ[t] - ϵ[t]^2 == c0

• What have you tried? Please post the code for the expression in the image. Jan 6, 2021 at 6:20
• What is the value of c0? Or are you looking for an analytical solution? Jan 6, 2021 at 11:37
• One point I don't understand: The integral equation holds for every t, especially t==0 . Because \[Epsilon][0]==0 it follows c0==0! If yes, the fixpointiteration yields \[Epsilon][t]==0 ! Jan 6, 2021 at 12:45
• @ Ulrich, My bad, you are right. [Epsilon][0] is not zero. I removed that from the post. c_0 is a positive number. Thank you! Jan 6, 2021 at 17:43
• @Klaas-Jan Thanks for correction. I'm still waiting and hoping for a numerical answer to your question. 'til now only "@AlexTrounev" tried a numerical answer Jan 11, 2021 at 8:57

You can transform the integral equation into a second order differential equation and solve with NDSolve.

Multiply the equation with E^(beta*t) and bring it to one side to get feq. Differentiate two times and eliminate the remaining integral. Get initial conditions.

feq = (-c0)*E^(t*\[Beta]) - \[Beta]*
Integrate[E^(s*\[Beta])*\[Epsilon][s],
{s, 0, t}] + \[Beta]*
Integrate[E^(s*\[Beta])*\[Epsilon][s]^2, {s, 0, t}] +
\[Beta]*\[Epsilon][t]*
Integrate[E^(s*\[Beta])*\[Epsilon][s]^2, {s, 0, t}] +
E^(t*\[Beta])*\[Epsilon][t] - E^(t*\[Beta])*\[Epsilon][t]^2;

dfeq = D[feq, t]

d2feq = D[dfeq, t]

eli[t_, \[Beta]_, c0_] = Simplify[Eliminate[{dfeq == 0, d2feq == 0},
Integrate[E^(s*\[Beta])*\[Epsilon][s]^2, {s, 0, t}]]]

{{ic1a[c0_]}, {ic1b[c0_]}} = Solve[feq == 0 /. t -> 0, \[Epsilon][0]]

ic2a[\[Beta]_, c0_] =
Solve[dfeq == 0 /. t -> 0 /. ic1a[c0], \[Epsilon]'[0]][[1, 1]]

ic2b[\[Beta]_, c0_] =
Solve[dfeq == 0 /. t -> 0 /. ic1b[c0], \[Epsilon]'[0]][[1, 1]]

\[Epsilon]sola[\[Beta]_, c0_] := \[Epsilon] /.
First@NDSolve[{eli[t, \[Beta], c0], Equal @@ ic1a[c0],
Equal @@ ic2a[\[Beta], c0]}, \[Epsilon], {t, 0, 20}]

\[Epsilon]solb[\[Beta]_, c0_] := \[Epsilon] /.
First@NDSolve[{eli[t, \[Beta], c0], Equal @@ ic1b[c0],
Equal @@ ic2b[\[Beta], c0]}, \[Epsilon], {t, 0, 20}]

Manipulate[{Plot[
Evaluate[(aa = \[Epsilon]sola[\[Beta], c0])[t]], {t,
0, .99*aa[[1, 1, 2]]}, PlotRange -> All, GridLines -> Automatic,
ImageSize -> 300] // Quiet,
Plot[Evaluate[(aa = \[Epsilon]solb[\[Beta], c0])[t]], {t,
0, .99*aa[[1, 1, 2]]}, PlotRange -> All, GridLines -> Automatic,
ImageSize -> 300] // Quiet}, {{\[Beta], .1}, -3, 3,
Appearance -> "Labeled"}, {{c0, -.5}, -3, 1/4,
Appearance -> "Labeled"}]


Test: Show that one sided equation is very close to zero. Takes a few minutes. (Even better with higher WorkingPrecision of NDSolve and NIntegrate).

ff[t_?NumericQ, \[Beta]_, c0_] =
feq /. \[Epsilon] -> (aa = \[Epsilon]sola[.1, -.5]); pl =
Plot[ff[t, .1, -.5] /. Integrate -> NIntegrate, {t, 0,
aa[[1, 1, 2]]}]


• I checked your solution with \[Beta] -> 1, c0 -> .2 but NDSolve gives a warning message "NDSolveValue::ndsz: At t == 0.3874777341671869, step size is effectively zero; singularity or stiff system suspected." and stops evaluation. Any idea why? Thanks! Jan 14, 2021 at 7:45
• Yes, i think for certain parameters the integral equation does not hold for larger t values. See also my comment to your answer. You may get the same termination at higher t for certain beta, c0 with my proposed change of nested equation. Jan 14, 2021 at 19:19

An iterative approach based on the fixed point existence.

$$\epsilon_{k+1}(t)=\Phi\left(\epsilon_{k}(t),t\right)$$

f = c0;
sols = {f}
n = 9;
For[k = 1, k <= n, k++,
Clear[ϵ];
ϵ[t_] := f;
f = (β*Integrate[ϵ[s]/E^((t - s)*β), {s, 0, t}]) - β*Integrate[ϵ[s]^2/E^((t - s)*β), {s, 0, t}] - β*ϵ[t]*Integrate[(ϵ[s]^2)/E^((t- s)*β), {s, 0, t}] + ϵ[t]^2 + c0;
AppendTo[sols, f]
]

parms = {c0 -> 0.2, β -> 1};
funcs = sols /. parms;
Plot[funcs, {t, 0, 5}]


As a checking the convergence

ϵ[t_] := funcs[[n+1]];
dif = (β*Integrate[ϵ[s]/E^((t - s)*β), {s, 0, t}]) - β*Integrate[ϵ[s]^2/E^((t - s)*β), {s, 0, t}] - β*ϵ[t]*Integrate[(ϵ[s]^2)/E^((t- s)*β), {s, 0, t}] + ϵ[t]^2 + c0 - ϵ[t];

Plot[dif, {t, 0, 5}, PlotStyle -> {Thick, Blue}, PlotRange -> All]


And now follows a layman solution (second order splines) to this problem.

fl[a_, b_, c_, delta_, t_, n_] := Sum[(UnitStep[t - k delta] - UnitStep[t - (k + 1) delta]) (a[k] + b[k] (t - k delta) + c[k] (t - k delta)^2), {k, 0, n}]

n = 40;
tmax = 10;
delta = tmax/n;
fl0 = fl[a, b, c, delta, t, n];
For[k = n - 1, k >= 0, k--,
fl0 = fl0 /. {a[k + 1] -> a[k] + delta b[k] + delta^2 c[k], b[k + 1] -> b[k] + 2 delta c[k]}]
vars = Join[{a[0], b[0]}, Table[c[k], {k, 0, n}]];

Clear[ϵ]
ϵ[t_] := fl0;
dif = (β*Integrate[ϵ[s]/E^((t - s)*β), {s, 0, t}]) - β*Integrate[ϵ[s]^2/E^((t - s)*β), {s, 0, t}] - β*ϵ[t]*Integrate[(ϵ[s]^2)/E^((t- s)*β), {s, 0, t}] - ϵ[t] + ϵ[t]^2 + c0;

parms = {c0 -> 0.20, β -> 0.2};
dif0 = dif /. parms;
points = Table[dif0, {t, 0, tmax, delta/2}];
npts = Length[points];
diag = Table[(npts - k + 1), {k, 1, npts}];
obj = points.DiagonalMatrix[diag].points;
sol = NMinimize[obj, vars, Method -> "DifferentialEvolution"];
sol[[1]]
et = fl0 /. sol[[2]];
Plot[et, {t, 0, tmax}, PlotRange -> All, PlotStyle -> {Thick, Blue}]


NOTE

Also at $$t=0$$ we have

ϵ[0]^2 - ϵ[0] + c0 == 0


so for $$c_0\lt \frac 14$$ two $$\epsilon[0]$$ are feasible, suggesting that two branches are allowed.

• Can we done with NIntegrate not Integrate ,because is very slow? Jan 6, 2021 at 12:39
• With Integrate we obtain a parametric solution. Jan 6, 2021 at 12:41
• try FixedPointList Jan 6, 2021 at 12:47
• @Cesareo It is nice code (+1). But you did not pay attention there are two branch of solution. Jan 6, 2021 at 15:40
• @Cesaro. Thanks, the last term on the rhs of "f" function. It should be +c_0 not -c_0, right? Jan 6, 2021 at 17:48

We can solve this equation using colocation method and Haar wavelets. There are two branch of solution for every $$c0 < 1/2$$. For some $$\beta, c0$$ we can get these branches with code

Get["NumericalDifferentialEquationAnalysis"];
plot[cc_, bb_, rr_, nn_] :=
Module[{c0 = cc, b = bb, r = rr, n = nn},
M = Sum[1, {j, 0, n, 1}, {i, 0, 2^j - 1, 1}] + 1;
dx = 1/M; A = 0; xl = Table[A + l*dx, {l, 0, M}];
xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, M + 1}];
psi1[x_] := WaveletPsi[HaarWavelet[], x];
psi2[x_] := WaveletPhi[HaarWavelet[], x];
psi1jk[x_, j_, k_] := psi1[j*x - k];
psi2jk[x_, j_, k_] := psi2[j*x - k];
psijk[x_, j_, k_] := Sqrt[j]*(psi1jk[x, j, k] + psi2jk[x, j, k]);
np = M; gw = GaussianQuadratureWeights[np, -1, 1];
points = gw[[All, 1]]; weights = gw[[All, 2]];
GaussInt[ff_, z_] :=
Sum[(ff /. z -> points[[i]])*weights[[i]], {i, 1, np}];
u[t_] := Sum[
a[j, k]*psijk[t, 2^j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}] +
a0 + 1/2 + r Sqrt[1/4 - c0];
u[t_] := Sum[
a[j, k]*psijk[t, 2^j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}] +
a0 + 1/2 + r Sqrt[1/4 - c0];
varM = Join[{a0},
Flatten[Table[a[j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]];
int1[t_] := (t/
2) GaussInt[(1 + z) u[t/2 (z + 1)] Exp[((-t/2 (1 + z) + t)*b)],
z](*s\[Rule]x/2 (1+z)*);
int2[t_] := (t/
2) GaussInt[(1 + z) u[
t/2 (z + 1)]^2 Exp[((-t/2 (1 + z) + t)*b)], z];
varM = Join[{a0},
Flatten[Table[a[j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]];
int1[t_] := (t/
2) GaussInt[(1 + z) u[t/2 (z + 1)] Exp[((-t/2 (1 + z) + t)*b)],
z](*s\[Rule]x/2 (1+z)*);
int2[t_] := (t/
2) GaussInt[(1 + z) u[
t/2 (z + 1)]^2 Exp[((-t/2 (1 + z) + t)*b)], z];
eq = Table[
u[t] - u[t]^2 - b int1[t] + b int2[t] + b u[t] int2[t] - c0 ==
0, {t, xcol}];
sol = FindRoot[eq, Table[{varM[[i]], r/10}, {i, Length[varM]}],
MaxIterations -> 1000];
unum = Table[ {xcol[[i]], Evaluate[u[xcol[[i]]] /. sol]}, {i,
Length[xcol]}]; ListLinePlot[unum]]


Here $$u=\epsilon, b=\beta$$, $$n$$ is number of colocation points and $$r=\pm 1$$ is the indicator of branch. For example, for $$c0=-0.5, \beta =0.1$$ we have

{plot[-.5, 0.1, -1, 4], plot[-.5, 0.1, 1, 4]}


• I tried with: c0 = -0.2; \[Beta] = 1;.The Plot is different than from answer from user @Cesareo ? Jan 6, 2021 at 17:27
• @MariuszIwaniuk First please pay attention that wavelets solution is defined on $0\le t\le 1$. Also this solution not stable with $\beta =1$. You may try to compare my example with Cesareo. Jan 6, 2021 at 22:09

Waiting some time for a straightforward numerical answer, here my attempt, which assumes  \[Beta] = 1 (without loss of generality) and predefined values c0==0.2,T0==5 :

c0 = 2/10; T = 5;
gip[eps_] :=Module[{x, t},
Interpolation[Table[{t,eps[t]^2 + c0 + Exp[-t] NIntegrate[eps[s] Exp[s], {s, 0, t}] -Exp[-t] (1 + eps[t]) NIntegrate[eps[s]^2 Exp[s], {s, 0, t}]}, {t, Subdivide[0, T, 25]}]]]


Function gip get's a pure function as input argument and returns an InterpolationFunction, which might be used iteratively.

With appropriate starting value eps[0]== 1/2 (1 - Sqrt[1 - 4 c0]) (second solution branch eps[0]== 1/2 (1 + Sqrt[1 - 4 c0]) omitted ) it follows

solm = NestList[gip, 1/2 (1 - Sqrt[1 - 4 c0]) &, 15] ;
Plot[Evaluate[Through[solm [t]]], {t, 0, T}, PlotRange -> {0, All}]


The solution doesn't match completely @cesaro's one (shows superimposed oscillations). Perhaps it shows a first step for a completely numerical solution.

To exclude the influence of simple Interpolation I also tried a numerical solution using NDSolveValue (substitution i1[t]==Integrate[eps[s]Exp[s],{s,0,t}] and i2[t]==Integrate[eps[s]^2 Exp[s],{s,0,t}]) :

ff = Function[{t},
NDSolveValue[{i1'[t] == #[t] Exp[t], i1[0] == 0,
i2'[t] == #[t]^2 Exp[t], i2[0] == 0},
#[t]^2 + c0 + Exp[-t] (i1[t] - (1 + #[t]) i2[t])
, {t, 0,T}, DependentVariables -> {i1, i2}
,Method -> {Automatic ,"DiscontinuityProcessing" -> False  },AccuracyGoal -> 10 ]] & ;


which gives the same result(I used only 7 iterations because of increased evaluation time):

sol = NestList[ff, 1/2 (1 - Sqrt[1 - 4 c0]) &, 7]; // AbsoluteTiming
Plot[Evaluate[Through[sol[t]]], {t, 0, T}, PlotRange -> {0, All}]


conclusion

The two numerical solutions, obtained with two independent methods, match very well. I'm quite convinced that the solutions describe the iteration of the given integral equation.

Hints for improvement are welcome!

• Dear @UlrichNeumann, up to now this solutions give the wrong derivative at t ==0 and bad equ test. If you solve the whole equation for eps[t], like sol = Solve[-(\[Beta]* Integrate[\[Epsilon][s]/E^((-s + t)*\[Beta]), {s, 0, t}]) + \[Beta]* Integrate[\[Epsilon][s]^2/E^((-s + t)*\[Beta]), {s, 0, t}] + \[Beta]*\[Epsilon][ t] Integrate[(\[Epsilon][s]^2)/E^((-s + t)*\[Beta]), {s, 0, t}] + \[Epsilon][t] - \[Epsilon][t]^2 == c0, \[Epsilon][t]]  you get two branches. Making nested iteration with that, your code gives very simular solutions to mine. Jan 14, 2021 at 19:15
• @Akku14 Thanks for your reply. I'll check your hint . Because the iterated integraleequation (in my approach) doesn't depend on the derivative eps'[t] I still have doubts. Jan 15, 2021 at 7:17
• @Akku14 I dont't think, that it's allowed to solve the Integral equation for \[Epsilon][t] , thereby ignoring the parts \[Epsilon][s] . This is only allowed for t==0 and gives the two initial conditions depending on c0. In my answer I only showed the solution for the first ic, the second follows similar. Jan 15, 2021 at 7:46
• To @UlrichNeumann You didn't understand me right. I never wanted to solve the int equation this way. Produce gip2 like c0 = -.5; T = 20; \[Beta] = .1; gip2[eps_] :=  Jan 15, 2021 at 7:56
• Module[{x, t}, Interpolation[ Table[{t, (1/ 2)*(1 + \[Beta]*Exp[-t \[Beta]] NIntegrate[E^((s)*\[Beta])* eps[s]^2, {s, 0, t}] - Sqrt[1 - 4*c0 - 4*\[Beta]*Exp[-t \[Beta]] NIntegrate[ E^((s )*\[Beta])*eps[s], {s, 0, t}] + 6*\[Beta]* Exp[-t \[Beta]] NIntegrate[E^((s )*\[Beta])*eps[s]^2, {s, 0, t}] + \[Beta]^2*Exp[-t \[Beta]] NIntegrate[E^((s)*\[Beta])*eps[s]^2, {s, 0, t}]^2])}, {t, 0, T, T/24}]]]  Jan 15, 2021 at 7:58