I have the following equation I am trying to plot,

$r''=f(r)+\sum_{n=1}^{\infty}a^n \frac{d^n}{dt^n}(f(r))$

and $r$ is a vector in $(x,y,z)$ that are time dependent. I wrote the code as follows,

Clear[x, y, z, t, a];
Dmax = Infinity;
tmax = 10;
r = {x[t], y[t], z[t]};
equation =  D[r, {t, 2}] + r/ (r^2 + 2^2)^(1/2) + 
 Sum[(a^n (D[r/ (r^2 + 2^2)^(1/2),   {t, n}])),{n,Dmax};

 p = Thread[(r /. t -> 0) - {1.0, 0, 0} == 0];
 v = Thread[(D[r, t] /. t -> 0) - {0.0, 1.0, 1.0} == 0];

 Equation2 = Join[Thread[equation == 0], p, v];

 solution1 = 
 NDSolve[Equation2 /. a -> 0.0, r, {t, 0, tmax}, 
 MaxSteps -> Infinity, 
 Method -> {"EquationSimplification" -> "Residual"}];

 solution2 = 
 NDSolve[Equation2 /. a -> 0.1, r, {t, 0, tmax}, 
 MaxSteps -> Infinity, 
 Method -> {"EquationSimplification" -> "Residual"}];

 b1[t_] = r /. solution1;
 b2[t_] = r /. solution2;

  plot3D1 = 
  ParametricPlot3D[{b1[t][[1, 1]], b1[t][[1, 2]], b1[t][[1, 3]]}, {t, 
  0, tmax}, PlotStyle -> Cyan, DisplayFunction -> Identity, 
  PlotLegends -> {"b1"}];

  plot3D2 = 
  ParametricPlot3D[{b2[t][[1, 1]], b2[t][[1, 2]], b2[t][[1, 3]]}, {t, 
  0, tmax}, PlotStyle -> {Blue, Dotted, Thickness -> 0.02}, 
  DisplayFunction -> Identity, PlotLegends -> {"b2"}];

  Show[plot3D1, plot3D2, PlotRange -> All, BoxRatios -> Automatic, 
  ImageSize -> 450, ViewPoint -> Above]

  Sum[(a^n (D[r/ (r^2 + 2^2)^(1/2), {t, n}])),{n,Dmax}

The issue I have is that NDSolve works for $Dmax=Infinity$ as well as for values of $2,3$, however it fails to solve for values equal and greater than $4$. How do I fix the problem? Why it works at infinity?

  • 2
    $\begingroup$ Your code doesn't run for Dmax = Infinity: NDSolve gives error ` NDSolve::ndord: Derivative order K[1] in term (x^(K[1]))[t] should be a non-negative machine-sized integer.` concerning the highest derivative in the ode. $\endgroup$ Aug 4, 2022 at 13:31
  • 1
    $\begingroup$ I'd recommend first getting an idea of what $\sum_{n=1}^{\infty}a^n\frac{d^nf(\vec{r}(t))}{dt^n}$ even means mathematically. Is there a closed-form expression for this term? If not, I'd be concerned about the numerical stability of calculating high-order numerical derivatives (as NDSolve does internally). $\endgroup$
    – Roman
    Aug 4, 2022 at 13:45

1 Answer 1


Not a complete solution; merely a simplification of the OP's right-hand side to get rid of high-order derivatives and replace them by a manageable formula.

Let's define a function $g(t)=f(\vec{r}(t))$ to simplify the notation a bit. We study the function $$ g_a(t) = \sum_{n=0}^{\infty}a^n\frac{d^ng(t)}{dt^n} = g(t) + \sum_{n=1}^{\infty}a^n\frac{d^ng(t)}{dt^n} $$ as used in the OP's differential equation, assuming $a\ge0$ in what follows. We note that $g_0(t)=\lim_{a\to0^+} g_a(t)=g(t)$.

Starting with the Fourier transform $$ g(t) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty}d\omega\, \hat{g}(\omega)e^{-i \omega t}\\ \hat{g}(\omega) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty}dt\, g(t)e^{i \omega t} $$ we compute $$ g_a(t) = \sum_{n=0}^{\infty}a^n\frac{d^n}{dt^n}\left( \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty}d\omega\, \hat{g}(\omega)e^{-i \omega t} \right)\\ = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty}d\omega\, \hat{g}(\omega) \left( \sum_{n=0}^{\infty}a^n\frac{d^n}{dt^n} e^{-i \omega t} \right) $$

Sum[a^n D[Exp[-I ω t], {t, n}], {n, 0, ∞}] // FullSimplify
(*    E^(-I ω t)/(1 + I ω a)    *)

$$ g_a(t) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty}d\omega\, \hat{g}(\omega) \left( \frac{e^{-i \omega t}}{1+i\omega a} \right)\\ = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty}d\omega\, \left( \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty}ds\, g(s)e^{i \omega s} \right) \left( \frac{e^{-i \omega t}}{1+i\omega a} \right)\\ = \int_{-\infty}^{\infty}ds\, g(s) \left( \frac{1}{2\pi} \int_{-\infty}^{\infty}d\omega\, \frac{e^{-i \omega (t-s)}}{1+i\omega a} \right) $$ With $\delta=t-s$ we can integrate the parenthesis using InverseFourierTransform:

1/Sqrt[2 π] InverseFourierTransform[1/(1 + I ω a), ω, δ, 
              Assumptions -> a > 0]
(*    E^(δ/a)*HeavisideTheta[-δ]/a    *)

and therefore $$ g_a(t) = \int_{-\infty}^{\infty}ds\, g(s) \left( \frac{e^{(t-s)/a} \theta(s-t)}{a} \right)\\ = \int_{t}^{\infty}ds\, \frac{g(s)e^{-(s-t)/a}}{a}\\ = \int_{0}^{\infty}d\tau\, \frac{g(t+\tau)e^{-\tau/a}}{a} $$ where $\tau=s-t$ was substituted. We see that the right-hand side of the OP's differential equation can be written as an exponentially decreasing average: the function $g_a(t)$ is an average of all $g(s)$ for $s\ge t$ with a weight that decreases exponentially as $s\to\infty$; the exponential decrease happens on a time-scale $a$. This exponential average is much easier to evaluate in practice than a sum over derivatives of very high order, especially when numerical data are involved (i.e. when $g(t)$ is known by numerical sampling).

The original differential equation is now $$ \frac{d^2}{dt^2} \vec{r}(t) = g_a(t) = \int_{0}^{\infty}d\tau\, \frac{e^{-\tau/a}}{a}f(\vec{r}(t+\tau)). $$ Solving this integro-differential equation is step 2 that would make this solution usable.

  • $\begingroup$ It is a very interesting approach that you are proposing and its very unique. I never though about this method that way. I will definately look at it and see how to adapt it to my problem. $\endgroup$
    – Aschoolar
    Aug 5, 2022 at 0:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.