L9maple

Luento ti 15.10.02

                          Diskreetit dynaamiset systeemit.

Alustukset

>    restart:

>    with(linalg):
with(LinearAlgebra): 

Warning, the protected names norm and trace have been redefined and unprotected




Warning, the assigned name GramSchmidt now has a global binding

>    alias(rref=ReducedRowEchelonForm): alias(ref=GaussianElimination):  alias(Diag=DiagonalMatrix,Id=IdentityMatrix):

>    with(plots): 

Warning, the name changecoords has been redefined

>    #plotdots:=data->display([plot(data,style=point,symbol=circle,symbolsize=17,symbol=circle),plot(data,color=blue)]);

>    plotdots:=data->display([plot(data,style=point),plot(data,color=blue)]);

plotdots := proc (data) options operator, arrow; display([plot(data,style = point), plot(data,color = blue)]) end proc

>    # Esim:

>    #plotdots([[0,0],[1,1],[1,0],[3,2]]);

>    traje:=proc(A,x0,n)
local x,k;
  x[0]:=x0:
  for k to n do x[k]:=A.x[k-1]: end do:
plotdots(map(convert,[seq(x[k],k=0..n)],list));
end:

>    trajedot:=proc(A,x0,n)
local x,k;
  x[0]:=x0:
  for k to n do x[k]:=A.x[k-1]: end do:
plot(map(convert,[seq(x[k],k=0..n)],list),style=point);
end:

>    #traje(A,<-3,3>,10);

 1.  Yleistä

Lay 5.6

      x[k+1] = A*x[k]

>   

>   

>   

 2. Ratkaisujen graafinen kuvaaminen

Ohjelmankehitystä ja Esim. 2

Jos et halua katsella ohjelmankehitystä, klikkaa itsesi tänne.

Katsotaan, miten trajektoreita olisi kätevää piirtää.

>    A:=<<.8,0>|<0,.64>>;

A := Matrix(%id = 4364152)

>    lambda[1]:=0.8; lambda[2]:=0.64;

lambda[1] := .8

lambda[2] := .64

>    v[1]:=<1,0>; v[2]:=<0,1>;

v[1] := Vector(%id = 12966956)

v[2] := Vector(%id = 12965876)

>    x[0]:=c[1]*'v[1]' + c[2]*'v[2]';

x[0] := c[1]*v[1]+c[2]*v[2]

>    x[k]:=c[1]*0.8^k*v[1] + c[2]*0.64^k*v[2];

x[k] := c[1]*.8^k*Vector(%id = 12966956)+c[2]*.64^k*Vector(%id = 12965876)

x[k] -> 0, kun k -> infinity

Tapa, jolla nollaa lähestyminen tapahtuu, on kiinnostava.

>    xnollat1:=[seq(<-3+j,3>,j=0..6)];

xnollat1 := [Vector(%id = 13097332), Vector(%id = 12176892), Vector(%id = 20232872), Vector(%id = 14805280), Vector(%id = 13781032), Vector(%id = 13097112), Vector(%id = 13096952)]

>    xnollat2:=seq(<-3+j,-3>,j=0..6);

xnollat2 := Vector(%id = 19832344), Vector(%id = 4328240), Vector(%id = 13096752), Vector(%id = 15308636), Vector(%id = 13485336), Vector(%id = 12302564), Vector(%id = 13101324)

>    x[0]:=xnollat1[1]:

>    n:=10: for k to n do x[k]:=A.x[k-1]: end do:

>    seq(x[k],k=0..n);

