# NDSolve issue with summation

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?

• 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. Aug 4, 2022 at 13:31
• 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). Aug 4, 2022 at 13:45

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.

• 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. Aug 5, 2022 at 0:16