Ingemar Nåsell, Institutionen för Matematik, KTH, 001029.

KURS I DIFFERENTIALEKVATIONER (DIFF OCH TRANS II, DEL 1).

OBS! Denna labb är skriven för Maple 6.

DATORLABORATION 3:

Fasporträtt för olinjära system.

I denna labb skall vi använda Maple för att rita fasporträttet för ett olinjärt system av diffekvationer i planet. Grundidéerna har redan behandlats i Datorlabb 2, som handlar om fasporträtt för linjära system.

Uppgiften här är att rita fasporträttet för följande ekvationssystem:

x' = y - x - xy,

y' = 1 + x - y.

Arbetsgången ser ut så här: Vi börjar med att bestämma de kritiska punkterna. Sedan linjariserar vi systemet kring var och en av de kritiska punkterna. Detta kräver att vi beräknar den jakobianska matrisen i varje kritisk punkt. Egenvärden och egenvektorer hos jakobianen bestämmer sedan typen av fasporträtt för det linjariserade systemet vid varje kritisk punkt, samt lutningarna av eventuella namngivna mångfalder. Under förutsättning att Linjariseringssatsen kan tillämpas gäller dessa resultat också för det givna icke-linjära systemet. Regler för att rita de s.k. namngivna mångfalderna vid varje kritisk punkt som har sådana gavs i Laboration nummer 2. Vi använder dessa regler här för att rita stabila och instabila mångfalder vid varje sadelpunkt, samt snabba och ev. också långsamma mångfalder vid varje nod. Den stomme i fasporträttet som vi får på detta sätt kompletteras sedan med ytterligare ett antal banor. Vi använder Maple för alla beräkningar.

Börja med att nollställa Maple och läsa in följande kommandon i minnet:

> restart;
with(DEtools,DEplot);
with(linalg,jacobian,eigenvects);
with(plots,display);

Högerleden i de två ekvationerna benämns F och G. För fortsättningen är det lämpligt att skriva någon konstant i flyttalsform.

> F:=y-x-x*y;

> G:=1.0+x-y;

Listan av dessa två högerled benämns f och mängden av diffekvationer benämns ekv:

> f:=[F, G];
ekv:={D(x)(t)=F,D(y)(t)=G};

Vi bestämmer först isoklinerna där derivatorna är lika med noll:

> isoklin1:=solve(F, y);
isoklin2:=solve(G, y);

Vi ritar isoklinerna i ett fönster vars storlek ev. får justeras senare. Det är viktigt att fönstret är så stort att alla kritiska punkter får plats.

> plot({isoklin1,isoklin2}, x=-3..3, y=-3..3, discont=true);

Vi ser att det finns två skärningspunkter mellan isoklinerna. Detta betyder att det

finns två kritiska punkter. Vi bestämmer dem:

> krit:=solve( {F, G}, {x, y} );

Vi använder P1 för att beteckna den kritiska punkt som ligger längst till vänster och har negativa koordinater, och P2 för att beteckna den andra punkten. Notera att Maple är inte konsekvent i att ordna de två punkterna. Den första komponenten av krit behöver inte vara P1. Om ordningen skulle vara omkastad i ditt fall så får du ändra numreringen i följande kommandon.

Det finns en svaghet med Maples representation av de kritiska punkterna, och det är att var och en av dem representeras som en mängd, med koordinaterna givna inom krullparenteser. Nu är ordningsföljden av elementen i en mängd oväsentlig. Detta är oacceptabelt för oss, eftersom vi måste kunna referera till x- och y-kordinaterna för de kritiska punkterna. Vi reparerar detta i nästa steg. Den preliminära representationen av de kritiska punkterna, i mängd-form, kallas P1prel och P2prel:

> P1prel:=krit[1];
P2prel:=krit[2];

Nu skriver vi om dessa mängder som listor, med x-koordinaten först. Observera att du måste läsa av resultaten ovan för att bestämma vilken komponent som är x-koordinaten. Notera att den första komponenten av P1prel, P1prel[1], är en ekvation där y-koordinaten för P1 står i högerledet. Vi plockar ut denna koordinat med kommandot rhs (right-hand side).

> P1:=[rhs(P1prel[2]), rhs(P1prel[1])];
P2:=[rhs(P2prel[1]), rhs(P2prel[2])];

Vi beräknar nu jakobianen av högerledet f till det givna ekvationssystemet:

> J:=jacobian( f, [x,y] );

Bestäm jakobianen i de två kritiska punkterna genom att sätta (subs=substituera) koordinaterna för respektive punkt i uttrycket för jakobianen:

> J1:=subs(P1prel, eval(J) );
J2:=subs(P2prel, eval(J) );

Egenvärden och egenvektorer till dessa två matriser bestäms med kommandot eigenvects :

