Ingemar Nåsell, Matematiska Institutionen, KTH, 001029.

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

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

DATORLABORATION 2:

Fasporträtt för linjära system.

I denna labb skall vi använda Maple för att rita fasporträtt för linjära system av diffekvationer i planet. Ekvationssystemet kan skrivas u' = Mu, där M en given 2x2-matris och u är en kolumnvektor med komponenterna x och y. Vi begränsar oss till ekvationssystem för vilka origo är en stabil egentlig nod eller en sadelpunkt. Dessa två typer av fasporträtt har en egenskap gemensam: matrisen M har två linjärt oberoende reella egenvektorer. Dessa två egenvektorer definierar två stycken s.k. mångfalder som är betydelsefulla för fasporträttet. För de linjära systemen är mångfalderna helt enkelt två räta linjer genom origo vars lutningar bestäms av egenvektorerna. För ickelinjära system är mångfalderna vanligtvis mera komplicerade kurvor genom de kritiska punkterna.

De namngivna mångfalderna är de viktigaste banorna i fasporträttet. Det är lätt att rita dem för ett linjärt system eftersom de är räta linjer. Den metod vi skall använda baserar sig dock inte på linjariteten. Den är litet mer komplicerad, men har den stora fördelen att den fungerar också för olinjära system.

Metoden går ut på att för varje mångfald låta Maple rita de två banor som tillsammans med den kritiska punkten (origo) utgör mångfalden ifråga. De två banorna ritas genom att välja två initialpunkter, en på var sin sida om den kritiska punkten. Initialpunkten väljs nära den kritiska punkten för alla mångfalder utom den stabila nodens långsamma mångfald, där initialpunkten väljs långt från den kritiska punkten. Riktningen från origo till initialpunkten bestäms av lutningen hos motsvarande egenvektor. Reglerna kan sammanfattas i följande tabell. Med avstånd menas här avstånd från den kritiska punkten till banornas initialpunkter.

Fasporträtt Mångfald Tidsriktning Avstånd

Sadelpunkt Instabil Framåt Litet

Sadelpunkt Stabil Bakåt Litet

Stabil eg. nod Snabb Bakåt Litet

Stabil eg. nod Långsam Framåt Stort

Notera att om du väljer en initialpunkt på sadelpunktens stabila mångfald och låter tiden gå framåt så kommer banan att närma sig origo när tiden växer. Om initialpunkten ligger mycket nära origo kommer du därmed att rita endast en mycket liten del av motsvarande bana. För att du skall kunna se varifrån denna bana kommer är det nödvändigt att låta tiden gå baklänges. Samma resonemang gäller för banorna på den stabila egentliga nodens snabba mångfald. Ytterligare kommentarer till reglerna i tabellen ovan ges i slutet av labben.

EXEMPEL 1 (stabil egentlig nod).

Börja med att nollställa Maple och att läsa in de kommandon i minnet som vi behöver i fortsättningen. Från paketet DEtools använder vi kommandot DEplot , från paketet linalg kommandona matrix , vector och eigenvects , samt från paketet plots kommandot display . Notera att om du trycker på RETURN efter att du har givit ett kommando till Maple så exekveras kommandot. Här skriver vi fyra kommandon i samma grupp. Detta åstadkommer du genom att trycka på SHIFT + RETURN efter var och en av de tre första raderna. När du sedan trycker på RETURN så exekveras alla tre kommandona samtidigt.

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

Vi ritar fasporträttet för det linjära ekvationssystemet u' = M1*u, där huvuddiagonalen i matrisen M1 har elementen -1 och -3, och de övriga två elementen är lika med 1. Matrisen M1 matas in som en lista av listor:

> M1:=matrix( [ [ -1, 1], [ 1, -3] ] );

Vi matar nu in vektorn u med komponenterna x och y:

> u:=vector([x, y]);

Produkten av matrisen M1 med kolumnvektorn u beräknas och benämnes M1u . Matrismultiplikation betecknas med &*. Kommandot evalm (evaluate matrix) behövs för att bestämma denna produkt.

> M1u:=evalm(M1&*u);

Vi bildar nu listan av de två differentialekvationer som bestäms av att derivatan av u med avseende på t är lika med M1u . Listan benämns ekv1 . Som alltid i Maple skrivs de två elementen i listan innanför hakparenteser.

> ekv1:=[D(x)(t)=M1u[1], D(y)(t)=M1u[2]];

