Duffing.eq.Bogolyubov_method

173 days ago by blinkenlichten

var('x t psi a w xd epsilon') 
        
\newcommand{\Bold}[1]{\mathbf{#1}}\left(x, t, \psi, a, w, \mbox{xd}, \epsilon\right)
\newcommand{\Bold}[1]{\mathbf{#1}}\left(x, t, \psi, a, w, \mbox{xd}, \epsilon\right)
def duffing_first_asymp(func): """ a function for evaluating first asymptotic approximation of solution for Duffing's equation via use of Krilov-bogolyubov method. Accepts 1 arguments: \param1 a function, which is dependent on 2 variables, for example f(u,ud)=u^2-ud^5/4, where u -- is a time-depenent function u(t) in Duffing's equation and ud -- it's derivative: ud(t)=du/dt """ a = var('a') epsilon = var('epsilon') psi=var('psi') t=var('t') u=var('u') ud=var('ud') def f1(a): return 1/pi*integrate(sin(psi)*func(a*cos(psi),-a*w*sin(psi)),psi,0,2*pi) def g1(a): return 1/pi*integrate(cos(psi)*func(a*cos(psi),-a*w*sin(psi)),psi,0,2*pi) def dadt(t): return -epsilon/(2*w)*f1(a) def dpsidt(t): return -epsilon/(2*w)*g1(a) print "Krilov-Bogolyubov first asymp. approximation for equation:\n u''+w^2*u=epsilon*f(u,u');\n f(u,ud)=",f(u,ud) print '\n results:' print 'da/dt(t)=' show(dadt(t)) print 'dpsi/dt(t)=' show(dpsidt(t)) print 'u(t)=a(t)*cos(psi(t))' 
       
print "problem 1: x''+3x+1/2x^3+1/2x^5=0" print "epsilon=1/2" def f(x,xd): return -(x^3+2/3*x^5) print "f(x,xd)=",f(x,xd) duffing_first_asymp(f) 
        
problem 1:  x''+3x+1/2x^3+1/2x^5=0
epsilon=1/2
f(x,xd)= -2/3*x^5 - x^3
Krilov-Bogolyubov first asymp. approximation for equation:
 u''+w^2*u=epsilon*f(u,u');
 f(u,ud)= -2/3*u^5 - u^3

 results:
da/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}0

dpsi/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(5 \, \pi a^{5} + 9 \, \pi a^{3}\right)} \epsilon}{24 \, \pi w}

u(t)=a(t)*cos(psi(t))
problem 1:  x''+3x+1/2x^3+1/2x^5=0
epsilon=1/2
f(x,xd)= -2/3*x^5 - x^3
Krilov-Bogolyubov first asymp. approximation for equation:
 u''+w^2*u=epsilon*f(u,u');
 f(u,ud)= -2/3*u^5 - u^3

 results:
da/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}0

dpsi/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(5 \, \pi a^{5} + 9 \, \pi a^{3}\right)} \epsilon}{24 \, \pi w}

u(t)=a(t)*cos(psi(t))
def krilov_bogolyubov_second_noxd(func,NMAX): """ A function for evaluating second asymptotic approximation of solution for Duffing equation u''+w^2*u=epsilon*f(u). where arbitrary function F(x,x') depends only on x: x->F(x) . by Krilov-bogolyubov method. Accepts 2 arguments: \param1 a function and its derivative, function is dependent on 1 variable only, for example f(u)=u^2, where u -- is a time-depenent function u(t) in Duffing's equation \param2 Number of maximal possible fourier series approximation. What is the output? It gives you these expressions da(t)/dt = epsilon*A1+epsilon^2*A2 ----> whereas a(t)=solution of this diff. eq.(you have to solve) dpsi(t)/dt = w+epsilon*B1+epsilon^2*B2 ----> whereas psi(t)=solution of this diff. eq.(you have to solve) u1(psi(t)) is calculated Finally aysmptotic aproximation of solution is given by formula: x(t)=a(t)*cos(psi(t))+u1(psi(t)) """ a = var('a') w = var('w') epsilon = var('epsilon') psi=var('psi') t=var('t') u=var('u') ud=var('ud') print "Krilov-Bogolyubov secondt approximation for equation:\n d^2 u/dt+w^2*u=epsilon*f(u);\n f(u)=",func(u) CN=[ ] for n in range(NMAX): CN.append(1/pi*integrate(func(a*cos(psi))*cos(n*psi),psi,0,2*pi)) if len(CN)<1: CN.append(0) CN[0]=1/(2*pi)*integrate(func(a*cos(psi)),psi,0,2*pi) print 'Cn[i]=' show(CN) GN = [ ] for n in range(NMAX): GN.append(-1*CN[n]) print 'gn=' show(GN) HN = [] for n in range(NMAX): HN.append(0) A1=0 B1=CN[1]/(2*w*a) print 'A1=' show(A1) print 'B1=' show(B1) def u1(a,psi): res=0 for n in range(NMAX): if n!=1: res+=CN[n]*cos(n*psi)/(n^2-1) return res/w^2 print 'u1(a,psi)=' show(u1(a,psi)) func_dx = diff(func) A2=0 B2=-1/(2*w)*(CN[1]/(2*w*a))^2 for n in range(NMAX): if n!=1 : B2+= 1/(2*w*pi*a)*CN[n]/(w^2*(n^2-1))*integrate(func_dx(a*cos(psi))*cos(psi)*cos(n*psi),psi,0,2*pi) print 'A2=' show (A2) print 'B2=' show (B2) def dadt(t): return 0 def dpsidt(t): return w+epsilon*B1+epsilon^2*B2 def psiint(t): return dpsidt(t).integrate(t) print '\nRESULTS:' print 'da/dt(t)=' show(dadt(t)) print 'dpsi/dt(t)=' show(dpsidt(t)) """resulting function x(t)""" def resu(a,psi): sum=0 for n in range(2,NMAX): sum+=epsilon/w^2*CN[n]*cos(n*psi)/(n^2-1) return a*cos(psi)-epsilon*CN[0]/w^2+sum print 'x(t)=' show(resu(a,psi)) print 'where a(t) and psi(t) are obtained from equations above' print 'psi(t)=' show(psiint(t)) return [A1,B1,A2,B2] 
       
