10
$\begingroup$

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 *)
Improve this question
$\endgroup$
0

3 Answers 3

Reset to default
15
$\begingroup$

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

enter image description here

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

enter image description here

Related post:

How to plot and solve the numerical solution of a integro-differential equation

Improve this answer
$\endgroup$
5
  • $\begingroup$ "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.) $\endgroup$ Commented Jun 3, 2015 at 13:29
  • $\begingroup$ 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. $\endgroup$ Commented Jun 3, 2015 at 13:48
  • $\begingroup$ @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… $\endgroup$ Commented Jun 4, 2015 at 3:51
  • $\begingroup$ 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/… $\endgroup$ Commented Jun 19, 2015 at 7:31
  • $\begingroup$ I meet the same question with you, I also need to solve the question, and i find the answer have something wrong $\endgroup$ Commented Jul 9, 2018 at 18:15
10
$\begingroup$

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

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

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} *)
Improve this answer
$\endgroup$
4
  • $\begingroup$ Very interesting. BTW the dae = y[t] == y0[t]; line isn't necessary. $\endgroup$ Commented Jun 3, 2015 at 12:30
  • $\begingroup$ As to the accuracy part: it's possible to improve the accuracy of FunctionInterpolation, see the edit of my answer. $\endgroup$ Commented Jun 3, 2015 at 12:31
  • 1
    $\begingroup$ @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. $\endgroup$ Commented Jun 3, 2015 at 13:26
  • $\begingroup$ 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 囧 $\endgroup$ Commented Jun 3, 2015 at 13:36
3
$\begingroup$

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

enter image description here

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

Improve this answer
$\endgroup$
5
  • 2
    $\begingroup$ It's faster if you Inactivate[.., Integrate] first and activate it again after DSolve[]. $\endgroup$ Commented Aug 8, 2016 at 21:24
  • 1
    $\begingroup$ A awkward truth is, if one modify the second argument of DSolve i.e. the x[t] to x, DSolve will fail... $\endgroup$ Commented Oct 6, 2016 at 11:44
  • $\begingroup$ I cannot reproduce your result with Mathematica 11.3 on Windows 10, 64 bit.. Suggestions? $\endgroup$ Commented Jul 10, 2018 at 5:43
  • $\begingroup$ @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. $\endgroup$ Commented Jul 10, 2018 at 7:48
  • 1
    $\begingroup$ How embarrassing. I even answered 134445, but I had forgotten about it. Thanks for reminding me. $\endgroup$ Commented Jul 10, 2018 at 13:08

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.