Kommandot DEplot är vårt huvudkommando för fortsättningen. Det tar tre (eller fler) argument. Det första argumentet specificerar systemet av diffekvationer, det andra ger en lista med namnen på variablerna, det tredje ger ett intervall av t-värden, och det fjärde ger en lista av initialvärden. Varje initialvärde skall ges i formen [x(t0)=x0, y(t0)=y0]. Resultatet av kommandot är att Maple ritar både ett riktningsfält och ett fasporträtt som innehåller de mot de givna initialvärdena svarande banorna. Om inga initialvärden ges så ritar Maple bara riktningsfältet. Ytterligare information om detta kommando finns under rubriken Help.

För att visa hur detta fungerar börjar vi med ett enda initialvärde [x(0)=1,y(0)=1]. En lista med detta initialvärde bildas genom att innesluta initialvärdet i hakparenteser. Vi lägger till två kommandon som bestämmer storleken av det fönster i vilket Maple ritar.

> init:=[[x(0)=1, y(0)=1]];
plot1:=DEplot( ekv1, [x,y], t=0..5, init, x=-1..1, y=-1..1 ):
plot1;

Vi ser att kommandot ritar både riktningsfältet och den bana som startar i punkten (1,1). Det är en smaksak om man vill ha med riktningsfältet eller inte. Vi väljer att rita fasporträtt utan riktningsfält i fortsättningen. Detta åstadkoms med tilläggskommandot arrows=none. Vi ser vidare att banan har en knyck i början. För att få en glattare bana minskar vi steglängden. Notera att du kan förenkla ditt skrivarbete genom att hämta kommandot på raden härovan genom att först markera och kopiera det och sedan klistra in det här nedan och slutligen lägga till de två sista argumenten.

> plot2:=DEplot( ekv1, [x,y], t=0..5, init, x=-1..1, y=-1..1, arrows=none, stepsize=0.1 ):
plot2;

För att bestämma typen av fasporträtt och för att rita banorna på de viktiga mångfalderna beräknar vi först matrisens egenvärden och egenvektorer. Detta intressanta arbete överlåter vi till Maple.

> eg1:=eigenvects(M1);

Numeriskt evaluering ger:

> evalf(eg1);

Detta visar att M1 har två distinkta negativa egenvärden, båda med algebraiska multipliciteten 1. Härav drar vi slutsatsen att origo är en stabil egentlig nod. Variabeln eg1 består av en sekvens av två listor. Absolutbeloppet av egenvärdet är minst i den första listan. Detta betyder att den första listan är associerad med den långsamma mångfalden och den andra listan med den snabba mångfalden. (Om ordningsföljden i dina beräkningar skulle vara omkastad så får du ändra benämningarna i de följande beräkningarna i enlighet med detta.) Vi ser att egenvektorn som hör till det första egenvärdet har positiv lutning. Detta betyder att den långsamma mångfalden är en rät linje genom origo med positiv lutning. Egenvektorn som hör till det andra egenvärdet har negativ lutning. Alltså är lutningen på den snabba mångfalden negativ.

Vi använder oss av uttrycket för eg1 för att bestämma initialvärden när vi ritar de två banorna på var och en av de två mångfalderna. För varje mångfald väljer vi två initialpunkter med koordinaterna (x0,y0) och (-x0,-y0). Vi väljer vektorn (x0,y0) så att den får samma lutning som motsvarande egenvektor och så att dess längd blir d (som vi sätter till ett lämpligt värde). Som vi sett ovan skall initialvärdena ges i formen [x(t0)=x0,y(t0)=y0], där t0 = 0. Vi skriver en liten procedur som bestämmer mängden av två initialvärden genom att plocka ut den nödvändiga informationen från motsvarande lista i eg1 . Vi benämner proceduren eg2init0 (på engelska är ordet 2=two homonymt med to; tänk på procedurnamnet som "from eg to init"). Nollan i slutet på procedurens namn är ditsatt för att namnet inte skall kollidera med namnet på en likartad procedur i Laboration 3. Indata till proceduren är den komponent (eg) av Maples svar på kommandot eigenvects som associeras med den mångfald som vi studerar, samt avståndet (d) från var och en av de kritiska punkterna och origo. I proceduren bestäms först u som den egenvektor associerad med mångfalden som finns i eg. Därefter beräknas längden lu av denna egenvektor. Sedan bestäms en egenvektor v som har längden d. Slutligen bestäms en sekvens av två listor av initialvärden, vart och ett i formen [x(0)=x0,y(0)=y0].

> eg2init0:=proc(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)=(-1)^k*v[1],y(0)=(-1)^k*v[2]], k=1..2)]
end;

