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?

  • 2
    $\begingroup$ Why do you think that $F=G$? It seems that they can't be equal because $F(0)=1$ and $G(0)=0$. $\endgroup$ – Tomi27 Dec 15 '17 at 6:03
  • $\begingroup$ @Tomi27 Yes, you are right! $\endgroup$ – Eufisky Dec 15 '17 at 6:29
  • $\begingroup$ 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? $\endgroup$ – LutzL Dec 15 '17 at 9:52
  • $\begingroup$ 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$. $\endgroup$ – LutzL Dec 15 '17 at 10:20
  • $\begingroup$ @LutzL I'm trying to solve this equation, and it is also ok if you have numerical method. $\endgroup$ – 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


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.


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.