# 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 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 at 13:45

Not an answer; too long for a comment; here's some preliminary mathematical work.

Let's take a function $$f(x)$$ and compute $$F_a(x)=\sum_{n=0}^{\infty} a^n \frac{d^n f(x)}{dx^n}$$. Starting with the Fourier transform

$$f(x) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} \hat{f}(k)e^{i k x}dk\\ \hat{f}(k) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} f(x)e^{-i k x}dx$$ Now we can apply the derivatives-transformation to the term $$e^{i k x}$$:

$$\sum_{n=0}^{\infty} a^n \frac{d^n e^{i k x}}{dx^n} = \frac{e^{i k x}}{1-i a k}$$

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


which makes the function $$F_a(x)$$ $$F_a(x) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} \hat{f}(k)\frac{e^{i k x}}{1-i a k}dk = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} \hat{f}(k)\hat{g}(k)e^{i k x}dk$$ in terms of a helper function $$\hat{g}_a(k)=\frac{1}{1-i a k}$$ that is the Fourier transform of $$g_a(x) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} \hat{g}_a(k)e^{i k x}dk = \frac{\sqrt{2\pi}}{a} e^{x/a} \theta(-x)$$ in terms of the Heaviside function $$\theta$$.

Assuming[a > 0, FourierTransform[1/(1 - I a k), k, x]]
(*    E^(x/a) Sqrt[2π] HeavisideTheta[-x]/a    *)


The function $$F_a(x)$$ is therefore the convolution of $$f(x)$$ with $$g_a(x)$$:

$$F_a(x) = (f\star g_a)(x) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} f(y)g_a(-x-y)dy\\ = \frac{1}{a}\int_{-\infty}^{\infty} f(y)e^{-(x+y)/a} \theta(x+y)dy = \frac{e^{-x/a}}{a}\int_{-x}^{\infty} f(y)e^{-y/a}dy$$ Maybe this derivation can be a starting point for a derivation of your infinite sum, to provide a way of evaluating it without recourse to derivatives of infinite order.

(PS: it's very likely that I got some signs wrong somewhere in the derivation.)

• 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 at 0:16