# Solution to differential equation $f^2(x) f''(x) = -x$ on [0,1]

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.

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 $$$$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)$$$$ where $$J, Y, K$$ are Bessel functions. Further define $$$$\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.$$$$ Then, the solution curve $$\{ (x, f(x)) \}_{x \in [0,1]}$$ is characterized by $$$$\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]}$$$$ where $$$$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.$$$$

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

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:

• $f(1)=0\$ would imply $\ 0 = -1.\$ – Wlod AA Mar 4 at 6:05
• @WlodAA We will have $f''(1) = -\infty$. – Seungki Min Mar 4 at 6:22
• I see. Thank you. – Wlod AA Mar 4 at 6:24
• 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". – Alexandre Eremenko Mar 4 at 11:51

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.

• 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? – Neil Strickland Mar 4 at 13:46
• Maple solves this using the Airy functions Ai and Bi, which are related to Bessel function of order $1/3, 2/3$. – Gerald Edgar Mar 4 at 16:45
• @Neil Strickland: Yes, I made a misprint while copying from the book. Now corrected (second definition of Z has I instead of J). – Alexandre Eremenko Mar 4 at 18:39

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

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.