Third order normal form for a center

Consider an ODE in the plane which has a center at the origin with eigenvalues i b and -i b.  Assuming that the linear part is already in real canonical form, we will show how to bring the nonlinear terms of order two and three into normal form.  If the ODE is xdot = f(x,y) and ydot = g(x,y) we will write the Taylor series in terms of the partial derivatives fxx, fxy etc.

xdot= -omega* y +fxx*x^2/2+fxy*x*y+fyy*y^2/2 + fxxx*x^3/6 + fxxy*x^2*y/2+fxyy*x*y^2/2+fyyy*y^3/6
ydot= omega*x + gxx*x^2/2+gxy*x*y+gyy*y^2/2 + gxxx*x^3/6 + gxxy*x^2*y/2+gxyy*x*y^2/2+gyyy*y^3/6

Convert the problem into complex notation to get the matrix to be diagonal.  Only the zdot equation is printed out.  The wdot equation is obtained by conjugating all the constants and switching z and w.

zdot =  xdot +I*ydot/.{x->(z+w)/2, y->(z-w)/(2*I)}//Expand
wdot =  xdot -I*ydot  /.{x->(z+w)/2, y->(z-w)/(2*I)}//Expand;

Next we find a quadratic change of variables which eliminates all the quadratic terms.  This involves solving the order two homological equation.  We set up the quadratic change of variables with undetermined coefficients
c[i,j] + I  d[i,j] where the second component is the conjugate of the first with the variables swapped.

quadraticterms = (Coefficient[zdot, z^2]/.{z->0,w->0})*z^2 + (Coefficient[zdot, z*w]/.{z->0,w->0})*z*w +(Coefficient[zdot,w^2]/.{z->0,w->0})*w^2

k2 = {Sum[(c[i,2-i]+I*d[i,2-i])*z^i*w^(2-i),{i,0,2}],Sum[(c[i,2-i]-I*d[i,2-i])*w^i*z^(2-i),{i,0,2}]}
Dk2=Map[D[k2,#]&,{z,w}]//Transpose;
A = DiagonalMatrix[{I*omega,-I*omega}];

homologicalequation = Dk2.A.{z,w}-A.k2 -quadraticterms//Simplify;
zhomeq = homologicalequation[[1]]

Solving for c[i,j] and d[i,j] gives us our quadratic transformation.

cdeqs =Map[Coefficient[zhomeq,#]&,{z^2,z*w,w^2}]//ComplexExpand;
realcdeqs = Re[cdeqs]~Join~Im[cdeqs]//ComplexExpand;
cdsol =  Solve[realcdeqs=={0,0,0,0,0,0},{c[2,0],c[1,1],c[0,2],d[2,0],d[1,1],d[0,2]}][[1]]//Simplify

Make the transformation and calculate the new ODE.  It’s a mess so none of this will be printed out.

k = {z,w}+k2/.cdsol//Simplify;
Dk = Map[D[k,#]&,{z,w}]//Transpose//Simplify;
Dkinv= Inverse[Dk]//Simplify;
newode = Dkinv.({zdot,wdot}/.{z->k[[1]],w->k[[2]]});

According to the theory, we could make a further coordinate change to eliminate all of the terms of degree three except the one involving the monomial z^2 w.  Here is its  coefficient c1 = a1+ I* b1

c1=SeriesCoefficient[newode[[1]],{z,0,2},{w,0,1}]//Simplify
a1= Re[c1]//ComplexExpand//Simplify
b1= Im[c1]//ComplexExpand//Simplify

The formulas can be conveniently split as the sum of the old coefficient for z^2 w and some new contributions coming from our quadratic change of variables.  These formulas can now be conveniently applied to analyze any center-type equilibrium.

oldc1= c1=SeriesCoefficient[zdot,{z,0,2},{w,0,1}]//Simplify;
a1old= Re[c1]//ComplexExpand//Simplify; a1new =a1-a1old//Simplify; 
b1old= Im[c1]//ComplexExpand//Simplify;b1new =b1-b1old//Simplify; 
Print["a1= ",a1old + a1new]
Print["b1= ",b1old + b1new]

Example

Here is a planar ODE with a center.  We will apply the formula to analyze the stabliity.  Consider a pendulum with an additional, nonlinear “friction” term:  xddot + 4 sin(x) +xdot^3 = 0.  If we let y =xdot be the velocity we get the two-dimensional system.

f= y;g= -4Sin[x] - y^3;
{f,g}

To get the linearized system at the origin into normal form we make a preliminary linear change of variables.
One such choice is x = x1, y = -2 y1.

P = {{1,0},{0,-2}}//Transpose
P.{x1,y1}
{f1,g1} =Inverse[P].({f,g}/.{x->x1,y->-2*y1})//Simplify

Now we can just compute the derivatives of the new f and g at the origin and subsitute into the formulas for a1, b1

sub = {fxx->D[f1,{x1,2}],fxy->D[f1,{x1,1},{y1,1}],fyy->D[f1,{y1,2}],gxx->D[g1,{x1,2}],gxy->D[g1,{x1,1},{y1,1}],gyy->D[g1,{y1,2}],
fxxx->D[f1,{x1,3}],fxxy->D[f1,{x1,2},{y1,1}],fxyy->D[f1,{x1,1},{y1,2}],fyyy->D[f1,{y1,3}],gxxx->D[g1,{x1,3}],gxxy->D[g1,{x1,2},{y1,1}],gxyy->D[g1,{x1,1},{y1,2}],gyyy->D[g1,{y1,3}],omega->2}/.{x1->0,y1->0}

a1/.sub
b1/.sub