Bords du domaine de stabilité pour l'équation de Mathieu (méthode de Lindstedt)

Une méthode pour résoudre x''(t)+k^2x(t)=f(t) lorsque f(t) est un polynome trigoometrique de degre j, sans termes en cos(kt) ou sin(kt)

> W:=proc(f,j,k);# procedure qui resout x"+k^2*x=f ,x(0)=0, x'(0)=0
> l:=array(0..j); m:=array(1..j);S:={seq(i, i=1..j)} minus {k}; ;
> for q from 1 to j do;
> l[q]:=coeff(f,cos(q*t)); m[q]:=coeff(f,sin(q*t));od; l[0]:=f-add(l[r]*cos(r*t)+m[r]*sin(r*t), r=1..j);
> add(1/(k^2-r^2)*l[r]*(cos(r*t)-cos(k*t))+1/(k^2-r^2)*m[r]*(sin(r*t)-sin(k*t)), r=S)+1/k^2*l[0]*(1-cos(k*t));
> end proc:
 

Ces programme trouve les courbes dans l'espace (x,y) pour lesquelles il existe des solutions 2Pi périodiques.

en suivant la méthode de l'exercice (J) du chapitre 3 du polycopié

> T:=proc(k,n);
> u:=array(0..n); x:=array(0..n);u[0]:=cos(k*t); x[0]:=k^2;
> v:=sum(u[j]*y^j, j=0..n); z:=sum(x[j]*y^j, j=0..n)+y*cos(2*t);
> P:=collect(expand(v*z),trig):
> for j from 1 to n do;
> F:=combine(coeff(P,y,j)-k^2*u[j]):
> u[j]:=W(-F,2*j*k+2,k);
> C:={coeff(F,cos(k*t)), coeff(F, sin(k*t))};
> V:=solve(C,x[j]);
> subs(V, x[j]);x[j]:=%;
> od;[eval(z-y*cos(2*t)),eval(v)];
> end proc:
 
 
 

> S:=proc(k,n):
> u:=array(0..n); x:=array(0..n);u[0]:=sin(k*t); x[0]:=k^2;
> v:=sum(u[j]*y^j, j=0..n); z:=sum(x[j]*y^j, j=0..n)+y*cos(2*t);
> P:=collect(expand(v*z),trig):
> for j from 1 to n do;
> F:=combine(coeff(P,y,j)-k^2*u[j]);
> u[j]:=W(-F,2*j*k+2,k);
> C:={coeff(F,cos(k*t)), coeff(F, sin(k*t))};
> V:=solve(C,x[j]);
> subs(V, x[j]);x[j]:=%;
> od;[eval(z-y*cos(2*t)),eval(v)];
> end proc:

Enfin nous tracons les courbes, limitant le domaine de stabilite, ou plutot un developpement limité de ces courbes. Notre programme donnant  x comme fonction de y, nous effectuons une rotation pour revenir aux coordonnées usuelles:

> with(plots):
> p:=plot([S(1,3)[1],T(1,3)[1],S(2,3)[1],T(2,3)[1]],y=-3..3):
> with(plottools): rotate(p,-Pi/2);
 
Voici le résultat