> eg1:=eigenvects(J1);
eg2:=eigenvects(J2);

Vi ser att jakobianen i den kritiska punkten P1 har egenvärdena -2.34 och 0.96 och att jakobianen i den kritiska punkten P2 har egenvärdena -0.79 och -2.82. Bägge de kritiska punkterna är alltså hyperboliska, och eftersom högerledet i det givna diffekvationssystemet har kontinuerliga derivator av alla ordningar så kan Linjariseringssatsen tillämpas. (Det räcker att högerledet har kontinuerliga derivator av ordning 2.) Vi drar slutsatsen att P1 är en sadelpunkt och att P2 är en stabil egentlig nod.

Associationen mellan komponenter av eg1 och eg2 och de namngivna mångfalderna vid de kritiska punkterna behövs för fortsättningen. Det är klokt att göra litet bokföring och med PoP notera följande:

eg1[1] associeras med den stabila mångfalden vid sadelpunkten P1.

eg1[2] associeras med den instabila mångfalden vid sadelpunkten P1.

eg2[1] associeras med den långsamma mångfalden vid den stabila egentliga noden P2.

eg2[2] associeras med den snabba mångfalden vid den stabila egentliga noden P2.

Om ordningsföljden mellan uppräkningen av egenvärdena skulle vara omkastad så får du ta hänsyn till detta genom lämplig omnumrering.

Dessa resultat skall användas när vi ritar banor på respektive mångfalder. Problemet har redan behandlats i en enklare tappning i Laboration 2. Enda skillnaden är att den kritiska punkten för det linjära problemet i Laboration 2 alltid var origo, medan vi här har två kritiska punkter med från noll skilda koordinater. För att få Maple att rita banor på de namngivna mångfalderna behöver vi bestämma initialpunkter som ligger nära den kritiska punkten på en rät linje genom den kritiska punkten med samma lutning som motsvarande egenvektor. För var och en av mångfalderna (utom den långsamma) bestämmer vi två sådana initialpunkter, en på vardera sidan om den kritiska punkten. Koordinaterna för dessa initialpunkter kan skrivas (x0,y0) = (Px,Py) + d*(xeg,yeg) och (x0,y0) = (Px,Py) - d*(xeg,yeg), där (Px,Py) är den kritiska punktens koordinater och (xeg,yeg) är en vektor av längden 1 med samma lutning som egenvektorn. Avståndet mellan den kritiska punkten och var och en av initialpunkterna är lika med d. Initialvärdena skall ges som en lista i formen [x(0)=x0,y(0)=y0].

Jobbet att bestämma en mängd av två sådana initialvärden för varje mångfald vid varje kritisk punkt utförs av proceduren eg2init här nedan. Den är mera generell än proceduren eg2init0 i Laboration 2, eftersom den inte förutsätter att origo är den kritiska punkten. Indata till proceduren består av den kritiska punktens koordinater (P), den komponent (eg) av Maples svar ovan på kommandot eigenvects som svarar mot den mångfald som studeras, samt avståndet (d) mellan den kritiska punkten och var och en av de två initialpunkterna. I proceduren bestäms u först som en egenvektor som associeras med den mångfald som vi studerar. Dess längd lu beräknas sedan. Därefter bestäms v som en egenvektor av längden d. I det sista steget i proceduren beräknas en sekvens av två listor av initialvärden. Varje lista ges i formen [x(0)=x0,y(0)=y0].

> eg2init:=proc(P,eg,d)
local u,lu,v,k;
u:=eg[3][1];
lu:=sqrt(u[1]^2+u[2]^2);
v:=evalf(evalm(d*u/lu));
[seq([x(0)=P[1]+(-1)^k*v[1],y(0)=P[2]+(-1)^k*v[2]],
k=1..2)];
end;

Vi använder nu proceduren eg2init för att bestämma initialvärden för banor på de tre mångfalder som är aktuella. Här använder vi oss av associationen som vi noterat ovan mellan komponenter av eg1 och eg2 och motsvarande mångfalder: eg1[1] svarar mot stabil mångfald, eg1[2] svarar mot instabil mångfald, eg2[2] svarar mot snabb mångfald. Vi väljer d=0.1.

> initstabil:=eg2init(P1, eg1[1], 0.1);
initinstabil:=eg2init(P1, eg1[2], 0.1);
initsnabb:=eg2init(P2, eg2[2], 0.1);

Vi börjar med att rita de två banorna på den instabila mångfalden vid sadelpunkten P1. Man får experimentera med tredje argumentet i kommandot DEplot som ger det tidsintervall över vilket banorna beräknas. Dessutom kan man eventuellt behöva ändra steglängden från standardvärdet (tf-t0)/20, där t0..tf är intervallet av t-värden. Vi definierar först en plotvariabel plotinstabil , och ritar sedan genom att ge kommandot plotinstabil . Notera att definitionen av plotvariablen bör avslutas med kolon och inte semikolon. Vi avstår från att rita riktningsfältet.