Vector(%id = 13097332), Vector(%id = 13624120), Vector(%id = 19357872), Vector(%id = 4301904), Vector(%id = 4369832), Vector(%id = 4438356), Vector(%id = 4444340), Vector(%id = 4462988), Vector(%id = 4...
Vector(%id = 13097332), Vector(%id = 13624120), Vector(%id = 19357872), Vector(%id = 4301904), Vector(%id = 4369832), Vector(%id = 4438356), Vector(%id = 4444340), Vector(%id = 4462988), Vector(%id = 4...
Vector(%id = 13097332), Vector(%id = 13624120), Vector(%id = 19357872), Vector(%id = 4301904), Vector(%id = 4369832), Vector(%id = 4438356), Vector(%id = 4444340), Vector(%id = 4462988), Vector(%id = 4...

>    map(convert,[%],list);

[[-3, 3], [-2.40000000000000034, 1.91999999999999993], [-1.92000000000000037, 1.22879999999999989], [-1.53600000000000048, .786431999999999910], [-1.22880000000000056, .503316479999999955], [-.98304000...
[[-3, 3], [-2.40000000000000034, 1.91999999999999993], [-1.92000000000000037, 1.22879999999999989], [-1.53600000000000048, .786431999999999910], [-1.22880000000000056, .503316479999999955], [-.98304000...
[[-3, 3], [-2.40000000000000034, 1.91999999999999993], [-1.92000000000000037, 1.22879999999999989], [-1.53600000000000048, .786431999999999910], [-1.22880000000000056, .503316479999999955], [-.98304000...
[[-3, 3], [-2.40000000000000034, 1.91999999999999993], [-1.92000000000000037, 1.22879999999999989], [-1.53600000000000048, .786431999999999910], [-1.22880000000000056, .503316479999999955], [-.98304000...
[[-3, 3], [-2.40000000000000034, 1.91999999999999993], [-1.92000000000000037, 1.22879999999999989], [-1.53600000000000048, .786431999999999910], [-1.22880000000000056, .503316479999999955], [-.98304000...

>    plotdots(%);

[Maple Plot]

>    jono[1]:=seq([x[k]],k=0..5);

jono[1] := [Vector(%id = 13097332)], [Vector(%id = 13624120)], [Vector(%id = 19357872)], [Vector(%id = 4301904)], [Vector(%id = 4369832)], [Vector(%id = 4438356)]
jono[1] := [Vector(%id = 13097332)], [Vector(%id = 13624120)], [Vector(%id = 19357872)], [Vector(%id = 4301904)], [Vector(%id = 4369832)], [Vector(%id = 4438356)]

Näiden kokeilujen jälkeen olemme valmiit kirjoittamaan pienen proseduurin.

>    traje:=proc(A,x0,n)
local x,k;
  x[0]:=x0:
  for k to n do x[k]:=A.x[k-1]: end do:
plotdots(map(convert,[seq(x[k],k=0..n)],list));
end:

>    traje(A,<-3,3>,10);

[Maple Plot]

>   

>    xnollat1;

[Vector(%id = 13097332), Vector(%id = 12176892), Vector(%id = 20232872), Vector(%id = 14805280), Vector(%id = 13781032), Vector(%id = 13097112), Vector(%id = 13096952)]

>    yla:=seq(traje(A,xnollat1[k],10),k=1..nops(xnollat1)):

>    ala:=seq(traje(A,xnollat2[k],10),k=1..nops(xnollat1)):

>    display(yla,ala,title="Ominaisarvot < 1");

[Maple Plot]

Origo on attraktiivinen kiintopiste.

Esim. 2, attraktiivinen kiintopiste

Otetaan sama uudelleen, ilman ohjelmankehitystä.

>    read("c:\\usr\\heikki\\v3\\v302.mpl"):

>    print(traje);

proc (A, x0, n) local x, k; x[0] := x0; for k to n do x[k] := A.x[k-1] end do; plotdots(map(convert,[seq(x[k],k = 0 .. n)],list)) end proc
proc (A, x0, n) local x, k; x[0] := x0; for k to n do x[k] := A.x[k-1] end do; plotdots(map(convert,[seq(x[k],k = 0 .. n)],list)) end proc
proc (A, x0, n) local x, k; x[0] := x0; for k to n do x[k] := A.x[k-1] end do; plotdots(map(convert,[seq(x[k],k = 0 .. n)],list)) end proc
proc (A, x0, n) local x, k; x[0] := x0; for k to n do x[k] := A.x[k-1] end do; plotdots(map(convert,[seq(x[k],k = 0 .. n)],list)) end proc

>    A:=<<.8,0>|<0,.64>>;

A := Matrix(%id = 12839120)

>    lambda[1]:=0.8: lambda[2]:=0.64:

>    v[1]:=<1,0>: v[2]:=<0,1>:

>    x[0]:=c[1]*'v[1]' + c[2]*'v[2]';

x[0] := c[1]*v[1]+c[2]*v[2]

>    x[k]:=c[1]*0.8^k*v[1] + c[2]*0.64^k*v[2];

x[11] := c[1]*Vector(%id = 5205500)+c[2]*Vector(%id = 5186920)

Pisteet

>    seq(X[k],k=0..10),"o o o";

X[0], X[1], X[2], X[3], X[4], X[5], X[6], X[7], X[8], X[9], X[10],

muodostavat dynaamisen systeemin trajektorin .

>    traje(A,<-3,3>,10);

[Maple Plot]

>    yla:=display(traje(A,<-3,3>,10),traje(A,<0,3>,10),traje(A,<1,3>,10),traje(A,<3,3>,10)):

>    ala:=display(traje(A,<-3,-3>,10),traje(A,<0,-3>,10),traje(A,<1,-3>,10),traje(A,<3,-3>,10)):

>    display(yla,ala,title="Ominaisarvot < 1, attraktori");

>   

[Maple Plot]

>    display(traje(A,<0,3>,10),traje(A,<0,-3>,10),traje(A,<3,0>,10),traje(A,<-3,0>,10),title="Ominaisvektorit");

[Maple Plot]

>    display(%,view=[-1..1,-1..1]); # Näin voi zoomata.

[Maple Plot]

Siis O on attraktiivinen kiintopiste , kun molemmat ominaisarvot positiivisia ja < 1.

Esim 3, repulsiivinen kiintopiste

>    read("c:\\usr\\heikki\\v3\\v302.mpl"):

>    A:=<<1.44,0>|<0,1.2>>;

A := Matrix(%id = 18676500)

>    display(traje(A,<.1,0>,10),traje(A,<-.1,0>,10),traje(A,<0,.1>,10),
traje(A,<0,-.1>,10),title="Ominaisvektorit, virtaus poispäin O:sta");

[Maple Plot]

>    display(%,traje(A,<.1,.1>,10),traje(A,<-.1,.1>,10),traje(A,<-.1,-.1>,10),traje(A,<.1,-.1>,10));

[Maple Plot]

>    #xnollat1:=[seq(.1*<-3+j,3>,j=0..6)]: #xnollat2:=seq(.1*<-3+j,-3>,j=0..6):xnollat3:=[.1*<-3,0>,.3*<3,0>]:

>   

>    #yla:=seq(traje(A,xnollat1[k],10),k=1..nops(xnollat1)):ala:=seq(traje(A,xnollat2[k],10),k=1..nops(xnollat1)):

>    #xaks:=seq(traje(A,xnollat3[k],10),k=1..nops(xnollat3)):

>    #display(yla,ala,traje(A,<.1,0>,10),traje(A,<-.1,0>,10),title="Ominaisarvot > 1");

Esim. 4, satulapiste

>    read("c:\\usr\\heikki\\v3\\v302.mpl"):

>    A:=<<2,0>|<0,0.5>>;

A := Matrix(%id = 13359528)

>    display(traje(A,<0.05,2>,10),traje(A,<-.05,2>,10),traje(A,<.1,2>,10),
traje(A,<-.1,2>,10),traje(A,<0.05,-2>,10),traje(A,<-.05,-2>,10),traje(A,<.1,-2>,10),
traje(A,<-.1,-2>,10),traje(A,<0.05,0>,10),traje(A,<-.05,0>,10),traje(A,<0,.2>,10),
traje(A,<0,-2>,10),view=[-5..5,-2..2],title="Ominaisarvot < 1, satulapiste");

[Maple Plot]

>    display(trajedot(A,<0.05,0>,10),trajedot(A,<-.05,0>,10),trajedot(A,<0,2>,10),
trajedot(A,<0,-2>,10),view=[-5..5,-2..2],title="Ominaisvektorit",axes=frame);

[Maple Plot]

x-akselilla lähdettiin läheltä O:a, y-akselilla kaukaa.  Edellisellä virtaus ulospäin, jälkimmäisellä sisäänpäin.
Voisi olla mukava sellainen traje -versio, jossa käytetään arrow :ta. Teepä harjoituksissa!

Esim. 6, Kompleksiset ominaisvektorit

>    A:=<<.8,-.1>|<.5,1>>;

A := Matrix(%id = 13006208)

>    (oa,ov):=Eigenvectors(A);

oa, ov := Vector(%id = 13584932), Matrix(%id = 14264352)
oa, ov := Vector(%id = 13584932), Matrix(%id = 14264352)

>    lambda:=oa[1]; w:=ov[1..-1,1];

lambda := .899999999999999912+.200000000000000010*I

w := Vector(%id = 15342784)

>    u:=map(Re,w);

u := Vector(%id = 15342824)

>    v:=map(Im,w);

v := Vector(%id = 15342864)

>    S:=<u|v>; SI:=S^(-1);

S := Matrix(%id = 19073072)

SI := Matrix(%id = 4514136)

>    alpha:=Re(lambda): beta:=Im(lambda):

>    R:=<<alpha,-beta>|<beta,alpha>>;

R := Matrix(%id = 4469548)

>    A=S.R.SI;

Matrix(%id = 13006208) = Matrix(%id = 4606876)

>    r:=abs(lambda);

r := .9219544457

>    phi:=argument(lambda);

phi := .2186689459

>    spiraali:=r*<<cos(phi),-sin(phi)>|<sin(phi),cos(phi)>>;

spiraali := Matrix(%id = 4896344)

>    traje(spiraali,<0,2.5>,50);# Uudet muuttujat peruskoordinaatistossa esitettynä.

[Maple Plot]

>    traje(A,<0,2.5>,50); # Iteroidaan alkuperäisellä matriisilla, kuva saataisiin
#edellisestä koordinaatistomuunnoksella, jonka välittää matriisi P.

[Maple Plot]

3. Täpläpöllöjen selviytymisongelma

>    restart:

>    with(linalg):
with(LinearAlgebra):

with(plots):

>    read("c:\\usr\\heikki\\v3\\v302.mpl"):

>    print(traje);

proc (A, x0, n) local x, k; x[0] := x0; for k to n do x[k] := A.x[k-1] end do; plotdots(map(convert,[seq(x[k],k = 0 .. n)],list)) end proc

>    A:=<<0,.18,0>|<0,0,.71>|<.33,0,.94>>;

A := Matrix(%id = 15101540)

>    A,<"juveniles","subadults","adults">;

Matrix(%id = 15101540), Matrix(%id = 15264816)

>    x[k]=<'j[k]','s[k]','a[k]'>;

x[k] = Vector(%id = 15078180)

>    (oa,ov):=Eigenvectors(A);

oa, ov := Vector(%id = 15392840), Matrix(%id = 14031208)
oa, ov := Vector(%id = 15392840), Matrix(%id = 14031208)
oa, ov := Vector(%id = 15392840), Matrix(%id = 14031208)
oa, ov := Vector(%id = 15392840), Matrix(%id = 14031208)
oa, ov := Vector(%id = 15392840), Matrix(%id = 14031208)
oa, ov := Vector(%id = 15392840), Matrix(%id = 14031208)

>    lambda[1]:=oa[1];

lambda[1] := -.217963698823998160e-1+.205918480469973492*I

>    lambda[3]:=oa[3];

lambda[3] := .983592739764799662+0.*I

>    Im(%);

0.

>    lambda[3]:=Re(lambda[3]);

lambda[3] := .983592739764799662

>    abs(lambda[1]);

.2070688348

Kaikkien ominaisarvojen itseisarvot < 1, paha merkki!

>    x[k] = c[1]*lambda[1]^k*v[1]+c[2]*lambda[2]^k*v[2]+c[3]*lambda[3]^k*v[3]

Jos lähdetään reaalisesta vektorista x[0] ,  niin x[1] = A*x[0]

Jne., siten tulokset ovat reaalisia, vaikka ominaisvektorikaavassa laskenta kulkee kompleksilukujen kautta.

Koska ominaisarvojen itseisarvot  < 1, niin x[k]  -> 0,  kun k -> infinity

Malli siis ennustaa täpläpöllöjen tuhoa.

Harjoitustehtävänä (LV6) voidaan katsoa, minkälaista "kuoleman spiraalia" pitkin kuljetaan.

Toisena harjoitustehtävänä katsotaan tilannetta, jossa nuorten pöllöjen pesäaluesiirtymävaiheen henkiinjäämisprosentti on 50, jolloin

A[2,1]=0.3 aiemman 0.18:n sijasta.(Nuorison henkiinjäämisprosentti oman pesän piirissä oli 60.)

>    A:=<<0,.3,0>|<0,0,.71>|<.33,0,.94>>;

A := Matrix(%id = 15833912)

>    A,<"juveniles","subadults","adults">;

Matrix(%id = 15833912), Matrix(%id = 14483068)

>   

>   

Tästä voit jatkaa harjoitustehtävässä