Solve singular fourth order ODE

Question

I am interested in solving the following differential equation : $$ \left(\frac{300 i}{r^{5}}-\frac{3}{r^{4}}\right) \psi(r)+\left(-\frac{45}{r^{6}}+\frac{3}{r^{3}}\right) \psi^{\prime}(r)+\left(\frac{21}{r^{5}}-\frac{3}{r^{2}}\right) \psi^{\prime \prime}(r)+\left(-\frac{3}{r^{4}}+\frac{2}{r}\right) \psi^{(3)}(r)+\psi^{(4)}(r)=0 $$ with the boundary conditions $\psi\to r$ as $r\to \infty$ and $\psi=0$ at $r=0$.

DSolve cannot solve this, so I must resort to NDSolve. I am not aware of a consistent way to regularize the ODE at the origin. I was thinking that near the origin, I can neglect the $\mathcal{O}(1/r^7)$ terms and exactly solve the ODE $$\left(\frac{300 i}{r^{5}}-\frac{3}{r^{4}}\right) \psi(r)+\left(\frac{3}{r^{3}}\right) \psi^{\prime}(r)+\left(-\frac{3}{r^{2}}\right) \psi^{\prime \prime}(r)+\left(\frac{2}{r}\right) \psi^{(3)}(r)+\psi^{(4)}(r)=0$$ and use the derivatives of this solution near $r\sim\epsilon$ where $\epsilon$ is some small number.

Experts, is there a consistent way to deal with singular ODEs? Any help is appreciated

Edit:

sol = NDSolve[{r^5 (((300 I)/r^5 - 3/r^4) \[Psi][
         r] + (-(45/r^6) + 3/r^3) Derivative[1][\[Psi]][
         r] + (21/r^5 - 3/r^2)  Derivative[2][\[Psi]][
         r] + (-(3/r^4) + 2/r) Derivative[3][\[Psi]][r] + 
       Derivative[4][\[Psi]][r]) == 0,  \[Psi][1] == 1,  
   Derivative[1][\[Psi]][$MinMachineNumber] == 0,  
   Derivative[2][\[Psi]][$MinMachineNumber] == 
    0, \[Psi][$MinMachineNumber] == 
    0 }, \[Psi], {r, $MinMachineNumber, 1}]

I added some code that does not work because the ODE at the origin is highly singular

Answer 1

eqn = r^5 (((300 I)/r^5 - 3/r^4) \[Psi][r] + (-(45/r^6) + 3/r^3) Derivative[1][\[Psi]][r] 
            + (21/r^5 - 3/r^2) Derivative[2][\[Psi]][r] + (-(3/r^4) + 2/r) Derivative[3][\[Psi]][r] + Derivative[4][\[Psi]][r]) == 0

Let's see if Mathematica can handle the asymptotic expansions as $r \to 0$ and $r \to \infty$:

smallsol[r_] = AsymptoticDSolveValue[eqn, \[Psi][r], r -> 0]
largesol[r_] = AsymptoticDSolveValue[eqn, \[Psi][r], r -> \[Infinity], GeneratedParameters -> B]

(* C[1] + E^(-(1/r^3)) r^3 C[2] + r^4 C[3] + r^6 C[4] *)
(* 6 r^3 B[1] + r B[2] + B[4]/r - 2 r B[3] Log[r] *)

The requirement that $\psi \to 0$ as $r \to 0$ requires that $C_1 = 0$; the requirement that $\psi \to r$ as $r \to \infty$ requires that $B_1 = B_3 = 0$ and $B_2 = 1$. This means that out of four degrees of freedom of the solution, three are determined by these boundary conditions. The remaining degree of freedom corresponds (I think) to the fact that this ODE is linear.

To parlay these asymptotic solutions into a set of boundary conditions, we can do the following:

Simplify[Table[(D[smallsol[r] - \[Psi][r], {r, n}] == 0) /. r -> eps, {n, 0, 3}] /. C[1] -> 0]
smallrBC = (List @@ Eliminate[%, {C[2], C[3], C[4]}])[[1]]

This creates a list of equations relating the constants $C_i$ to the derivatives of $\psi$ at $r = \epsilon$, sets $C_1 = 0$, and then finds the relations that the derivatives of $\psi$ must satisfy with $C_1$ constrained to be zero. Similarly, we have at "infinity"

Simplify[Table[(D[largesol[r] - \[Psi][r], {r, n}] == 0) /. r -> inf, {n, 0, 3}] /. {B[1] -> 0, B[3] -> 0, B[2] -> 1}]
largerBC = (List @@ Eliminate[%, {B[4]}])[[1 ;; 3]]

With all that, will Mathematica solve it for us? The answer is "sort of":

bdvals = {eps -> 10^-2, inf -> 10^2}
sol = NDSolve[
        Join[{eqn}, smallrBC, largerBC] /. bdvals, \[Psi][r], {r, inf, eps} /. bdvals, 
        Method -> {"Shooting", 
                   "StartingInitialConditions" -> {\[Psi][inf] == inf} /. bdvals}]
Plot[Evaluate[ReIm[\[Psi][r] /. sol]], {r, eps, inf} /. bdvals]

