input_case nil$ n:=6$ % Classical nonlinear oscillator % ------------------------------ % particle with unit mass in the potential x^2/2 + V(x) % where V(x) = for i:=1:infinity sum c(i)*x^(i+2) % perturbation theory up to the order n in the amplitude a % the variable tau in the program means omega * time % where omega^2 = 1 + u. Equation of motion: % (1+u) * df(x,tau,2) + x = R % where R = - df(V(x),x) operator b,c; weight a=1$ wtlevel n+1$ % general form of the solution x:=a*cos(tau) +for i:=2:n+1 sum a^i* for j:=if evenp(i) then 0 else 3 step 2 until i sum b(i,j)*cos(j*tau)$ % expansion in harmonics R:=(-for i:=1:n sum c(i)*(i+2)*x^(i+1) where {cos(~x)^2=>(1+cos(2*x))/2, cos(~x)*cos(~y)=>(cos(x+y)+cos(x-y))/2})$ u:=(z*R where z*cos(tau)=>1)$ % separation of the resonant term u:=-(u where z=>0)/a$ % condition of its cancellation determines u R:=R-u*df(x,tau,2)$ % equation of motion is now df(x,tau,2)+x=R R:=(R where cos(~j*tau)=>z^j)$ % for separation of harmonics clear a; % small variables can't be used in coeff and factor l:=coeff(R,a)$ clear R; l:=rest(rest(l))$ % separating powers of a for i:=2:n+1 do << Ri:=first(l); l:=rest(l); if evenp(i) then j0:=0 else << j0:=3; Ri:=Ri/z^3 >>; li:=coeff(sub(z=sqrt(z),Ri),z); % separating harmonics for j:=j0 step 2 until i do << b(i,j):=first(li)/(1-j^2); li:=rest(li) >> >>; factor a; on revpri,rat; % printing results 1+u; % omega^2 factor cos; coeff(x,a); % particle motion showtime; off nat,pri,echo; out "class.res"; u:=u; x:=x; write "end"$ shut "class.res"; end;