3743.sagews

var('x t psi a w xd epsilon')

\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))

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]

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]

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

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

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=B2= <html><div class="math">\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