1
$\begingroup$

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?

$\endgroup$
2
  • 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 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 at 13:45

1 Answer 1

1
$\begingroup$

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

$\endgroup$
1
  • $\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 at 0:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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