print "========================================" print "|| equation1: x''+3x+1/2x^3+1/2x^5=0 " print "========================================" print "epsilon=1/2" f(x)=(x^3+2/3*x^5) print "f(x,xd)=",f(x) show(f(x)) krilov_bogolyubov_second_noxd(f,10) 
        
========================================
|| equation1: x''+3x+1/2x^3+1/2x^5=0    
========================================
epsilon=1/2
f(x,xd)= 2/3*x^5 + x^3

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{2}{3} \, x^{5} + x^{3}

Krilov-Bogolyubov secondt approximation for equation:
 d^2 u/dt+w^2*u=epsilon*f(u);
 f(u)= 2/3*u^5 + u^3
Cn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, \frac{5 \, \pi a^{5} + 9 \, \pi a^{3}}{12 \, \pi}, 0, \frac{5 \, \pi a^{5} + 6 \, \pi a^{3}}{24 \, \pi}, 0, \frac{1}{24} \, a^{5}, 0, 0, 0, 0\right]

gn=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, -\frac{5 \, \pi a^{5} + 9 \, \pi a^{3}}{12 \, \pi}, 0, -\frac{5 \, \pi a^{5} + 6 \, \pi a^{3}}{24 \, \pi}, 0, -\frac{1}{24} \, a^{5}, 0, 0, 0, 0\right]

A1=

