# Second-Order Ordinary Differential Equation System

If we have second-order ordinary differential equation system $$\begin{cases} F''(x)=(F(x))^3+F(x)(G(x))^2\\ G''(x)=2G(x)(F(x))^2 \end{cases}$$ and it satisfies boundary conditions $$F(0)=G'(0)=1,\qquad F'(0)=G(0)=0.$$ Find $F\left(\frac{\pi}{3}\right)$.

It seems that $F=G$, we have $$F''=2F^3\Rightarrow F'F''=2F'F^3\Rightarrow \frac{1}{2}\left( F' \right) ^2=\frac{1}{2}F^4+\frac{1}{2}C_1,$$ that is $F'=\pm \sqrt{F^4+C_1}$, So $$\int_0^x{\frac{1}{\sqrt{F^4+C_1}}dt}=\pm x+C_2.$$ But I think it may have some problems with the boundary conditions.

Edit: In fact, $F(x)=1/\cos x,G(x)=\tan x$, but how can we prove it?

• Why do you think that $F=G$? It seems that they can't be equal because $F(0)=1$ and $G(0)=0$. – Tomi27 Dec 15 '17 at 6:03
• @Tomi27 Yes, you are right! – Eufisky Dec 15 '17 at 6:29
• What exactly is the context of this question? Is it to try out some exotic symbolic solution method or is it to explore the possibilities and limits of numerical integrators? – LutzL Dec 15 '17 at 9:52
• You can get $$F'(x)^2+\frac12G'(x)^2-\frac12F(x)^4-F(x)^2G(x)^2$$ as a constant of the dynamics. Numerical integration shows that the result is more or less exactly $2$. – LutzL Dec 15 '17 at 10:20
• @LutzL I'm trying to solve this equation, and it is also ok if you have numerical method. – Eufisky Dec 15 '17 at 11:42

Using the python scipy integrators with the script

from scipy.integrate import odeint
def odefunc(y,t):
F,dF,G,dG = y;
return [ dF, F**3+F*G**2, dG, 2*G*F**2 ]
y0 = [1., 0., 0., 1.]
sol = odeint(odefunc, y0, [0,np.pi/3],atol=1e-14, rtol=1e-13,mxstep=10000)
print "%.14f"%sol[-1,0]


gives the value $F(\pi/3)$ as

2.00000000000097


Other software with numerical integrators should work in a similar way.

Or modify to

x = np.linspace(0,1,11)*np.pi/3
sol = odeint(ode, y0, x,atol=1e-14, rtol=1e-13,mxstep=10000)
def E(y):
F,dF,G,dG = y;
return dF**2 + 0.5*dG**2 - 0.5*F**4 - (F*G)**2
for xx,y in zip(x,sol):
print "x=%.2f*pi/3, F(x)=%.14f, E(x)=%.14f"%(xx*3/np.pi, y[0], E(y))


to get the result

x=0.00*pi/3, F(x)=1.00000000000000, E(x)=0.00000000000000
x=0.10*pi/3, F(x)=1.00550827956352, E(x)=0.00000000000001
x=0.20*pi/3, F(x)=1.02234059486504, E(x)=0.00000000000002
x=0.30*pi/3, F(x)=1.05146222423829, E(x)=0.00000000000003
x=0.40*pi/3, F(x)=1.09463627850608, E(x)=0.00000000000005
x=0.50*pi/3, F(x)=1.15470053837930, E(x)=0.00000000000013
x=0.60*pi/3, F(x)=1.23606797749988, E(x)=0.00000000000047
x=0.70*pi/3, F(x)=1.34563272960657, E(x)=0.00000000000225
x=0.80*pi/3, F(x)=1.49447654986494, E(x)=0.00000000000307
x=0.90*pi/3, F(x)=1.70130161670461, E(x)=0.00000000000392
x=1.00*pi/3, F(x)=2.00000000000092, E(x)=0.00000000001329


where the energy function as theoretically conserved quantity indicates the error level of the numerical result.