Numerically solving a non-linear integro-differential equation

I am interested in solving the non-linear integro-differential equation $$\frac{\mathrm{d} y}{\mathrm{d} \tau} = - \int_{0}^{\tau} \mathrm{d} s \, \big[ y(\tau -s) \big]^{2} y(s) ; \quad y (0) = 1 .$$ I believe that there is no known analytical closed form solution(?) So I wanted to rely on numerical methods, and wrote

eqn = y'[tau] == - Integrate[y[tau-s] y[tau-s] y[s],{s,0,tau}];
init = y[0] == 1;
sol = NDSolveValue[{eqn,init},y,{tau,0,10}];


From the error messages, it seems that NDSolveValue cannot handle such an integro-differential equation. Would there be any easy work-arounds to solve this equation numerically using Mathematica?

• This ndsolve-integro-differential-equation is similar question. unfortunately the nice trick shown there by Carl Woll does not work here since your integrand has $y$ that depends on both delayed input and also not delayed input. If your $y$ inside the integrand depended only on delayed input, then Carl's trick would have worked by replacing $t=\tau-s$ as shown in the above. So my guess is that this can not be solved by NDSolve as it stands but I could be wrong. Commented May 12, 2023 at 16:39
• With MellinTransform I have only a general solution:y[tau] = -E^(I tau) C[1] - E^(-I tau) C[2]. Commented May 12, 2023 at 16:48
• @MariuszIwaniuk what is C[1] and C[2] in the above? constant of integrations? But this is first order ode, should the solution not only have one constant of integration? Commented May 12, 2023 at 16:50
• @Nasser I used example from help pages of MellinTransform. Try:eqn = y'[tau] == -Integrate[y[tau - s]^2*y[s], {s, 0, tau}]; F = MellinTransform[eqn, tau, q] /. y[0] -> 1; R = RSolveValue[F, MellinTransform[y[tau], tau, q], q]; InverseMellinTransform[R, q, tau] // Expand Commented May 12, 2023 at 17:03
• @MariuszIwaniuk: Unfortunately, eqn /. {y[tau] -> E^(I tau), y'[tau] -> D[E^(I tau), tau], y[s] -> E^(I s), y[tau - s] -> E^(I (tau - s))} results in I E^(I tau) == I E^(I tau) (-1 + E^(I tau)) which is not true. Hope I am not mistaken. Commented May 12, 2023 at 18:53

Here's a very inefficient way to solve it (because of re-interpolation and recalculations in the NIntegrate). I'm no expert on differential equation solving, but thought it would be worth sharing a very naïve approach to this because it's surprisingly brief and requires little expert knowledge:

next[pts_, d_] :=
With[{i = Interpolation[Normal@pts, InterpolationOrder -> 1],
p = pts["Part", -1]},
With[{x0 = p[[1]], y0 = p[[2]]},
With[{g = -NIntegrate[i[x0 - s]^2 i[s], {s, 0, x0}]},
pts["Append", {x0 + d, y0 + d*g}]
]]]

points = CreateDataStructure["DynamicArray"];
points["Append", {0, 1}];
With[{n = 100},
result = Nest[next[#, 10./n] &, points, n];
ListLinePlot[Normal@result, PlotRange -> {{0, 10}, {-0.2, 1.1}},
Frame -> True]]


• This is very nice code (+1). Commented May 13, 2023 at 15:20
• I can recommend option InterpolationOrder -> 0 and n=200 then it consider with my solution. Commented May 13, 2023 at 16:02

This equation can be solved numerically using method described in our paper. First note, that integral $$\int_{0}^{\tau} \mathrm{d} s \, \big[ y(\tau -s) \big]^{2} y(s)$$ can be mapped on (-1,1) by substitution $$s=\tau/2 (1+z)$$, therefore

int[tau_] :=
tau/2 Integrate[y[tau/2 (1 - z)]^2 y[tau/2 (1 + z)], {z, -1, 1}];


In turn int can be integrated using Gauss quadrature, for this we call

Get["NumericalDifferentialEquationAnalysis"];


To convert integrodifferential equation into system of algebraic equations we use the Euler wavelets as follows

UE[m_, t_] := EulerE[m, t];
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2)  UE[m, 2^k t - 2 n + 1], (n - 1)/2^(k - 1) <= t <
n/2^(k - 1)}, {0, True}}];
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 4; M0 = 7; With[{k = k0, M = M0},
nn = Length[Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dx = 1/(nn); xl = Table[l*dx, {l, 0, nn}]; xcol =
Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Psi[x_] := Psijk /. t1 -> x;
int1[x_] := Int1 /. t1 -> x; var = Array[a, {nn}];
y[t_] := var . int1[t] + a0 ; dy[t_] := var . Psi[t];

np = nn; g = GaussianQuadratureWeights[np, -1, 1]; points =
g[[All, 1]];
weights = g[[All, 2]];
Int[ff_, z_] :=
Sum[(ff /. z -> points[[i]])*weights[[i]], {i, 1, np}];
intNum[tau_] := tau/2 Int[y[tau/2 (1 - z)]^2 y[tau/2 (1 + z)], z];
T=10;eqs = Table[dy[t]/T + T intNum[t] == 0, {t, xcol}];


To solve system eqs with initial condition y[0]==1 we use

 eqAll = Join[eqs, {y[0] == 1}]; varAll = Join[{a0}, var]; vv =
Table[{varAll[[i]], 1/10}, {i, Length[varAll]}]; sol =
FindRoot[eqAll, vv];


Visualization

Plot[y[t/T] /. sol, {t, 0, T}, PlotRange -> All,
FrameLabel -> {"\[Tau]", "y"}, Frame -> True]


Update 1 We can compare this solution with that proposed by flinty at n=200 (red points) as follows

next[pts_, d_] :=
With[{i = Interpolation[Normal@pts, InterpolationOrder -> 0],
p = pts["Part", -1]},
With[{x0 = p[[1]], y0 = p[[2]]},
With[{g = -NIntegrate[i[x0 - s]^2 i[s], {s, 0, x0}]},
pts["Append", {x0 + d, y0 + d*g}]]]]

points = CreateDataStructure["DynamicArray"];
points["Append", {0, 1}];
With[{n = 200}, result = Nest[next[#, 10./n] &, points, n];
pl = ListPlot[Normal@result, PlotRange -> {{0, 10}, {-0.2, 1.1}},
Frame -> True, PlotStyle -> Red]]

Show[Plot[y[t/T] /. sol, {t, 0, T}, PlotRange -> All,
FrameLabel -> {"\[Tau]", "y"}, Frame -> True], pl]


Update 2. We also can compare our solution with solution proposed by Ulrich as follows (please, pay attention that his code has been updated with int and initial data)

nl = NestList[
Function[{fa},
Block[{int, ip},
int[tau_?NumericQ] :=
Block[{s},
NIntegrate[
fa[tau - s] fa''[tau - s] fa[s] +
fa'[tau - s] fa'[tau - s] fa[s], {s, 0, tau},
Method -> {{Automatic, "LocalAdaptive"}[[-1]],
"SymbolicProcessing" -> 0}]];
NDSolveValue[{f'''[tau] == - f'[tau] - 2 int[tau], f[0] == 1,
f'[0] == 0, f''[0] == -1}, f, {tau, 0, 10}]]], Exp[-#^2] &, 10];
pl1=Plot[nl[[-1]][t], {t, 0, 10}, PlotRange -> All,
PlotStyle -> {Red, Dashed}]
Show[Plot[y[t/T] /. sol, {t, 0, T}, PlotRange -> All,
FrameLabel -> {"\[Tau]", "y"}, Frame -> True], pl1]


• I tried yours code it dosen't work properly? a0 ? Commented May 13, 2023 at 14:47
• @MariuszIwaniuk Thank you. There is a typo in the code with psi definition. Please try new version. Commented May 13, 2023 at 15:19
• Ok Works fine. Thanks Commented May 13, 2023 at 15:23
• @AlexTrounev, Sir, sorry to bother you here! Could you please have a look at my recent post? It is an interesting problem with a working code. Thank you for the previously invaluable help!
– lxy
Commented May 17, 2023 at 11:33

In addition to the interesting answers from Alex Trounev and flinty here I show an iterative solution procedure of the differentiated integro-differential equation using NestList, NIntegrate and NDSolveValue

with bc y[0]==1, y'[0]==0

nl = NestList[
Function[{fa},
Block[{int, ip },
int[tau_?NumericQ] := Block[{s},
NIntegrate[ fa [tau - s] fa '[tau - s] fa [s], {s, 0, tau },
Method -> {{Automatic, "LocalAdaptive"}[[-1]],"SymbolicProcessing"-> 0} ]];
NDSolveValue[{ f'' [ tau] == -f[tau]  - 2 int[tau],
f[0] == 1,f'[0] == 0}, f, {tau, 0 , 10} ]]
], 1 &, 10];

Plot[nl[[-1]][t], {t, 0, 10} ,PlotRange -> {-.1, 1}, PlotStyle -> {Red, Dashed}]


which fits quite well with the solutions found in the other answers.

Unfortunately at tau~7 we see numerical problems in the iteration process, which I don't understand yet.

• It is a nice attempt (+1). We can compute right solution with this code. See Update 2 to my answer. Commented May 17, 2023 at 14:25
• @AlexTrounev clever idea to differentiate one more time. Thank you! Commented May 17, 2023 at 15:11
• @AlexTrounev clever idea to differentiate one more time. Thank you! Still I'm wondering why my attempt suddenly diverges. Commented May 17, 2023 at 15:26
• It is not wondering that solution with your code diverges. But it is unexplained that my solution with updated code converges. Normally solution with iterative method diverges in a case of nonlinear integrodifferential equations. Commented May 17, 2023 at 16:01

Here is another solution using the trapezoidal rule for the time integral in the rhs, and a second-order predictor-corrector for the differential equation:

(*Computes Integrate[-f[\[Tau]-s]^2 f[\[Tau]],{\[Tau],0,T}] using the \
trapezoidal rule.*)
computeIntegral[] := Module[{n, int},
int = dt (0.5  tab[[-1]]^2 tab[[1]] + 0.5 tab[[-1]] tab[[1]]^2);
int += dt Sum[tab[[-1 - i + 1]]^2 tab[[i]], {i, 2, Length@tab - 1}];
-int
];
(*Time-integration using a second-order predictor-corrector.*)
update[] := Module[{int, pred, intnew, corr},
int = computeIntegral[];
pred = tab[[-1]] + dt int;
AppendTo[tab, pred];
intnew = computeIntegral[];
corr = tab[[-2]] + 0.5 dt (int + intnew);
tab[[-1]] = corr;
];


The solution is then obtained via

dt = 0.1;
tab = {1.0};
Do[
update[];
, {i, 1, 100}];
(*.*)
tabt = Table[i dt, {i, 0, 100}];
ListPlot[{tabt, tab} // Transpose, PlotRange -> Full]
`