Ingemar Nåsell, Institutionen för Matemaitk, KTH, 001123.

Fasporträttet för den ickelinjära pendeln.

> restart; with(DEtools,DEplot): with(linalg,vector): with(plots,display):

Efter skalning kan ekvationssystemet skrivas så här:

> ekv1:=D(x)(t)=y:
ekv2:=D(y)(t)=-sin(x)-2*alpha*y:

Rita riktningsfältet med :

> alpha:=0.1:

> DEplot([ekv1,ekv2],[x,y],t=0..1,x=-9..9,y=-4..4);

Jag har sadelpunkter vid de kritiska punkterna ( ) och ( ) . Jag kallar dem P1 och P2:

> P1:=vector([-evalf(Pi),0]): P2:=vector([evalf(Pi),0]):

Egenvärdena vid varje sadelpunkt är och . Det positiva egenvärdet kallas rinstab och det negativa rstab. Med får jag de numeriska värdena:

> rinstab:=-alpha+sqrt(1+alpha^2): rstab:=-alpha-sqrt(1+alpha^2):

Motsvarande egenvektorer har lutningen rinstab respektive rstab. Jag bestämmer en egenvektor med längden 1 för vart och ett av de två egenvärdena:

> eginstab:=evalm(vector([1,rinstab])/sqrt(1+rinstab^2)):

> egstab:=evalm(vector([1,rstab])/sqrt(1+rstab^2)):

Nu bestämmer jag 4 initialpunkter nära var och en av de två sadelpunkterna med riktning bestämd av egenvektorn och med avståndet 0.1 från sadelpunkten:

> initP1instab1:=evalf(evalm(P1+0.1*eginstab)):
initP1instab2:=evalf(evalm(P1-0.1*eginstab)):
initP1stab1:=evalf(evalm(P1+0.1*egstab)):
initP1stab2:=evalf(evalm(P1-0.1*egstab)):

> initP2instab1:=evalf(evalm(P2+0.1*eginstab)):
initP2instab2:=evalf(evalm(P2-0.1*eginstab)):
initP2stab1:=evalf(evalm(P2+0.1*egstab)):
initP2stab2:=evalf(evalm(P2-0.1*egstab)):

Initialpunkter för att rita fasporträtt med DEplot skall ges i formen [x(0)=x0,y(0)=y0]. Jag bildar två sekvenser av sådana initialvärden, en för de fyra initialvärden som behövs för att rita banorna på de instabila mångfalderna, och en för de fyra initialvärden som behövs för att rita banorna på de stabila mångfalderna.

> initinstab:=[x(0)=initP1instab1[1], y(0)=initP1instab1[2]],
[x(0)=initP1instab2[1], y(0)=initP1instab2[2]],
[x(0)=initP2instab1[1], y(0)=initP2instab1[2]],
[x(0)=initP2instab2[1], y(0)=initP2instab2[2]]:

> initstab:=[x(0)=initP1stab1[1], y(0)=initP1stab1[2]],
[x(0)=initP1stab2[1], y(0)=initP1stab2[2]],
[x(0)=initP2stab1[1], y(0)=initP2stab1[2]],
[x(0)=initP2stab2[1], y(0)=initP2stab2[2]]:

Nu bildar jag två plotstrukturer, en som visar de 4 banorna på de instabila mångfalderna, och en med de 4 banorna på de stabila mångfalderna:

> plotinstab:=DEplot([ekv1,ekv2],[x,y],t=0..5,[initinstab], x=-9..9, y=-4..4, arrows=none, linecolor=green, stepsize=0.1):

> plotstab:=DEplot([ekv1,ekv2],[x,y],t=-10..0,[initstab], x=-9..9, y=-4..4, arrows=none, linecolor=red, stepsize=0.1):

Lägg in ytterligare 3 banor med startpunkter för x=-10:

> init1:=[x(0)=-10,y(0)=1.6]: init2:=[x(0)=-10,y(0)=3.0]: init3:=[x(0)=-10,y(0)=4.4]:

> plot3:=DEplot([ekv1,ekv2],[x,y],t=0..20,[init1,init2,init3], x=-9..9, y=-4..4, arrows=none, stepsize=0.1, obsrange=false):

> display({plotinstab,plotstab,plot3}, thickness=3, title="Dämpad olinjär pendel med alfa=0.1");