Med hjälp av proceduren eg2init0 bestämmer vi nu de två listorna av initialvärden på de två mångfalderna. Vi tar initialpunkterna långt från origo (d=2) på den långsamma mångfalden och nära origo (d=0.05) på den snabba mångfalden. Här måste du använda observationen ovan att den långsamma mångfalden associeras med eg1[1] och att den snabba mångfalden associeras med eg1[2].

> `initlångsam`:=eg2init0(eg1[1], 2);
initsnabb:=eg2init0(eg1[2], 0.05);

De två banorna på den långsamma mångfalden ritas först. Man kan få pröva sig fram till en lämplig längd på tidsintervallet. Tiden skall gå framåt. Det första kommandot nedan definierar en plotstruktur, medan kommandot på den andra raden beordrar plottning. Avsluta det första kommandot med ett kolon för att förhindra utskrift.

> `plotlångsam`:=DEplot( ekv1, [x,y], t=0..5, `initlångsam`, x=-1..1, y=-1..1, arrows=none ):
`plotlångsam`;

Som väntat ser vi två linjesegment med positiv lutning. De två banorna på den snabba mångfalden ritas genom att låta tiden gå baklänges. Här går det fortare (varför?) så det räcker med en kortare integrationstid.

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

Som väntat ser vi här två linjesegment med negativ lutning.

Med hjälp av kommandot display ritar vi alla 4 banorna samtidigt:

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

Vi kompletterar nu fasporträttet med att rita ytterligare ett antal banor, förslagsvis 12 eller 16. För att få ett estetiskt tillfredsställande resultat kan man få experimentera med olika initialvärden. Vi väljer 12 initialpunkter fördelade på följande fyra sekvenser av punkter:

> inita:= [x(0)=1,y(0)=1], [x(0)=0.5,y(0)=1], [x(0)=0,y(0)=1];
initb:=[x(0)=-1,y(0)=-1], [x(0)=-0.5,y(0)=-1], [x(0)=0,y(0)=-1];
initc:=[x(0)=1,y(0)=-1], [x(0)=1,y(0)=-0.5], [x(0)=1,y(0)=0];
initd:=[x(0)=-1,y(0)=1], [x(0)=-1,y(0)=0.5], [x(0)=-1,y(0)=0];

Vi bildar en lista som innehåller dessa tolv punkter. Genom att avsluta med ett kolon i stället för ett semikolon förhindras utskrift.

> init12:=[inita, initb, initc, initd]:

Nu kan vi rita de 12 banor som motsvarar dessa 12 initialpunkter:

> plot12:=DEplot( ekv1, [x,y], t=0..5, init12, stepsize=0.1, arrows=none ):
plot12;

Nu återstår endast att rita hela fasporträttet med de 4 banor på de snabba och långsamma mångfalderna som vi ritade först och de 12 banorna ovan. Slutresultatet får vi genom att använda kommandot display . När du senare gör en utskrift så är det lämpligt att du ersätter Mitt namn i kommandot nedan med ditt namn, så att du kan skilja din utskrift från dina kamraters. Namnet i kommandot omges av bakåtlutande apostrofer.

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

Notera den goda överensstämmelsen mellan de två uppsättningarna banor! Inga banor korsar varandra.

Notera också huvudegenskapen hos den stabila noden: Alla banor som inte startar på den snabba mångfalden närmar sig origo i riktning av den långsamma mångfalden.

Det sista steget i produktionen av fasporträttet består i att göra en utskrift och att för hand förse banorna med pilar som visar i vilken riktning man rör sig på respektive bana när tiden växer.

EXEMPEL 2 (sadelpunkt).

I detta exempel skall vi plotta fasporträttet för ett linjärt system av diffekvationer u' = M2*u, där matrisen M2 ser ut som följer:

> M2:=matrix( [[-1, 4], [4, -3]] );

I analogi med behandlingen i Exempel 1 inför vi vektorn u, beräknar u' = M2*u, och inför ekv2 som beteckning för ekvationssystemet:

> u:=vector([x, y]);
M2u:=evalm(M2&*u);
ekv2:=[D(x)(t)=M2u[1], D(y)(t)=M2u[2]];

Vi börjar med att bestämma egenvärden och egenvektorer för matrisen M2:

> eg2:=eigenvects(M2);

> evalf(eg2);

