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

  • 1
    $\begingroup$ What data do you have? Values for $c(t)$ and $t_1$? $\endgroup$ Sep 8, 2017 at 15:26
  • $\begingroup$ Also, what is the regularity of f? Is it, e.g. Continuous, bounded derivative? $\endgroup$
    – Zach Boyd
    Sep 9, 2017 at 18:13
  • $\begingroup$ Thanks for the reply. I have none. I just have discrete values for $f(t-t_1)$ for each $t_1$ from $0$ up to $t$. In the end I want to make the steps of $t_1$ very small and 'smear out' my discrete values to get a continuous function for $f(t-t_1)$. Then I want to solve the integro-differential equation given. The functional dependence may be very messy, so solving it with the Laplace-transform is not my first choice (I would need the inverse Laplace-transform at some point, and that may be difficult because of many roots of the function). Thus I want to know if there are numberical approaches. $\endgroup$
    – Martin
    Sep 9, 2017 at 18:14
  • $\begingroup$ Why can't we just differentiate both sides and solve the resultingequation using standard RK methods? $\endgroup$
    – K.defaoite
    May 30, 2021 at 21:38

1 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 $$


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)$.

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

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]

enter image description here

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)} $$


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

enter image description here


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .