NDSolve issue with summation

Question

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?

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.