Vi drar slutsatsen att matrisen M2 har ett positivt och ett negativt egenvärde. Origo är alltså en sadelpunkt för vårt ekvationssystem. Den första listan i eg2 är associerad med den instabila mångfalden eftersom egenvärdet i den listan är positivt. Den andra listan är associerad med den stabila mångfalden. Detta betyder att eg2[1] innehåller informationen om den instabila mångfalden och att eg2[2] innehåller informationen om den stabila mångfalden. Ordningsföljden är inte garanterad, utan kan vara omkastad. Om den är det för dig så får du göra den nödvändiga omnumreringen. Precis som i Exempel 1 använder vi proceduren eg2init0 för att bestämma de mot de två mångfalderna svarande mängderna av initialvärden. I båda fallen väljer vi initialpunkter nära origo (d=0.05).

> initinstabil:=eg2init0(eg2[1], 0.05);
initstabil:=eg2init0(eg2[2], 0.05);

Banorna på den instabila mångfalden ritas genom att låta tiden gå framåt. Som förut får man experimentera med längden av tidsintervallet.

> plotinstabil:=DEplot( ekv2, [x,y], t=0..2, initinstabil, x=-1..1, y=-1..1, arrows=none ):
plotinstabil;

Banorna på den stabila mångfalden närmar sig origo när tiden växer. För att se dem från våra initialpunkter nära origo låter vi tiden gå baklänges.

> plotstabil:=DEplot( ekv2, [x,y], t=-0.5..0, initstabil, x=-1..1, y=-1..1, arrows=none ):
plotstabil;

Kan du förklara varför det räcker med ett kortare tidsintervall här?

Vi kompletterar fasporträttet med ytterligare cirka 12 banor. Här gäller det att välja initialvärden förhållandevis nära den stabila mångfalden på omkretsen av vårt fönster. Man får experimentera med olika val av initialvärden tills man blir nöjd.

> inita:=[x(0)=1,y(0)=-1], [x(0)=1,y(0)=-0.8], [x(0)=1,y(0)=-0.6];
initb:=[x(0)=-1,y(0)=1], [x(0)=-1,y(0)=0.8], [x(0)=-1,y(0)=0.6];
initc:=[x(0)=0.6,y(0)=-1], [x(0)=0.4,y(0)=-1], [x(0)=0.2,y(0)=-1];
initd:=[x(0)=-0.6,y(0)=1], [x(0)=-0.4,y(0)=1], [x(0)=-0.2,y(0)=1];
init2:=[inita, initb, initc, initd]:

> plot12s:=DEplot( ekv2, [x,y], t=0..1, init2, x=-1..1, y=-1..1, arrows=none ):
plot12s;

Nu återstår finalen, som består i att använda kommandot display för att rita hela fasporträttet:

> display({plotstabil, plotinstabil, plot12s}, title=`Mitt namn`);

Det slutgiltiga fasporträttet visar konsistens på det viset att inga banor korsar varandra.

Gör en utskrift av figuren. Hämta resultatet vid skrivaren och förse banorna med pilar som visar i vilken riktning man rör sig längs varje bana när tiden ökar!

Tolkning av fasporträtt: Om högerledet i vårt ekvationssystem kan skrivas som en gradient så kan man tolka banorna i fasporträttet som kartbilder av de vägar som vattnet kommer att rinna efter när det regnar på ett bergslandskap. En beskrivning av detta bergslandskap kan hjälpa till att förklara varför de kombinationer av initialvärden och tidsaxelriktningar som vi har valt ger oss banorna på de olika mångfalderna. Förklaringen fungerar kvalitativt även om högerledet inte är en gradient.

De två banorna på den snabba mångfalden vid en stabil egentlig nod ger kartbilderna av två små bäckar som rinner snabbt utför varsin brant lutande bergssida rakt mot origo. Banorna på den långsamma mångfalden motsvaras å andra sidan av två långsamt rinnande floder som också är på väg mot origo. Alla andra banor än dessa fyra motsvaras av bäckar som först snabbt närmar sig en av de två floderna och sedan så småningom närmar sig denna. (För att analogin med vattenströmmarna skall fungera måste vi göra det orealistiska antagandet att allt vatten tappas av i origo.)

En sadelpunkt enligt vår definition motsvaras av en sadelpunkt i bergslandskapet. De två banorna på den stabila mångfalden motsvaras av två bäckar som rinner rakt mot origo. Dessa bäckar utgör samtidigt vattendelare i bergslandskapet. De två banorna på den instabila mångfalden motsvaras av två floder som rinner åt var sitt håll bort från origo. Alla andra banor än dessa fyra motsvaras av bäckar som så småningom närmar sig en av floderna.