\newcommand{\Bold}[1]{\mathbf{#1}}0

B1=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{5 \, \pi a^{5} + 9 \, \pi a^{3}}{24 \, \pi a w}

u1(a,psi)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{a^{5} \cos\left(5 \, \psi\right) + \frac{3 \, {\left(5 \, \pi a^{5} + 6 \, \pi a^{3}\right)} \cos\left(3 \, \psi\right)}{\pi}}{576 \, w^{2}}

A2=

\newcommand{\Bold}[1]{\mathbf{#1}}0

B2=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{5 \, a^{8}}{27648 \, w^{3}} + \frac{{\left(25 \, \pi a^{4} + 18 \, \pi a^{2}\right)} {\left(5 \, \pi a^{5} + 6 \, \pi a^{3}\right)}}{9216 \, \pi^{2} a w^{3}} - \frac{{\left(5 \, \pi a^{5} + 9 \, \pi a^{3}\right)}^{2}}{1152 \, \pi^{2} a^{2} w^{3}}


RESULTS:
da/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}0

dpsi/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{27648} \, {\left(\frac{5 \, a^{8}}{w^{3}} + \frac{3 \, {\left(25 \, \pi a^{4} + 18 \, \pi a^{2}\right)} {\left(5 \, \pi a^{5} + 6 \, \pi a^{3}\right)}}{\pi^{2} a w^{3}} - \frac{24 \, {\left(5 \, \pi a^{5} + 9 \, \pi a^{3}\right)}^{2}}{\pi^{2} a^{2} w^{3}}\right)} \epsilon^{2} + w + \frac{{\left(5 \, \pi a^{5} + 9 \, \pi a^{3}\right)} \epsilon}{24 \, \pi a w}

x(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{a^{5} \epsilon \cos\left(5 \, \psi\right)}{576 \, w^{2}} + a \cos\left(\psi\right) + \frac{{\left(5 \, \pi a^{5} + 6 \, \pi a^{3}\right)} \epsilon \cos\left(3 \, \psi\right)}{192 \, \pi w^{2}}

where a(t) and psi(t) are obtained from equations above
psi(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{27648} \, {\left({\left(\frac{5 \, a^{8}}{w^{3}} + \frac{3 \, {\left(25 \, \pi a^{4} + 18 \, \pi a^{2}\right)} {\left(5 \, \pi a^{5} + 6 \, \pi a^{3}\right)}}{\pi^{2} a w^{3}} - \frac{24 \, {\left(5 \, \pi a^{5} + 9 \, \pi a^{3}\right)}^{2}}{\pi^{2} a^{2} w^{3}}\right)} \epsilon^{2} + 27648 \, w + \frac{1152 \, {\left(5 \, \pi a^{5} + 9 \, \pi a^{3}\right)} \epsilon}{\pi a w}\right)} t

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, \frac{5 \, \pi a^{5} + 9 \, \pi a^{3}}{24 \, \pi a w}, 0, \frac{5 \, a^{8}}{27648 \, w^{3}} + \frac{{\left(25 \, \pi a^{4} + 18 \, \pi a^{2}\right)} {\left(5 \, \pi a^{5} + 6 \, \pi a^{3}\right)}}{9216 \, \pi^{2} a w^{3}} - \frac{{\left(5 \, \pi a^{5} + 9 \, \pi a^{3}\right)}^{2}}{1152 \, \pi^{2} a^{2} w^{3}}\right]
========================================
|| equation1: x''+3x+1/2x^3+1/2x^5=0    
========================================
epsilon=1/2
f(x,xd)= 2/3*x^5 + x^3

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{2}{3} \, x^{5} + x^{3}

Krilov-Bogolyubov secondt approximation for equation:
 d^2 u/dt+w^2*u=epsilon*f(u);
 f(u)= 2/3*u^5 + u^3
Cn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, \frac{5 \, \pi a^{5} + 9 \, \pi a^{3}}{12 \, \pi}, 0, \frac{5 \, \pi a^{5} + 6 \, \pi a^{3}}{24 \, \pi}, 0, \frac{1}{24} \, a^{5}, 0, 0, 0, 0\right]

gn=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, -\frac{5 \, \pi a^{5} + 9 \, \pi a^{3}}{12 \, \pi}, 0, -\frac{5 \, \pi a^{5} + 6 \, \pi a^{3}}{24 \, \pi}, 0, -\frac{1}{24} \, a^{5}, 0, 0, 0, 0\right]

A1=

\newcommand{\Bold}[1]{\mathbf{#1}}0

B1=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{5 \, \pi a^{5} + 9 \, \pi a^{3}}{24 \, \pi a w}

u1(a,psi)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{a^{5} \cos\left(5 \, \psi\right) + \frac{3 \, {\left(5 \, \pi a^{5} + 6 \, \pi a^{3}\right)} \cos\left(3 \, \psi\right)}{\pi}}{576 \, w^{2}}

A2=

\newcommand{\Bold}[1]{\mathbf{#1}}0

B2=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{5 \, a^{8}}{27648 \, w^{3}} + \frac{{\left(25 \, \pi a^{4} + 18 \, \pi a^{2}\right)} {\left(5 \, \pi a^{5} + 6 \, \pi a^{3}\right)}}{9216 \, \pi^{2} a w^{3}} - \frac{{\left(5 \, \pi a^{5} + 9 \, \pi a^{3}\right)}^{2}}{1152 \, \pi^{2} a^{2} w^{3}}


RESULTS:
da/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}0

dpsi/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{27648} \, {\left(\frac{5 \, a^{8}}{w^{3}} + \frac{3 \, {\left(25 \, \pi a^{4} + 18 \, \pi a^{2}\right)} {\left(5 \, \pi a^{5} + 6 \, \pi a^{3}\right)}}{\pi^{2} a w^{3}} - \frac{24 \, {\left(5 \, \pi a^{5} + 9 \, \pi a^{3}\right)}^{2}}{\pi^{2} a^{2} w^{3}}\right)} \epsilon^{2} + w + \frac{{\left(5 \, \pi a^{5} + 9 \, \pi a^{3}\right)} \epsilon}{24 \, \pi a w}

x(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{a^{5} \epsilon \cos\left(5 \, \psi\right)}{576 \, w^{2}} + a \cos\left(\psi\right) + \frac{{\left(5 \, \pi a^{5} + 6 \, \pi a^{3}\right)} \epsilon \cos\left(3 \, \psi\right)}{192 \, \pi w^{2}}

where a(t) and psi(t) are obtained from equations above
psi(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{27648} \, {\left({\left(\frac{5 \, a^{8}}{w^{3}} + \frac{3 \, {\left(25 \, \pi a^{4} + 18 \, \pi a^{2}\right)} {\left(5 \, \pi a^{5} + 6 \, \pi a^{3}\right)}}{\pi^{2} a w^{3}} - \frac{24 \, {\left(5 \, \pi a^{5} + 9 \, \pi a^{3}\right)}^{2}}{\pi^{2} a^{2} w^{3}}\right)} \epsilon^{2} + 27648 \, w + \frac{1152 \, {\left(5 \, \pi a^{5} + 9 \, \pi a^{3}\right)} \epsilon}{\pi a w}\right)} t

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, \frac{5 \, \pi a^{5} + 9 \, \pi a^{3}}{24 \, \pi a w}, 0, \frac{5 \, a^{8}}{27648 \, w^{3}} + \frac{{\left(25 \, \pi a^{4} + 18 \, \pi a^{2}\right)} {\left(5 \, \pi a^{5} + 6 \, \pi a^{3}\right)}}{9216 \, \pi^{2} a w^{3}} - \frac{{\left(5 \, \pi a^{5} + 9 \, \pi a^{3}\right)}^{2}}{1152 \, \pi^{2} a^{2} w^{3}}\right]
print "=================================" print "|| equation2: x''+3x+1/2x^2=0 " print "=================================" print "epsilon=1/2" f(x)=x^2 print "f(x,xd)=",f(x) show(f(x)) krilov_bogolyubov_second_noxd(f,10) 
        
=================================
|| equation2: x''+3x+1/2x^2=0    
=================================
epsilon=1/2
f(x,xd)= x^2

\newcommand{\Bold}[1]{\mathbf{#1}}x^{2}

Krilov-Bogolyubov secondt approximation for equation:
 d^2 u/dt+w^2*u=epsilon*f(u);
 f(u)= u^2
Cn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[\frac{1}{2} \, a^{2}, 0, \frac{1}{2} \, a^{2}, 0, 0, 0, 0, 0, 0, 0\right]

gn=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{1}{2} \, a^{2}, 0, -\frac{1}{2} \, a^{2}, 0, 0, 0, 0, 0, 0, 0\right]

A1=

\newcommand{\Bold}[1]{\mathbf{#1}}0

B1=

\newcommand{\Bold}[1]{\mathbf{#1}}0

u1(a,psi)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{a^{2} \cos\left(2 \, \psi\right) - 3 \, a^{2}}{6 \, w^{2}}

A2=

\newcommand{\Bold}[1]{\mathbf{#1}}0

B2=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{5 \, a^{2}}{12 \, w^{3}}


RESULTS:
da/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}0

dpsi/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{5 \, a^{2} \epsilon^{2}}{12 \, w^{3}} + w

x(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{a^{2} \epsilon \cos\left(2 \, \psi\right)}{6 \, w^{2}} + a \cos\left(\psi\right) - \frac{a^{2} \epsilon}{2 \, w^{2}}

where a(t) and psi(t) are obtained from equations above
psi(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{12} \, {\left(\frac{5 \, a^{2} \epsilon^{2}}{w^{3}} - 12 \, w\right)} t

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, 0, 0, -\frac{5 \, a^{2}}{12 \, w^{3}}\right]
=================================
|| equation2: x''+3x+1/2x^2=0    
=================================
epsilon=1/2
f(x,xd)= x^2

\newcommand{\Bold}[1]{\mathbf{#1}}x^{2}

Krilov-Bogolyubov secondt approximation for equation:
 d^2 u/dt+w^2*u=epsilon*f(u);
 f(u)= u^2
Cn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[\frac{1}{2} \, a^{2}, 0, \frac{1}{2} \, a^{2}, 0, 0, 0, 0, 0, 0, 0\right]

gn=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{1}{2} \, a^{2}, 0, -\frac{1}{2} \, a^{2}, 0, 0, 0, 0, 0, 0, 0\right]

A1=

\newcommand{\Bold}[1]{\mathbf{#1}}0

B1=

\newcommand{\Bold}[1]{\mathbf{#1}}0

u1(a,psi)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{a^{2} \cos\left(2 \, \psi\right) - 3 \, a^{2}}{6 \, w^{2}}

A2=

\newcommand{\Bold}[1]{\mathbf{#1}}0

B2=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{5 \, a^{2}}{12 \, w^{3}}


RESULTS:
da/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}0

dpsi/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{5 \, a^{2} \epsilon^{2}}{12 \, w^{3}} + w

x(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{a^{2} \epsilon \cos\left(2 \, \psi\right)}{6 \, w^{2}} + a \cos\left(\psi\right) - \frac{a^{2} \epsilon}{2 \, w^{2}}

where a(t) and psi(t) are obtained from equations above
psi(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{12} \, {\left(\frac{5 \, a^{2} \epsilon^{2}}{w^{3}} - 12 \, w\right)} t

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, 0, 0, -\frac{5 \, a^{2}}{12 \, w^{3}}\right]
def krilov_bogolyubov_second_nox(func,NMAX): """ A function for evaluating second asymptotic approximation of solution for Duffing equation u''+w^2*u=epsilon*f(u'). where arbitrary function F(u,u') depends only on u': u'->F(u') . by Krilov-bogolyubov method. Accepts 2 arguments: \param1 : a function, function is dependent on 1 variable only, for example f(ud)=ud^2, where ud -- is a derivative of time-depenent function u(t) in Duffing's equation \param2 Number of maximal possible fourier series approximation What is the output? It gives you these expressions da(t)/dt = epsilon*A1+epsilon^2*A2 ----> whereas a(t)=solution of this diff. eq.(you have to solve) dpsi(t)/dt = w+epsilon*B1+epsilon^2*B2 ----> whereas psi(t)=solution of this diff. eq.(you have to solve) u1(psi(t)) is calculated Finally aysmptotic aproximation of solution is given by formula: x(t)=a(t)*cos(psi(t))+u1(psi(t)) """ a = var('a') w = var('w') epsilon = var('epsilon') psi=var('psi') t=var('t') u=var('u') ud=var('ud') FN=[ ] for n in range(NMAX): FN.append(1/pi*integrate(func(-a*w*sin(psi))*cos(n*(psi+pi/2.0)),psi,0,2*pi)) if len(FN)<1: FN.append(0) FN[0]=1/(2*pi)*integrate(func(-a*w*sin(psi)),psi,0,2*pi) print 'Fn[i]=' show(FN) GN = [ ] HN = [] for n in range(NMAX): GN.append(FN[n]*cos(n*pi/2)) HN.append(-FN[n]*cos(n*pi/2)) print 'gn[i]=' show(GN) print 'hn[i]=' show(HN) B1=0 A1=FN[1]/(2*w) print 'B1=' show(B1) print 'A1=' show(A1) def u1(a,psi): res=0 for n in range(NMAX): if n!=1: res+=FN[n]*cos(n*(psi+pi/2))/(n^2-1) return -1/w^2*res print 'u1(a,psi)=' show(u1(a,psi)) A2=0 sum=0 for n in range(2,NMAX): if n!=1: sum += n*(FN[n])^2/(n^2-1) B2=FN[1]/(8*w^3*a)*diff(FN[1],a)-(FN[1])^2/(4*w^3*a^2)-1/(2*w^3*a^2)*sum print 'A2=' show(A2) print 'B2=' show(B2) def dadt(t): return epsilon/(2*w)*FN[1] def dpsidt(t): return w+epsilon*B1+epsilon^2*B2 def psiint(t): return dpsidt(t).integrate(t) print "Krilov-Bogolyubov secondt approximation for equation:\n d^2 u/dt+w^2*u=epsilon*f(ud);\n f(ud)=",func(ud),' ud=du/dt' print '\nRESULTS:' print 'da/dt(t)=' show(dadt(t)) print 'dpsi/dt(t)=' show(dpsidt(t)) """resulting function x(t)""" def resu(a,psi): sum=0 for n in range(2,NMAX): sum+=epsilon/w^2*FN[n]*cos(n*psi)/(n^2-1) return a*cos(psi)-epsilon*FN[0]/w^2+sum print 'x(t)=' show(resu(a,psi)) print 'where a(t) and psi(t) are obtained from equations above' return 0 
       
print "========================================" print "|| equation3: x''+2x+1/2x'^3+1/2x'^5=0 " print "========================================" print "epsilon=1/2" f(xd)=(xd^3+2/3*xd^5) print "f(x,xd)=",f(xd) show(f(xd)) krilov_bogolyubov_second_nox(f,NMAX=4) 
        
========================================
|| equation3: x''+2x+1/2x'^3+1/2x'^5=0    
========================================
epsilon=1/2
f(x,xd)= 2/3*xd^5 + xd^3

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{2}{3} \, \mbox{xd}^{5} + \mbox{xd}^{3}

Fn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, \frac{5 \, \pi a^{5} w^{5} + 9 \, \pi a^{3} w^{3}}{12 \, \pi}, 0, \frac{5 \, \pi a^{5} w^{5} + 6 \, \pi a^{3} w^{3}}{24 \, \pi}\right]

gn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, 0, 0, 0\right]

hn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, 0, 0, 0\right]

B1=

\newcommand{\Bold}[1]{\mathbf{#1}}0

A1=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{5 \, \pi a^{5} w^{5} + 9 \, \pi a^{3} w^{3}}{24 \, \pi w}

u1(a,psi)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{{\left(5 \, \pi a^{5} w^{5} + 6 \, \pi a^{3} w^{3}\right)} \cos\left(\frac{3}{2} \, \pi + 3 \, \psi\right)}{192 \, \pi w^{2}}

A2=

\newcommand{\Bold}[1]{\mathbf{#1}}0

B2=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(25 \, \pi a^{4} w^{5} + 27 \, \pi a^{2} w^{3}\right)} {\left(5 \, \pi a^{5} w^{5} + 9 \, \pi a^{3} w^{3}\right)}}{1152 \, \pi^{2} a w^{3}} - \frac{{\left(5 \, \pi a^{5} w^{5} + 6 \, \pi a^{3} w^{3}\right)}^{2}}{3072 \, \pi^{2} a^{2} w^{3}} - \frac{{\left(5 \, \pi a^{5} w^{5} + 9 \, \pi a^{3} w^{3}\right)}^{2}}{576 \, \pi^{2} a^{2} w^{3}}

Krilov-Bogolyubov secondt approximation for equation:
 d^2 u/dt+w^2*u=epsilon*f(ud);
 f(ud)= 2/3*ud^5 + ud^3   ud=du/dt

RESULTS:
da/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(5 \, \pi a^{5} w^{5} + 9 \, \pi a^{3} w^{3}\right)} \epsilon}{24 \, \pi w}

dpsi/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{9216} \, {\left(\frac{8 \, {\left(25 \, \pi a^{4} w^{5} + 27 \, \pi a^{2} w^{3}\right)} {\left(5 \, \pi a^{5} w^{5} + 9 \, \pi a^{3} w^{3}\right)}}{\pi^{2} a w^{3}} - \frac{3 \, {\left(5 \, \pi a^{5} w^{5} + 6 \, \pi a^{3} w^{3}\right)}^{2}}{\pi^{2} a^{2} w^{3}} - \frac{16 \, {\left(5 \, \pi a^{5} w^{5} + 9 \, \pi a^{3} w^{3}\right)}^{2}}{\pi^{2} a^{2} w^{3}}\right)} \epsilon^{2} + w

x(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}a \cos\left(\psi\right) + \frac{{\left(5 \, \pi a^{5} w^{5} + 6 \, \pi a^{3} w^{3}\right)} \epsilon \cos\left(3 \, \psi\right)}{192 \, \pi w^{2}}

where a(t) and psi(t) are obtained from equations above
\newcommand{\Bold}[1]{\mathbf{#1}}0
========================================
|| equation3: x''+2x+1/2x'^3+1/2x'^5=0    
========================================
epsilon=1/2
f(x,xd)= 2/3*xd^5 + xd^3

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{2}{3} \, \mbox{xd}^{5} + \mbox{xd}^{3}

Fn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, \frac{5 \, \pi a^{5} w^{5} + 9 \, \pi a^{3} w^{3}}{12 \, \pi}, 0, \frac{5 \, \pi a^{5} w^{5} + 6 \, \pi a^{3} w^{3}}{24 \, \pi}\right]

gn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, 0, 0, 0\right]

hn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, 0, 0, 0\right]

B1=

\newcommand{\Bold}[1]{\mathbf{#1}}0

A1=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{5 \, \pi a^{5} w^{5} + 9 \, \pi a^{3} w^{3}}{24 \, \pi w}

u1(a,psi)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{{\left(5 \, \pi a^{5} w^{5} + 6 \, \pi a^{3} w^{3}\right)} \cos\left(\frac{3}{2} \, \pi + 3 \, \psi\right)}{192 \, \pi w^{2}}

A2=

\newcommand{\Bold}[1]{\mathbf{#1}}0

B2=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(25 \, \pi a^{4} w^{5} + 27 \, \pi a^{2} w^{3}\right)} {\left(5 \, \pi a^{5} w^{5} + 9 \, \pi a^{3} w^{3}\right)}}{1152 \, \pi^{2} a w^{3}} - \frac{{\left(5 \, \pi a^{5} w^{5} + 6 \, \pi a^{3} w^{3}\right)}^{2}}{3072 \, \pi^{2} a^{2} w^{3}} - \frac{{\left(5 \, \pi a^{5} w^{5} + 9 \, \pi a^{3} w^{3}\right)}^{2}}{576 \, \pi^{2} a^{2} w^{3}}

Krilov-Bogolyubov secondt approximation for equation:
 d^2 u/dt+w^2*u=epsilon*f(ud);
 f(ud)= 2/3*ud^5 + ud^3   ud=du/dt

RESULTS:
da/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(5 \, \pi a^{5} w^{5} + 9 \, \pi a^{3} w^{3}\right)} \epsilon}{24 \, \pi w}

dpsi/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{9216} \, {\left(\frac{8 \, {\left(25 \, \pi a^{4} w^{5} + 27 \, \pi a^{2} w^{3}\right)} {\left(5 \, \pi a^{5} w^{5} + 9 \, \pi a^{3} w^{3}\right)}}{\pi^{2} a w^{3}} - \frac{3 \, {\left(5 \, \pi a^{5} w^{5} + 6 \, \pi a^{3} w^{3}\right)}^{2}}{\pi^{2} a^{2} w^{3}} - \frac{16 \, {\left(5 \, \pi a^{5} w^{5} + 9 \, \pi a^{3} w^{3}\right)}^{2}}{\pi^{2} a^{2} w^{3}}\right)} \epsilon^{2} + w

x(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}a \cos\left(\psi\right) + \frac{{\left(5 \, \pi a^{5} w^{5} + 6 \, \pi a^{3} w^{3}\right)} \epsilon \cos\left(3 \, \psi\right)}{192 \, \pi w^{2}}

where a(t) and psi(t) are obtained from equations above
\newcommand{\Bold}[1]{\mathbf{#1}}0
print "=================================" print "|| equation4: x''+2x+1/2x'^2=0 " print "=================================" print "epsilon=1/2" f4(xd)=-(xd)^2 print "f(x,xd)=",f4(xd) show(f4(xd)) krilov_bogolyubov_second_nox(f4,NMAX=4) 
        
=================================
|| equation4: x''+2x+1/2x'^2=0    
=================================
epsilon=1/2
f(x,xd)= -xd^2

\newcommand{\Bold}[1]{\mathbf{#1}}-\mbox{xd}^{2}

Fn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{1}{2} \, a^{2} w^{2}, 0, -\frac{1}{2} \, a^{2} w^{2}, 0\right]

gn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{1}{2} \, a^{2} w^{2}, 0, \frac{1}{2} \, a^{2} w^{2}, 0\right]

hn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[\frac{1}{2} \, a^{2} w^{2}, 0, -\frac{1}{2} \, a^{2} w^{2}, 0\right]

B1=

\newcommand{\Bold}[1]{\mathbf{#1}}0

A1=

\newcommand{\Bold}[1]{\mathbf{#1}}0

u1(a,psi)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{a^{2} w^{2} \cos\left(\pi + 2 \, \psi\right) - 3 \, a^{2} w^{2}}{6 \, w^{2}}

A2=

\newcommand{\Bold}[1]{\mathbf{#1}}0

B2=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{12} \, a^{2} w

Krilov-Bogolyubov secondt approximation for equation:
 d^2 u/dt+w^2*u=epsilon*f(ud);
 f(ud)= -ud^2   ud=du/dt

RESULTS:
da/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}0

dpsi/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{12} \, a^{2} \epsilon^{2} w + w

x(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{6} \, a^{2} \epsilon \cos\left(2 \, \psi\right) + \frac{1}{2} \, a^{2} \epsilon + a \cos\left(\psi\right)

where a(t) and psi(t) are obtained from equations above
\newcommand{\Bold}[1]{\mathbf{#1}}0
=================================
|| equation4: x''+2x+1/2x'^2=0    
=================================
epsilon=1/2
f(x,xd)= -xd^2

\newcommand{\Bold}[1]{\mathbf{#1}}-\mbox{xd}^{2}

Fn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{1}{2} \, a^{2} w^{2}, 0, -\frac{1}{2} \, a^{2} w^{2}, 0\right]

gn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{1}{2} \, a^{2} w^{2}, 0, \frac{1}{2} \, a^{2} w^{2}, 0\right]

hn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[\frac{1}{2} \, a^{2} w^{2}, 0, -\frac{1}{2} \, a^{2} w^{2}, 0\right]

B1=

\newcommand{\Bold}[1]{\mathbf{#1}}0

A1=

\newcommand{\Bold}[1]{\mathbf{#1}}0

u1(a,psi)=

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{a^{2} w^{2} \cos\left(\pi + 2 \, \psi\right) - 3 \, a^{2} w^{2}}{6 \, w^{2}}

A2=

\newcommand{\Bold}[1]{\mathbf{#1}}0

B2=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{12} \, a^{2} w

Krilov-Bogolyubov secondt approximation for equation:
 d^2 u/dt+w^2*u=epsilon*f(ud);
 f(ud)= -ud^2   ud=du/dt

RESULTS:
da/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}0

dpsi/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{12} \, a^{2} \epsilon^{2} w + w

x(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{6} \, a^{2} \epsilon \cos\left(2 \, \psi\right) + \frac{1}{2} \, a^{2} \epsilon + a \cos\left(\psi\right)

where a(t) and psi(t) are obtained from equations above
\newcommand{\Bold}[1]{\mathbf{#1}}0
def auto_oscillation_system_equation(func,NMAX): """ a function for evaluating second approximation of equation: x''+w^2*x = epsilon*f(x)*dx/dt via Krilov-Bogolyubov method. Accepts 2 arguments: \param1 : a function, function is dependent on 1 variable only, for example f(x)=x^2, \param2 Number of maximal possible fourier series approximation What is the output? It gives you these expressions da(t)/dt = epsilon*A1+epsilon^2*A2 ----> whereas a(t)=solution of this diff. eq.(you have to solve) dpsi(t)/dt = w+epsilon*B1+epsilon^2*B2 ----> whereas psi(t)=solution of this diff. eq.(you have to solve) u1(psi(t)) is calculated Finally aysmptotic aproximation of solution is given by formula: x(t)=a(t)*cos(psi(t))+u1(psi(t)) """ a = var('a') w = var('w') epsilon = var('epsilon') psi=var('psi') t=var('t') u=var('u') x=var('x') ud=var('ud') IFunc = integrate(func(x),x) FN=[ 0] for n in range(1, NMAX): FN.append(1/pi*integrate(IFunc(a*cos(psi))*cos(n*psi),psi,0,2*pi)) if len(FN)<1: FN.append(0) print 'Fn[i]=' show(FN) B1=0 A1=FN[1]/(2) print 'B1=' show(B1) print 'A1=' show(A1) def u1(a,psi): sum=0 for n in range(2,NMAX): sum+=n*FN[n]*sin(n*psi)/(n^2-1) return epsilon/w*sum print 'u1(a,psi)=' show(u1(a,psi)) A2=0 sum=0 for n in range(2,NMAX): if n!=1: sum += n^2*(FN[n])^2/(n^2-1) B2=-FN[1]/(8*w*a)*diff(FN[1],a)-1/(2*w*a^2)*sum print 'A2=' show(A2) print 'B2=' show(B2) def dadt(t): return epsilon/2*FN[1] def dpsidt(t): return w def psiint(t): return dpsidt(t).integrate(t) print "Krilov-Bogolyubov secondt approximation for equation:\n x''+w^2*x = epsilon*f(x)*dx/dt;\n f(x)=",func(x) print '\nRESULTS:' print 'da/dt(t)=' show(dadt(t)) print 'dpsi/dt(t)=' show(dpsidt(t)) """resulting function x(t)""" def resu(a,psi): sum=0 for n in range(2,NMAX): sum+=n*FN[n]*sin(n*psi)/(n^2-1) return a*cos(psi)+epsilon/w*sum print 'x(t)=' show(resu(a,psi)) print 'where a(t) and psi(t) are obtained from equations above' return 0 
       
print "================================================" print "|| equation5: x''+2x-1/2x(1-x^2-1/2*x^4)x'=0 " print "================================================" print "epsilon=1/2" f5(x)=(1-x^2-1/2*x^4) print "f(x,xd)=",f5(x) show(f5(x)) auto_oscillation_system_equation(f5,NMAX=10) 
        
================================================
|| equation5: x''+2x-1/2x(1-x^2-1/2*x^4)x'=0    
================================================
epsilon=1/2
f(x,xd)= -1/2*x^4 - x^2 + 1

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{2} \, x^{4} - x^{2} + 1

Fn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, -\frac{\pi a^{5} + 4 \, \pi a^{3} - 16 \, \pi a}{16 \, \pi}, 0, -\frac{3 \, \pi a^{5} + 8 \, \pi a^{3}}{96 \, \pi}, 0, -\frac{1}{160} \, a^{5}, 0, 0, 0, 0\right]

B1=

\newcommand{\Bold}[1]{\mathbf{#1}}0

A1=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{\pi a^{5} + 4 \, \pi a^{3} - 16 \, \pi a}{32 \, \pi}

u1(a,psi)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{{\left(a^{5} \sin\left(5 \, \psi\right) + \frac{3 \, {\left(3 \, \pi a^{5} + 8 \, \pi a^{3}\right)} \sin\left(3 \, \psi\right)}{\pi}\right)} \epsilon}{768 \, w}

A2=

\newcommand{\Bold}[1]{\mathbf{#1}}0

B2=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{a^{10} + \frac{3 \, {\left(3 \, \pi a^{5} + 8 \, \pi a^{3}\right)}^{2}}{\pi^{2}}}{49152 \, a^{2} w} + \frac{{\left(\pi a^{5} + 4 \, \pi a^{3} - 16 \, \pi a\right)} {\left(16 \, \pi - 5 \, \pi a^{4} - 12 \, \pi a^{2}\right)}}{2048 \, \pi^{2} a w}

Krilov-Bogolyubov secondt approximation for equation:
 x''+w^2*x = epsilon*f(x)*dx/dt;
 f(x)= -1/2*x^4 - x^2 + 1

RESULTS:
da/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{{\left(\pi a^{5} + 4 \, \pi a^{3} - 16 \, \pi a\right)} \epsilon}{32 \, \pi}

dpsi/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}w

x(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}a \cos\left(\psi\right) - \frac{{\left(a^{5} \sin\left(5 \, \psi\right) + \frac{3 \, {\left(3 \, \pi a^{5} + 8 \, \pi a^{3}\right)} \sin\left(3 \, \psi\right)}{\pi}\right)} \epsilon}{768 \, w}

where a(t) and psi(t) are obtained from equations above
\newcommand{\Bold}[1]{\mathbf{#1}}0
================================================
|| equation5: x''+2x-1/2x(1-x^2-1/2*x^4)x'=0    
================================================
epsilon=1/2
f(x,xd)= -1/2*x^4 - x^2 + 1

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{2} \, x^{4} - x^{2} + 1

Fn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, -\frac{\pi a^{5} + 4 \, \pi a^{3} - 16 \, \pi a}{16 \, \pi}, 0, -\frac{3 \, \pi a^{5} + 8 \, \pi a^{3}}{96 \, \pi}, 0, -\frac{1}{160} \, a^{5}, 0, 0, 0, 0\right]

B1=

\newcommand{\Bold}[1]{\mathbf{#1}}0

A1=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{\pi a^{5} + 4 \, \pi a^{3} - 16 \, \pi a}{32 \, \pi}

u1(a,psi)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{{\left(a^{5} \sin\left(5 \, \psi\right) + \frac{3 \, {\left(3 \, \pi a^{5} + 8 \, \pi a^{3}\right)} \sin\left(3 \, \psi\right)}{\pi}\right)} \epsilon}{768 \, w}

A2=

\newcommand{\Bold}[1]{\mathbf{#1}}0

B2=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{a^{10} + \frac{3 \, {\left(3 \, \pi a^{5} + 8 \, \pi a^{3}\right)}^{2}}{\pi^{2}}}{49152 \, a^{2} w} + \frac{{\left(\pi a^{5} + 4 \, \pi a^{3} - 16 \, \pi a\right)} {\left(16 \, \pi - 5 \, \pi a^{4} - 12 \, \pi a^{2}\right)}}{2048 \, \pi^{2} a w}

Krilov-Bogolyubov secondt approximation for equation:
 x''+w^2*x = epsilon*f(x)*dx/dt;
 f(x)= -1/2*x^4 - x^2 + 1

RESULTS:
da/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{{\left(\pi a^{5} + 4 \, \pi a^{3} - 16 \, \pi a\right)} \epsilon}{32 \, \pi}

dpsi/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}w

x(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}a \cos\left(\psi\right) - \frac{{\left(a^{5} \sin\left(5 \, \psi\right) + \frac{3 \, {\left(3 \, \pi a^{5} + 8 \, \pi a^{3}\right)} \sin\left(3 \, \psi\right)}{\pi}\right)} \epsilon}{768 \, w}

where a(t) and psi(t) are obtained from equations above
\newcommand{\Bold}[1]{\mathbf{#1}}0
print "================================================" print "|| equation Van der Pol: x''+x-epsilon*(1-x^2)x'=0 " print "================================================" print "epsilon=1/2" f_Pol(x)=(1-x^2) print "f(x,xd)=",f_Pol(x) show(f_Pol(x)) auto_oscillation_system_equation(f_Pol,NMAX=10) 
        
================================================
|| equation Van der Pol: x''+x-epsilon*(1-x^2)x'=0    
================================================
epsilon=1/2
f(x,xd)= -x^2 + 1

\newcommand{\Bold}[1]{\mathbf{#1}}-x^{2} + 1

Fn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, -\frac{\pi a^{3} - 4 \, \pi a}{4 \, \pi}, 0, -\frac{1}{12} \, a^{3}, 0, 0, 0, 0, 0, 0\right]

B1=

\newcommand{\Bold}[1]{\mathbf{#1}}0

A1=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{\pi a^{3} - 4 \, \pi a}{8 \, \pi}

u1(a,psi)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{a^{3} \epsilon \sin\left(3 \, \psi\right)}{32 \, w}

A2=

\newcommand{\Bold}[1]{\mathbf{#1}}0

B2=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{a^{4}}{256 \, w} + \frac{{\left(\pi a^{3} - 4 \, \pi a\right)} {\left(4 \, \pi - 3 \, \pi a^{2}\right)}}{128 \, \pi^{2} a w}

Krilov-Bogolyubov secondt approximation for equation:
 x''+w^2*x = epsilon*f(x)*dx/dt;
 f(x)= -x^2 + 1

RESULTS:
da/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{{\left(\pi a^{3} - 4 \, \pi a\right)} \epsilon}{8 \, \pi}

dpsi/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}w

x(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{a^{3} \epsilon \sin\left(3 \, \psi\right)}{32 \, w} + a \cos\left(\psi\right)

where a(t) and psi(t) are obtained from equations above
\newcommand{\Bold}[1]{\mathbf{#1}}0
================================================
|| equation Van der Pol: x''+x-epsilon*(1-x^2)x'=0    
================================================
epsilon=1/2
f(x,xd)= -x^2 + 1

\newcommand{\Bold}[1]{\mathbf{#1}}-x^{2} + 1

Fn[i]=

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, -\frac{\pi a^{3} - 4 \, \pi a}{4 \, \pi}, 0, -\frac{1}{12} \, a^{3}, 0, 0, 0, 0, 0, 0\right]

B1=

\newcommand{\Bold}[1]{\mathbf{#1}}0

A1=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{\pi a^{3} - 4 \, \pi a}{8 \, \pi}

u1(a,psi)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{a^{3} \epsilon \sin\left(3 \, \psi\right)}{32 \, w}

A2=

\newcommand{\Bold}[1]{\mathbf{#1}}0

B2=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{a^{4}}{256 \, w} + \frac{{\left(\pi a^{3} - 4 \, \pi a\right)} {\left(4 \, \pi - 3 \, \pi a^{2}\right)}}{128 \, \pi^{2} a w}

Krilov-Bogolyubov secondt approximation for equation:
 x''+w^2*x = epsilon*f(x)*dx/dt;
 f(x)= -x^2 + 1

RESULTS:
da/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{{\left(\pi a^{3} - 4 \, \pi a\right)} \epsilon}{8 \, \pi}

dpsi/dt(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}w

x(t)=

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{a^{3} \epsilon \sin\left(3 \, \psi\right)}{32 \, w} + a \cos\left(\psi\right)

where a(t) and psi(t) are obtained from equations above
\newcommand{\Bold}[1]{\mathbf{#1}}0