9
$\begingroup$

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?

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – Nasser
    Commented May 12, 2023 at 16:39
  • $\begingroup$ With MellinTransform I have only a general solution:y[tau] = -E^(I tau) C[1] - E^(-I tau) C[2]. $\endgroup$ Commented May 12, 2023 at 16:48
  • $\begingroup$ @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? $\endgroup$
    – Nasser
    Commented May 12, 2023 at 16:50
  • $\begingroup$ @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 $\endgroup$ Commented May 12, 2023 at 17:03
  • 1
    $\begingroup$ @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. $\endgroup$
    – user64494
    Commented May 12, 2023 at 18:53

4 Answers 4

5
$\begingroup$

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

enter image description here

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

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]

Figure 1

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]

Figure 2

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]

Figure 3

$\endgroup$
4
  • $\begingroup$ I tried yours code it dosen't work properly? a0 ? $\endgroup$ Commented May 13, 2023 at 14:47
  • $\begingroup$ @MariuszIwaniuk Thank you. There is a typo in the code with psi definition. Please try new version. $\endgroup$ Commented May 13, 2023 at 15:19
  • $\begingroup$ Ok Works fine. Thanks $\endgroup$ Commented May 13, 2023 at 15:23
  • $\begingroup$ @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! $\endgroup$
    – lxy
    Commented May 17, 2023 at 11:33
5
$\begingroup$

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

enter image description here

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

enter image description here

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.

$\endgroup$
4
  • $\begingroup$ It is a nice attempt (+1). We can compute right solution with this code. See Update 2 to my answer. $\endgroup$ Commented May 17, 2023 at 14:25
  • $\begingroup$ @AlexTrounev clever idea to differentiate one more time. Thank you! $\endgroup$ Commented May 17, 2023 at 15:11
  • $\begingroup$ @AlexTrounev clever idea to differentiate one more time. Thank you! Still I'm wondering why my attempt suddenly diverges. $\endgroup$ Commented May 17, 2023 at 15:26
  • $\begingroup$ 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. $\endgroup$ Commented May 17, 2023 at 16:01
2
$\begingroup$

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]

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.