13
$\begingroup$

I'd like to solve a differential equation $$ f^2(x) f''(x)=-x $$ where $f(x)$ is defined on $[0,1]$ and has a boundary condition $f(0)=f(1)=0$.

I somehow found out that the solution is fairly close to $f(x) = x^{1/3} \phi^{2/3}(\Phi^{-1}(1-x))$ where $\phi$ and $\Phi$ are pdf and cdf of a standard normal distribution, but it fails to solve the differential equation exactly.


Thank for all comments!

Based on the solution structure of Emden–Fowler Equation, I was able to identify the values of constants that satisfy the boundary conditions. The followings are the details:

Define \begin{equation} Z_R(\tau) \triangleq \sqrt{3} J_{1/3}(\tau) - Y_{1/3}(\tau) , \quad Z_L(\tau) \triangleq - \frac{2}{\pi} K_{1/3}(\tau) \end{equation} where $J, Y, K$ are Bessel functions. Further define \begin{equation} \bar{\tau} \triangleq \inf\{ \tau > 0; Z_R(\tau) = 0 \} \approx 2.3834 , \quad a \triangleq \frac{1}{ \bar{\tau}^{4/3} Z_R'(\bar{\tau})^2 } \approx 0.2910 , \quad b \triangleq a \left( \frac{9}{2} \right)^{1/3} \approx 0.1763. \end{equation} Then, the solution curve $\{ (x, f(x)) \}_{x \in [0,1]}$ is characterized by \begin{equation} \left\{ \left( x_R(\tau), y_R(\tau) \right) \right\}_{\tau \in [0, \bar{\tau}]} \bigcup \left\{ \left( x_L(\tau), y_L(\tau) \right) \right\}_{\tau \in [0, \infty]} \end{equation} where \begin{equation} x_R(\tau) \triangleq a \tau^{-2/3}\left[ \left( \tau Z_R'(\tau) + \frac{1}{3} Z_R(\tau) \right)^2 + \tau^2 Z_R(\tau)^2 \right] , \quad y_R(\tau) \triangleq b \tau^{2/3} Z_R(\tau)^2. \end{equation}

