Non-linear integro-differential equation with infinite domain of integration

Question

$$I(t)=\frac{dy(t)}{dt}+\int_{-\infty}^t dt'\left[\sin\left(\frac{y(t)-y(t')}{2}\right)f(t-t')-\sin\left(\frac{y(t)+y(t')}{2}\right)g(t-t')\right]$$ with $y(t\leq 0)=0$ and $I(t)=I_0\Theta(t)$ with $\Theta(t)$ being the Heaviside step function. $f$ and $g$ are known, well-behaved functions, which we can take as $f(\tau)=g(\tau)=\frac{\tau}{(1+\tau^2)^2}$.

I am currently solving it self-consistently in Matlab starting with a guess for $y(t)$, and using the derivative term to do a finite difference update of $y(t)$. That is, $$y_{new}(t+\Delta t)=\Delta t\left(I(t)-\int_{-\infty}^t dt'\left[\sin\left(\frac{y_{old}(t)-y_{old}(t')}{2}\right)f(t-t')-\sin\left(\frac{y_{old}(t)+y_{old}(t')}{2}\right)g(t-t')\right]\right)$$ I was wondering if there is possibly a neater way to solve this on Mathematica.

Answer 1

The integro-differential equation also can be solved using the Delay Differential Equation feature of NDSolve. This involves discretizing the integral, for instance by

tmax = 50;
fi = FunctionInterpolation[t^2/(1 + t^2)/2, {t, 0, tmax}, MaxRecursion -> 10];

where t^2/(1 + t^2)/2 is the integral of f, determined in my earlier answer. The resulting InterpolationFunction contains values of the integral (84 in total) at

cs = Flatten@fi["Coordinates"]
(* {0., 0.000529258, 0.00211033, ... 41.6667, 45.8333, 50.} *)

with the values themselves,

fi["ValuesOnGrid"]
(* {0., 0.000529258, 0.00211033. ... 0.499712, 0.499762, 0.4998} *)

from which weights can be computed,

ws = -MovingAverage[Subtract @@@ Partition[Join[{%[[1]]}, %, {%[[-1]]}], 2, 1], 2]
(* {0.000264629, 0.00105517, 0.00209703, ... 0.0000586998, 0.0000439571, 0.0000189916} *)

Not surprisingly, Total@ws equals 0.4998 to high accuracy. There are, of course, several ways to discretize the integral, but this one is straightforward and efficient. Applying this discretization to the differential equation given near the beginning of my earlier answer then yields after a small amount of algebra,

eq = D[y[t], t] == i0 + 2 Cos[y[t]/2] 
    Sum[ws[[i]] Sin[y[t - cs[[i]]]/2], {i, Length[cs]}];

Solving eq with i0 = 1 then yields the same result obtain in my earlier answer in a few seconds, rather than the 4480 seconds needed previously. In fact, the ten solutions for i0 = Range[.1, 1., .1] can be obtained in 19 seconds.

Do[i0 = j/10; 
  eq = D[y[t], t] == i0 + 2 Cos[y[t]/2] 
    Sum[ws[[i]] Sin[y[t - cs[[i]]]/2], {i, Length[cs]}];
  s[j] = NDSolveValue[{eq, y[t /; t <= 0] == 0}, y[t], {t, 0, tmax}], {j, 10}]

Plot[Evaluate@Array[s, 10], {t, 0, tmax}, ImageSize -> Large, 
  AxesLabel -> {t, y}, LabelStyle -> {15, Bold, Black}, 
  PlotRange -> All, PlotLegends -> Placed[Range[.1, 1., .1], {.1, .6}]]

Answer 2

The integro-differential equation can be solved with only a modest amount of code, as follows: First note that the integrand can be simplified because, as stated in the question, f[r] == g[r].

int = (Sin[(y[t] - y[tp])/2] - Sin[(y[t] + y[tp])/2]) f[t - tp];
Simplify[int]
(* -2 Cos[y[t]/2] f[t - tp] Sin[y[tp]/2] *)

and -2 Cos[y[t]/2] can be moved from within the integral. Additionally, because y[tp] vanishes for tp <= 0, the lower bound on the integral can be increased from -Infinity to 0. Consequently, the integro-differential equation reduces to

D[y[t], t] == i0 + 2 Cos[y[t]/2] 
  NIntegrate[Sin[y[tp]/2] (t - tp)/(1 + (t - tp)^2)^2, {tp, 0, t}]

Unfortunately, NDSolve is unable to handle such equations. It turns out, though, that the equation can be solved iteratively. For instance, with i0 = 1,

Clear[s, h];
n = 14; tmax = 50; i0 = 1;
s[0] = FunctionInterpolation[i0 t, {t, 0, tmax}];
Do[
  h[t_?NumericQ] := NIntegrate[Sin[s[i - 1][tp]/2] 
    (t - tp)/(1 + (t - tp)^2)^2, {tp, 0, t}]; 
  s[i] = NDSolveValue[{D[y[t], t] == i0 + 2 Cos[y[t]/2] h[t], y[0] == 0}, 
    y, {t, 0, tmax}], {i, 1, n}];
Plot[Evaluate@Array[s[#][t] &, n + 1, 0], {t, 0, tmax}, ImageSize -> Large, 
  AxesLabel -> {t, y}, LabelStyle -> {15, Bold, Black}]

Successive iterations lead to progressively lower curves, with good convergence obtained with ten iterations. Values at t = tmax illustrate the convergence.

Array[s[#][tmax] &, n + 1, 0]
(* 50., 48.5208, 47.4234, 45.481, 43.3451, 42.2148, 41.175, 39.5509, 
   37.8209, 36.9781, 36.6383, 36.4922, 36.4351, 36.4158, 36.41} *)

The computation takes about an hour.

Addendum: Asymptotically constant solutions.

Suppose that y[t] is approximately constant for large t. Then, Sin[y[t]/2] can be removed from within the integral, leaving

Integrate[(t - tp)/(1 + (t - tp)^2)^2, {tp, 0, t}, Assumptions -> t > 0]
{* t^2/(2 + 2 t^2) *}

which reduces to 1/2 for large t, and the integro-differential equation collapses to i0 == -Sin[y/2]/2. Hence, asymptotically flat solutions exist if and only if -1/2 < i0 < 1/2, and the asymptotic values should be, for i0 = {1/10, 2/10, 3/10, 4/10, 5/10},

N@Table[Min@SolveValues[{Sin[x]/2 == -j, Pi < x < 2 Pi}, x], {j, 1/10, 5/10, 1/10}]
(* {3.34295, 3.55311, 3.78509, 4.06889, 4.71239} *)

The earlier code block gives corresponding converged solutions at tmax = 50 of

(* {3.34306, 3.55333, 3.78547, 4.06955, 4.61028} *)

Evidently, the i0 = 1/2 solution is not yet in the asymptotic regime at tmax = 50. Corresponding converged solutions are