> plotinstabil:=DEplot( ekv, [x,y], t=0..8, initinstabil, x=-3..3, y=-3..3, arrows=none):
plotinstabil;

Övertyga dig om att banorna utgår från punkten P1 med den lutning som du kan se från eg1[1]! Notera att den bana som går uppåt och åt höger närmar sig P2 när tiden växer. Detta betyder att denna bana på sadelpunktens instabila mångfald sammanfaller med en av banorna på den stabila nodens långsamma mångfald. En sådan bana kallas heteroklinisk, eftersom den sammanbinder två kritiska punkter.

Nästa steg är att rita banorna på den stabila mångfalden vid sadelpunkten. Vi låter nu tiden gå baklänges. Tidsintervallet kan göras kortare. Det går snabbare här än på den heterokliniska banan, som ju är den som dikterade tidsintervallet för banorna på den instabila mångfalden. Som ovan bildar vi först en plotvariabel som vi sedan gör en utskrift av.

> plotstabil:=DEplot( ekv, [x,y], t=-2..0, initstabil, x=-3..3, y=-3..3, arrows=none ):
plotstabil;

Banorna på den snabba mångfalden vid den stabila egentliga noden P2 ritar vi genom att låta tiden gå baklänges. Tidsintervallets längd får man bestämma genom att experimentera med olika värden.

> plotsnabb:=DEplot( ekv, [x,y], t=-1.9..0, initsnabb, x=-3..3, y=-3..3, arrows=none ):
plotsnabb;

Som nämnts ovan har vi redan ritat den bana på den långsamma mångfalden som närmar sig P2 nerifrån och från vänster. Den återstående banan på den egentliga nodens långsamma mångfald ritar vi genom att ta en initialpunkt långt bort från P2 i den första kvadranten. Vi öppnar fönstret så att initialpunkten syns. Eftersom fönstret inte har samma storlek som för de stabila och instabila mångfalderna är detta en preliminär version.

> `plotlångsamprel`:=DEplot( ekv, [x,y], t=0..20, {[x(0)=6,y(0)=6]}, x=-3..6, y=-3..6, stepsize=0.1, arrows=none ):
`plotlångsamprel`;

Detta stämmer med vad vi kan förvänta oss. Vi ser att denna bana närmar sig noden i P2 med den långsamma mångfaldens lutning. Vi ritar den del av banan som faller innanför det fönster som vi har använt tidigare. Vi skriver obsrange=false för att få Maple att rita banan även fast den startat utanför fönstret.

> `plotlångsam`:=DEplot( ekv, [x,y], t=0..20, {[x(0)=6,y(0)=6]}, x=-3..3, y=-3..3, stepsize=0.1, arrows=none, obsrange=false ):
`plotlångsam`;

Vi bildar mängden av plotvariabler och ritar de banor vi har:

> display( {plotinstabil, plotstabil, plotsnabb, `plotlångsam`} );

De sju banorna här bildar stommen till fasporträttet. Vi noterar att den stabila mångfalden vid sadelpunkten P1 ger viktig kvalitativ information. Alla initialpunkter som ligger ovanför den ger banor som närmar sig P2 när tiden växer, medan alla initialpunkter under den ger banor som avlägsnar sig från vårt fönster. Man säger att den stabila mångfalden vid P1 ger en avgränsning av attraktionsdomänen för P2.

Vi kompletterar denna stomme till fasporträttet med ytterligare några banor. För att kunna välja lämpliga initialvärden har vi nytta av att se riktningsfältet:

> rikt:=DEplot( ekv, [x,y], t=0..1, x=-3..3, y=-3..3):
rikt;

Vi väljer fem initialvärden för fem kompletterande banor och stoppar in dem i följande lista:

> init5:=[[x(0)=-3,y(0)=1.5], [x(0)=-1,y(0)=-3], [x(0)=-1,y(0)=3], [x(0)=1.5,y(0)=-3], [x(0)=3,y(0)=2]];

Vi ritar nu dessa fem banor i samma figur.

> plot5:=DEplot( ekv, [x,y], t=0..4, init5, x=-3..3, y=-3..3, stepsize=0.1, arrows=none ):
plot5;

Finalen består i att rita alla banor vi hittills sett i samma figur. Ersätt Mitt namn med Ditt så att du kan skilja ditt resultat från dina kamraters vid skrivaren.

> display({plotinstabil,plotstabil,plotsnabb,`plotlångsam`,plot5}, title=`Mitt namn`);

Skicka slutresultatet till utskrift, hämta det vid skrivaren, och markera för hand med pilar den riktning i vilken rörelsen på de olika banorna sker när tiden ökar.