Solving integro-differential equation - numerically

Question

I want to solve a integro-differential equation numerically.

The equation is given by :

$\dot{c}(t)=-\int_0^t \mathrm{d}t_1f(t-t_1)c(t_1)$

Hereby, $f(t-t_1)$ will be given a realisation of some random numbers, e.g. $f(t-t_1)$ originally was a rondom variable, and I want to solve the integer differential equation for one realisation. So the function will most certainly be a messy polynomial, where the solution cannot be obtained easily with the Laplace-transform used to solve the equation.

How is the best procedure to numerically solve this (or any) integer differential equation?

Thanks already Martin

Answer 1

To solve numerically this problem is not more difficult than solving a linear system. Using the collocation method, using a family of interpolation functions $\phi_k(t), k = 1,\cdots, n$ over a fixed grid $t_1,\cdots,t_n$ we can postulate $c(t) = \sum_{k=1}^n\alpha_k\phi_k(t)$. With those interpolation functions we can also model $f(t) = \sum_{k=1}^nf_k\phi_k(t)$ and then the problem reduces to determine $\alpha_k$ in

$$ \sum_{k=1}^n\alpha_k\phi'_k(t)+\int_0^t\left(\sum_{k=1}^nf_k\phi_k(t-\tau)\right)\left(\sum_{j=1}^n\alpha_j\phi_j(\tau)\right)d\tau = 0 $$

over the grid $t_1,\cdots, t_n$ so we have to solve the linear system

$$ \sum_{k=1}^n\alpha_k\phi'_k(t_{\nu})+\int_0^{t_{\nu}}\left(\sum_{k=1}^nf_k\phi_k(t_{\nu}-\tau)\right)\left(\sum_{j=1}^n\alpha_j\phi_j(\tau)\right)d\tau = 0,\ \ \ \text{for}\ \ \nu = 1,\cdots, n $$

NOTE

We need non null initial conditions on $c(t)$ to avoid the null solution.

Follows a MATHEMATICA script using as interpolation functions Lagrange polynomials. In red $f(t)$ and in blue $c(t)$.

Clear["Global`*"]
lagrangeInterpolation[xvalues_, yvalues_, var_] := Module[{subs, complements, nominators, denominators}, 
  subs = Subsets[xvalues, {Length@xvalues - 1}];
  complements = Flatten[Complement[xvalues, #] & /@ subs];
  nominators = Times @@ (# - var) & /@ subs;
  denominators = Table[Times @@ (subs[[i]] - complements[[i]]), {i, 1, Length@subs}];
Reverse[nominators/denominators].yvalues]

n = 10;
tf = 5;
c0 = 1;
tvalues = Range[0, tf, tf/(n - 1)];
yvalues = Table[Sin[2 Pi/(n-1) j], {j, 0, n-1}] // N;
ftmtau = Expand[lagrangeInterpolation[tvalues, yvalues, t - tau]];
cs = Table[c[j], {j, 1, n}];
ct = Expand[lagrangeInterpolation[tvalues, cs, t]];
dct = D[ct, t];
ctau = Expand[lagrangeInterpolation[tvalues, cs, tau]];
intt = Integrate[ctau ftmtau, {tau, 0, t}];
res = dct + intt;
equsc1 = (ct /. {t -> 0}) - c0;
equs = Table[res /. {t -> tvalues[[k]]}, {k, 2, n}];
equstot = Join[{equsc1}, equs];

solred = Solve[equstot == 0, cs][[1]];
ctau0 = ctau /. solred
ft0 = Expand[lagrangeInterpolation[tvalues, yvalues, t]];
gr1 = Plot[ft0, {t, 0, tf}, PlotStyle -> Red];
gr2 = Plot[ctau0, {tau, 0, tf}, PlotStyle -> Blue];
Show[gr1, gr2, PlotRange -> All]

The theoretical result can be obtained by solving

$$ s C(s) + C(s)F(s) = 1,\ \ F(s) = \frac{2 \pi }{5 \left(s^2+\frac{4 \pi ^2}{25}\right)} $$

giving

$$ C(s) = \frac{25 s^2+4 \pi ^2}{25 s^3+4 \pi ^2 s+10 \pi } $$

Follows the Laplace inverse plot in blue and $f(t)$ in red