\begin{equation} x_L(\tau) \triangleq a \tau^{-2/3}\left[ \left( \tau Z_L'(\tau) + \frac{1}{3} Z_L(\tau) \right)^2 - \tau^2 Z_L(\tau)^2 \right] , \quad y_L(\tau) \triangleq b \tau^{2/3} Z_L(\tau)^2. \end{equation}


In addition to this analytic solution, I also obtained a numerical solution by repeatedly computing $$ f_{k+1}(x) \gets \left[ \left( f_k(x-2h) + f_k(x+2h) \right) + 4 \left(f_k(x-h)+f_k(x+h)\right) + \frac{8 x h^2}{f_k^2(x)} \right] \big/ 10 $$ on the grid $x \in \{2h,3h,\ldots,1-3h,1-2h\}$ for small $h$ with an initialization $f_0(x) \triangleq 0.5(1-(1-2x)^2)$.

The following figure shows these solutions:

enter image description here

$\endgroup$
  • $\begingroup$ $f(1)=0\ $ would imply $\ 0 = -1.\ $ $\endgroup$ – Wlod AA Mar 4 at 6:05
  • 2
    $\begingroup$ @WlodAA We will have $f''(1) = -\infty$. $\endgroup$ – Seungki Min Mar 4 at 6:22
  • $\begingroup$ I see. Thank you. $\endgroup$ – Wlod AA Mar 4 at 6:24
  • 6
    $\begingroup$ This is a case of the Emden-Fowler equation, see eqworld.ipmnet.ru/en/solutions/ode/ode0302.pdf, which arises in astrophysics and was subject of intensive study. E. Hille wrote that "it has fantastically complicated singularities near x=0". $\endgroup$ – Alexandre Eremenko Mar 4 at 11:51
24
$\begingroup$

Surprisingly, this case of the Emden-Fowler equation is explicitly solvable: see formula (2.3.27) in A. Polyanin and V. Zaitsev, Handbook of exact solutions of ordinary differential equations, Chapman & Hill, 2003.

I copy the formula, without verifying it. Let $$Z=C_1J_{1/3}(\tau)+C_2Y_{1/3}(\tau),$$ or $$Z=C_1I_{1/3}(\tau)+C_2K_{1/3}(\tau),$$ where $J,Y$ are Bessel and $I$, $K$ are modified Bessel functions. Then
$$x=a\tau^{-2/3}[(\tau Z^\prime+(1/3)Z)^2\pm\tau^2Z^2],\quad y=b\tau^{2/3}Z^2$$ satisfy $d^2y/dx^2=Axy^{-2}$ with $A=-(9/2)(b/a)^3.$

For the $+$ sign in $\pm$ take the first formula for $Z$, and for the $-$ the second one.

Remark. Emden-Fowler equation appears for the first time in the famous book by R. Emden, Gaskugeln (1907) and since then frequently arises in the study of stars and black holes.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ The first version with $Y$ seems to be correct, but I think it does not lead to a solution of the original problem. However, the version with $K$ instead of $Y$ does not seem to solve the differential equation (with either version of the $\pm$ sign). Perhaps there is some transcription error? $\endgroup$ – Neil Strickland Mar 4 at 13:46
  • 1
    $\begingroup$ Maple solves this using the Airy functions Ai and Bi, which are related to Bessel function of order $1/3, 2/3$. $\endgroup$ – Gerald Edgar Mar 4 at 16:45
  • 2
    $\begingroup$ @Neil Strickland: Yes, I made a misprint while copying from the book. Now corrected (second definition of Z has I instead of J). $\endgroup$ – Alexandre Eremenko Mar 4 at 18:39
8
$\begingroup$

I tried the following approach. Put $y=f(x)$ and $t=2x-1$ so the differential equation becomes $y^2\ddot{y}=-(t+1)/8$ with boundary conditions $y=0$ at $t=\pm 1$. We can then write $y=\sum_ia_it^i$. The differential equation gives a recurrence relation expressing all the coefficients $a_i$ in terms of $a_0$ and $a_1$. We can then truncate the power series to a given order $d$ and solve numerically for the boundary conditions. This seems to work in a well-behaved way, with good convergence at the endpoints and a result that is stable when we increase $d$. It looks like $a_0=0.450$ and $a_1=0.120$ to $3$ decimal places. Maple code is as follows:

with(plots):
Digits := 50:
d := 50:
y := add(a[i] * t^i,i=0..d):
sol0 := solve([coeffs(rem(expand(y^2 * diff(y,t,t) + (t+1)/8),t^(d-1),t),t)] 
              {seq(a[i],i=2..d)}):
y0 := expand(subs(sol0,y)):
sol1 := fsolve({subs(t= 1,y0),subs(t=-1,y0)},{a[0]=0.45,a[1]=0.1}):
aa[d] := subs(sol1,[a[0],a[1]]);
y1 := subs(sol1,y0);
y1x := subs(t = 2*x-1,y1):
Phi := unapply((1 + erf(x))/2,x):
phi := unapply(diff(Phi(x),x),x):
display(
 plot(y1x,x=0..1,colour=red),
 plot(x^(1/3) * phi(RootOf(1-x-Phi(_Z)))^(2/3),x=0..1,colour=blue)
);

This generates the following picture. The power series solution is in red and the function $x^{1/3}\phi(\Phi^{-1}(1-x))^{2/3}$ is in blue.

enter image description here

The coefficients $a_{2i}$ lie on a nice smooth curve, and the coefficients $a_{2i+1}$ lie on a similar curve shifted down slightly. Logs of the absolute values can be displayed as follows:

display(
  listplot([seq(log(-coeff(y1,t,2*i)),i=3..(d-1)/2)],style=point,colour=red),
  listplot([seq(log(-coeff(y1,t,2*i+1)),i=3..(d-1)/2)],style=point,colour=blue)
);

One could probably get further by finding an exact or approximate formula for these curves.

enter image description here

| cite | improve this answer | |
$\endgroup$

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.