Solving Integral Equation -numerical solution

Question

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

Answer 1

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]]}]

Answer 2

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.

Answer 3

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]}

Answer 4

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!