# How to solve the differential equation with Duhamel's integral?

How do I solve a differential equation with Duhamel's integral? I tried to solve it with NDSolve, but failed:

NDSolve[
{0.01 - 6.25 x[t] + (1.2 Integrate[x'[t - τ]/Sqrt[τ], {τ, 0, t}])/10^7 == 16 x''[t],
x[0] == 0, x'[0] == 0}, x, {t, 0, 10}]
(* The input is returned unevaluated with a warning *)


NDSolve currently can't handle this kind of differential equation, LaplaceTransform is your friend. Since in this case inverse Laplace transform can't be done analytically by InverseLaplaceTransform, you need the help of numerical Laplace inversion package in addition:

eq = {0.01 - 6.25 x[t] + (1.2 Integrate[x'[t - τ]/Sqrt[τ], {τ, 0, t}])/10^7 == 16 x''[t]};
ic = {x[0] == 0, x'[0] == 0};

f[s_] = With[{l = LaplaceTransform},
Solve[l[Rationalize[eq, 0], t, s] /. Rule @@@ ic, l[x[t], t, s]]][[1, 1, -1]]

sol = FunctionInterpolation[GWR[f, t], {t, $MachineEpsilon, 10}] Plot[sol[t], {t, 0, 10}]  Remark: 1. eq is Rationalized because the first argument of GWR should not contain approximate numbers. 2. $MachineEpsilon is chosen as the lower limit of FunctionInterpolation because GWR has some difficulty in calculating the inversion at $t=0$.

The figure of sol already looks nice, but if you're not satisfied with the current accuracy of FunctionInterpolation, you can improve it by adding derivatives of the expression to its first argument. This is mentioned in the Scope of the document of FunctionInterpolation (Why it's not placed in a more conspicuous place?):

gwr[t_?NumericQ] := GWR[f, t]

solOrder3 =
FunctionInterpolation[
D[gwr@t, {t, #}] & /@ Range[0, 3] // Evaluate, {t, $MachineEpsilon, 10}] (* Notice Chop is necessary in this case *) Plot[solOrder3[t] // Chop, {t, 0, 10}] (* Error check *) Plot[Subtract @@@ eq /. x -> solOrder3 // Abs // Evaluate, {t, ##}, PlotPoints -> 20, MaxRecursion -> 2, PlotRange -> All] & @@@ {{0., 0.002}, {0.002, 1}, {1, 10}} // GraphicsRow  Related post: How to plot and solve the numerical solution of a integro-differential equation • "if you're not satisfied with the current accuracy of FunctionInterpolation, you can improve it by adding derivatives of the expression to its first argument" - right, this is why Hermite interpolation is a good idea for approximating analytic functions; the more information you have, (hopefully) the better your approximant is. (This is also probably why Michael's solution was more accurate than your previous one, since the resulting InterpolatingFunction[] has derivative information.) – J. M. is away Jun 3 '15 at 13:29 • You could use the properties of the Laplace transforms of derivatives to create the higher order solOrder3 instead of numeric differentiation. It will be faster. You could also use the ODE to solve for the higher order derivatives, since, at the initial condition, the integration term is zero. – Michael E2 Jun 3 '15 at 13:48 • @MichaelE2 Er…sorry, but I just can't figure out how to calculate the n-th derivative at$t=0$where$n>2$.$x''(0)$can be obtained by Solve[eq /. t -> 0, x''[0]]/. x[0] -> 0 (* => x''[0] == 0.000625 *), but when it comes to D[eq, t], a$\frac{x'(0)}{\sqrt{t}}$term involves in. One may argue that since$x'(0)=0$so this term is (probably) zero at$t=0$, too, but when it comes to D[eq, {t, 2}], a$\frac{x''(0)}{\sqrt{t}}$term involves in… – xzczd Jun 4 '15 at 3:51 • Sorry for bothering you again. I changed my model into a more sophisticated one, and LaplaceTransform method seems cannot work. Is there some method to solve it? This is the question:mathematica.stackexchange.com/questions/86244/… – Scott Wang Jun 19 '15 at 7:31 • I meet the same question with you, I also need to solve the question, and i find the answer have something wrong – Hukai Jul 9 '18 at 18:15 One can use Picard-type iteration to get the solution: Using an approximation to x'[t] (in the integral), we can integrate the ODE to obtain a new approximation. Remarkably, it converges in just two steps. My original thought was to step through the integration using the tools from tutorial/NDSolveStateData to build an interpolation of x'[t] at each step for use in the integration term; that proved too difficult to manage (or perhaps I had set it up in way that made it difficult). The approximation of x'[t] is represented by xp[t]. We start with the initial guess for it to be xp[t] == 0.01 t, which corresponds to extrapolating from the initial conditions (by inspection -- one might solve the ODE for x''). (Actually, starting with xp[t] == 0 works nearly as well and makes the first iteration faster.) We put the integration factor in separate black-box function y0. Adding the dummy algebraic equation y[t] == y0[t] to the system helps with the accuracy. ClearAll[xp, y0, t, x, y]; xp = 0.01 # &; y0[t_?NumericQ] := NIntegrate[xp[t - τ]/Sqrt[τ], {τ, 0, t}]; ode = 0.01 - 6.25 x[t] + 1.2 y0[t] / 10^7 == 16 x''[t]; dae = y[t] == y0[t]; ics = {x[0] == 0, x'[0] == 0}; {sol[10.]} = NDSolve[{ode, ics, dae}, x, {t, 0, 10}]; xp = x' /. sol[10.]; (* iterate with next approximation to x' *) {sol["Final"]} = NDSolve[{ode, ics, dae}, x, {t, 0, 10}];  Let's compare with the solution produced by the numerical Laplace method used by xzczd's answer. In what follows, we'll use sol["Laplace"] = x -> FunctionInterpolation[GWR[f, t], {t,$MachineEpsilon, 10}]


where f and GWR are as in the other answer.

The solutions are roughly the same:

Plot[{x[t] /. sol["Laplace"], x[t] /. sol["Final"]}, {t, 0, 10}]


We can compare how well the solutions track the ODE. The main reason that the Laplace method appears much worse is due to FunctionInterpolation. It does, however, appear to be a better approximation at small values of t. The function opODE gives the residual of a given solution sol at time t of the OP's ODE, with NIntegrate in place of Integrate.

opODE[t_?NumericQ, sol_] := Hold[
0.01 - 6.25 x[t] + (1.2 NIntegrate[x'[t - τ]/Sqrt[τ], {τ, 0, t}])/10^7 - 16 x''[t]
] /. sol // ReleaseHold;

GraphicsRow[
Plot[{opODE[t, sol["Laplace"]], opODE[t, sol["Final"]]},
{t, ##}, PlotPoints -> 20, MaxRecursion -> 2, PlotRange -> All
] & @@@ {{0., 0.001}, {0.001, 1}, {1, 10}}
]


NDSolve is a better alternative to FunctionInterpolation for constructing an accurate interpolation. Oddly the Laplace method shows a similar erratic behavior near t == 0 as the NDSolve-iteration method. The function GWR of the numerical Laplace inversion package needs the argument t to be numeric, but does not protect it with ?NumericQ; hence the wrapper gwr below. With this method of interpolation the numerical Laplace method seems comparable.

gwr[t_?NumericQ] := GWR[f, t];
{sol["Laplace"]} = NDSolve[{x[t] == gwr[t], y'[t] == 1,
y[$MachineEpsilon] ==$MachineEpsilon},
x, {t, $MachineEpsilon, 10}]; GraphicsRow[ Plot[{opODE[t, sol["Laplace"]], opODE[t, sol["Final"]]}, {t, ##}, PlotPoints -> 20, MaxRecursion -> 2, PlotRange -> All ] & @@@ {{0., 0.002}, {0.002, 1}, {1, 10}} ]  Presumably x -> gwr (from @xzczd) produces a highly accurate solution, but it takes a long time to evaluate. For instance, opODE[0.01, x -> gwr] // AbsoluteTiming opODE[9.95, x -> gwr] // AbsoluteTiming (* {27.7618, -1.13159*10^-13 - 7.74942*10^-15 I} {35.6744, -1.36696*10^-15 - 6.18062*10^-26 I} *)  • Very interesting. BTW the dae = y[t] == y0[t]; line isn't necessary. – xzczd Jun 3 '15 at 12:30 • As to the accuracy part: it's possible to improve the accuracy of FunctionInterpolation, see the edit of my answer. – xzczd Jun 3 '15 at 12:31 • @xzczd Yes, that's what I meant at the end of my post, only I didn't take the time to investigate completely. For me, without the dae, the first step taken by NDSolve is about 0.01 (as I recall) with a max. error of about 0.03, so that the solution does not settle down as quickly. That's why I was showing the {0., 0.002} interval, where this method is a little shaky, why I included dae, and why I made it separate code that can be easily taken out of the system. – Michael E2 Jun 3 '15 at 13:26 • Oh, "Adding the dummy algebraic equation y[t] == y0[t] to the system helps with the accuracy." I missed this sentence, seems that I'm a little tired today 囧 – xzczd Jun 3 '15 at 13:36 Answer from Mathematica 11.0.0 for Microsoft Windows (64-bit) (July 28, 2016) ieqn = 0.01 - 6.25*x[t] + 1.2* Integrate[x'[t - s]/Sqrt[s], {s, 0, 10}]/10^7 == 16*x''[t]; ic = {x[0] == 0, x'[0] == 0}; sol = DSolve[{ieqn, ic}, x[t], t] $\{\{x(t)\to 0.0016\, -0.0016 \cos (0.625 t)\}\}\$

Plot[x[t] /. sol, {t, 0, 10}]


MMA can solve integro-differential equations now.

EDITED 10.07.2018:

It seems that in MMA 11.3 DSolve can't solve and return unevaluated.Probably it's a bug(regression).

• It's faster if you Inactivate[.., Integrate] first and activate it again after DSolve[]. – Michael E2 Aug 8 '16 at 21:24
• A awkward truth is, if one modify the second argument of DSolve i.e. the x[t] to x, DSolve will fail... – xzczd Oct 6 '16 at 11:44
• I cannot reproduce your result with Mathematica 11.3 on Windows 10, 64 bit.. Suggestions? – bbgodfrey Jul 10 '18 at 5:43
• @bbgodfrey. Probably it's related to my question: mathematica.stackexchange.com/questions/134445/…. I send e-mail to Wolfram Technical Support and them patched the bug. They did not predict that it would affect to other equations,its regression. – Mariusz Iwaniuk Jul 10 '18 at 7:48
• How embarrassing. I even answered 134445, but I had forgotten about it. Thanks for reminding me. – bbgodfrey Jul 10 '18 at 13:08