This code does return an error, though: "The equations derived from the boundary conditions are numerically ill-conditioned. The boundary conditions may not be sufficient to uniquely define a solution." What this indicates is that Mathematica was not able to reliably find a solution that satisfies the boundary conditions at "zero" and "infinity" due to the fact that both $r = 0$ and $r = \infty$ are singular, and so the coefficients in our boundary conditions (particularly the one at $r = \epsilon$) are difficult to solve numerically. (The problem is even worse if we let Mathematica "shoot from zero" rather than "shoot from infinity"; in that case, it doesn't return a solution at all.)

Answer 2

This is an incomplete answer, but perhaps you can combine it with the answer of @MichaelSeifert. It is more manual in a way, but perhaps allows you to understand where the numerical issues come from. So this being a linear homogeneous equation, for all $s,t>0$ there is a $4 \times 4$ matrix $P(s,t)$ such that for all solutions of the ODE one has $$ \begin{pmatrix} \psi(s) \\ \psi'(s) \\ \psi''(s) \\ \psi'''(s) \end{pmatrix} \;=\; P(s,t)\begin{pmatrix} \psi(t) \\ \psi'(t) \\ \psi''(t) \\ \psi'''(t) \end{pmatrix} $$ and I am calling this $P$ for propagator. One can construct this matrix numerically in Mathematica as follows (not sure if there is a more direct command, I simply call NDSolve a few times):

eqn=r^5 (((300 I)/r^5-3/r^4) \[Psi][r]+(-(45/r^6)+3/r^3) Derivative[1][\[Psi]][r]+(21/r^5-3/r^2) Derivative[2][\[Psi]][r]+(-(3/r^4)+2/r) Derivative[3][\[Psi]][r]+Derivative[4][\[Psi]][r])==0;
sol[s_,t_][vec0_]:={\[Psi][s],\[Psi]'[s],\[Psi]''[s],\[Psi]'''[s]}/.First[NDSolve[Join[{eqn},Thread[{\[Psi][t],\[Psi]'[t],\[Psi]''[t],\[Psi]'''[t]}==vec0]],\[Psi],{r,t,s},AccuracyGoal->12,PrecisionGoal->12]];
P[s_,t_]:=Transpose[Map[sol[s,t],IdentityMatrix[4]]];

For example, P[1/2,2] evaluates to:

{{-0.500673-4.04689 I,-0.600881+2.27858 I,0.550522 -0.884086 I,-0.370662+0.327502 I},
 {8.81853 +5.23293 I,-3.71234-5.23891 I,1.09816 +2.77047 I,-0.221186-1.33511 I},
 {-28.4567+24.5997 I,18.2139 -8.63023 I,-7.76148+0.821715 I,3.56598 +1.28922 I},
 {-89.7965-85.4907 I,2.74802 +85.5374 I,8.55136 -49.2012 I,-3.68842+24.3657 I}}

I have not checked this carefully. Note that the exact propagator satisfies $P(s,t)P(t,u) = P(s,u)$ for all $s,t,u>0$ and one can check for various values of $s,t,u>0$ that the numerically constructed $P(s,t)$ satisfies this more or less, try P[1/2,2].P[2,1/2] for example. If you need more precision, you can adjust the settings, but there are of course limits with 64 bit numbers at least.

Now With[{eps=10^(-2),inf=10^2}, P[eps,inf]] evaluates to:

{{144.35 +656.623 I,4760.13 -49174.1 I,-854332.+2.53709*10^6 I,5.3008*10^7-1.07601*10^8 I},
 {-82.0461+17.4275 I,6130.6 +637.386 I,-315763.-108830. I,1.3377*10^7+6.70836*10^6 I},
 {-8195.79+1631.87 I,609936. +71291.1 I,-3.13188*10^7-1.12433*10^7 I,1.32413*10^9+6.85308*10^8 I},
 {3484.19 -28965.3 I,-898781.+1.94636*10^6 I,7.13089*10^7-9.17075*10^7 I,-3.71148*10^9+3.64932*10^9 I}}

This will certainly contain significant numerical errors. This matrix is very ill-conditioned, its singular values are {5.41708*10^9,150.532,0.0015424,2.38302*10^-8}. It is of course a mathematical fact that the true propagator $P(s,t)$ exists and is invertible for all $s,t>0$.

Anyhow, in principle, this matrix gives you 4 linear homogeneous equations relating the solution and its first three derivatives at eps and at inf. You could solve this together with the (approximate) equations smallrBC and largerBC in the answer of @MichaelSeifert, that is just linear algebra and can easily be analyzed. The question is whether you find any choices of eps and inf and of the solver settings where all errors are under control: the numerical errors in the propagator, and the approximations involved in smallrBC and largerBC.

Two possible improvements:

• You could improve the approximations on the $r$-intervals $(0,\text{eps}]$ and $[\text{inf},\infty)$ say write down more terms of a converging expansion of the solution and so on, perhaps AsymptoticDSolveValue can help you here (I did not know about this command!). If you have good expansions, then you do not have to make eps quite so small and you do not have to make inf quite so big.
• You could improve the conditioning of the differential equations itself. The known leading asymptotic behaviour (on either side) can be used to change to other variables where the conditioning is good.

Hope this helps.