L12maple.mws

Lineaariset differentiaaliyhtälösysteemit. Ke 30.10-ti 5.11.02

Alustukset

 > with(LinearAlgebra):with(linalg):

```Warning, the assigned name GramSchmidt now has a global binding
```

```Warning, the name adjoint has been redefined
```

```Warning, the previous binding of the name GramSchmidt has been removed and it now has an assigned value
```

 > with(plots): with(DEtools):

```Warning, the name adjoint has been redefined
```

Seosprobleema

 > restart:

 > with(LinearAlgebra):with(linalg):

```Warning, the previous binding of the name GramSchmidt has been removed and it now has an assigned value
```

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

 > A:=<<-1/10,1/10>|<3/40,-1/5>>;

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

 > x:=[ov[1..-1,1],ov[1..-1,2]];

 > Yh:=C1*exp(lambda[1]*t)*x[1]+C2*exp(lambda[2]*t)*x[2];

HY-yrite: Yksinkertaisinta on yrittää vakiovektoria y=v;

 > yp:=;

 > b:=<3/2,3>;

 > map(diff,yp,t)='A.yp'+b;

 > Yp:=LinearSolve(A,-b);

 >

 > Y:=Yh+Yp;

 > diff(Y,t);

 > A.Y+b;

Tarkistaminen vektorimuodossa ei käykään helposti, Maplen tietorakenteet tekevät tenän.

No, jatketaan ja jätetään tarkistus loppuun.

 > AE:=subs(t=0,Y)=<25,15>;AE:=map(eval,%);

 > evalm(lhs(AE));

 > Ceet:=LinearSolve(<<1,-2>|<3/2,1>>,<25-42,15-36>);

 > Y:=subs(C1=Ceet[1],C2=Ceet[2],Y);

 > y12:=evalm(Y);

 > y1:=y12[1];y2:=y12[2];

 > plot([y1,y2],t=0..100);

 > ;

 > A.+b;

Kävihän se tarkistus toki näin!

2x2-systeemejä, faasitasoja

Katsotaan ensin, mikä on "aikatason" ja faasitason suhde ja mikä on vastaava Maple-syntaksi.

Merkitään havainnollisuuden vuoksi faasitason pisteitä (x(t) , y(t)).

Olk.  ,  .

 > x:=cos(t): y:=sin(t):

 > plot([x,y],t=0..2*Pi,title="Aikakuva");

 > plot([x,y,t=0..2*Pi],title="Faasikuva");

 >

Esim. 1 (KRE s. 163)

 > restart: with(LinearAlgebra): with(linalg): with(plots):

```Warning, the previous binding of the name GramSchmidt has been removed and it now has an assigned value
```

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

```Warning, the name changecoords has been redefined
```

 > A:=<<-3,1>|<1,-3>>;

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

 > x1:=ov[1..-1,1]; x2:=ov[1..-1,2];

 > Y:=C1*exp(lambda[1]*t)*x1+C2*exp(lambda[2]*t)*x2;

 > C1:=1: C2:=.5: plot([Y[1],Y[2]],t=0..2,title="Ratkaisufunktiot y1(t)ja y2(t) aikakuvassa",color=[red,blue]);

Faasikuvan piirtämistä ajatellen kannattaa katsoa taas ominaisvektorikantaa, Sekä käsin, että Maplella on helpompaa, jos otetaan käyräparametriksi

 > s=exp(-2*t);

 > C1:='C1': C2:='C2':

 > Ys:=C1*s^2*x1+C2*s*x2;

 > YY:=evalm(Ys);

 > C1:=arvo1: C2:=arvo2: 'plot'([YY[1],YY[2],s=a..b]); # Huomaa param. piirron syntaksi!

 > C1:=-2: C2:=1: plot([YY[1],YY[2],s=0..1]);kuva1:=%:

 >

 > C1:=2: C2:=-1: kuva2:=plot([YY[1],YY[2],s=0..1]):

 > display(kuva1,kuva2);

 > C1:='C1': C2:='C2':

 > parvi:=seq(seq([YY[1],YY[2],s=0..1],C2=-2..2),C1=-2..2);

 > plot([parvi]);

 > display(%,arrow([x1,x2]));

 > kuva:=%:

 > with(DEtools):

```Warning, the name adjoint has been redefined
```

 > DEplot([D(y1)(t)=-3*y1(t)+y2(t),D(y2)(t)=y1(t)-3*y2(t)], [y1(t),y2(t)],t=-5..2,y1=-4..4,y2=-4..4,color=grey);

 > display(kuva,%,title="Trajektoreita ja suuntakenttä");

2x2-faasitasotyypit

KRE:ssä on perusesimerkkejä ss. 164 - 168.

Käydään tässä EN:n sivujen mukaan ja käytetään EN-nimityksiä

[EN]: Eirola-Nevanlinna: Diff. yht. sys. L3-luentoja

Esim. 2.2 Lähde, ominaisarvot > 0

 > restart: with(LinearAlgebra): with(linalg): with(plots):

```Warning, the previous binding of the name GramSchmidt has been removed and it now has an assigned value
```

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

```Warning, the name changecoords has been redefined
```

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

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

 > v1:=ov[1..-1,1]; v2:=ov[1..-1,2];

 > C:='C1': C:='C2':

 > y:=t->C1*exp(lambda[1]*t)*v1+C2*exp(lambda[2]*t)*v2;

 > Y:=y(t);

 > C1:=1: C2:=.5: plot([Y[1],Y[2]],t=0..2,title="Ratkaisufunktiot y1(t)ja y2(t) aikakuvassa",color=[red,blue]);

 > C1:='C1': C2:='C2':

 > Ys:=map(eval,subs(t=log(s),y(t)));

 > YY:=evalm(Ys);

 > parvi:=seq(seq([YY[1],YY[2],s=0..4],C2=-2..2),C1=-2..2):

 > plot([parvi]);

 > display(%,arrow([v1,v2]),view=[-5..5,-5..5]);

 > kuva:=%:

 > with(DEtools):

```Warning, the name adjoint has been redefined
```

 > DEplot([D(y1)(t)=y1(t)+y2(t),D(y2)(t)=2*y2(t)], [y1(t),y2(t)],t=-5..2,y1=-5..5,y2=-5..5,color=grey);

 > display(kuva,%,title="Lähde,lambda[1] > 0, lambda[2] > 0");

Esim. 2.3, Nielu, ominaisarvot < 0

Periaatteessa samanlainen kuin lähde edellä, mutta nyt ominaisarvot negatiiviset. Trajektirit virtaavat origoon.

 > A:=<<-2,-1>|<-1,-2>>;

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

 > v1:=ov[1..-1,1]; v2:=ov[1..-1,2];

 > C1:='C1': C2:='C2':

 > y:=t->C1*exp(lambda[1]*t)*v1+C2*exp(lambda[2]*t)*v2;

 > subs(t=log(s),y(t)): Ys:=eval(%); # Siirrytään logaritmiseen parametriin. Tämä ei ole välttämätöntä, mutta helpottaa erit. käsin piirtoa, myös Maple-piirrossa skaala pysyy paremmin hanskassa.

 > # 0 < s < 1, kun t > 0.

 > YY:=evalm(Ys);

 > parvi:=seq(seq([YY[1],YY[2],s=0..4],C2=-2..2),C1=-2..2):

 > plot([parvi]);

 > display(%,arrow([v1,v2]),view=[-5..5,-5..5]);

 > kuva:=%:

 > with(DEtools):

 > display(kuva,DEplot([D(y1)(t)=-2*y1(t)-y2(t),D(y2)(t)=-y1(t)-2*y2(t)], [y1(t),y2(t)],t=-5..2,y1=-5..5,y2=-5..5,color=grey),title="Nielu");

Nielu on vahvasti stabiili. Kaikki trajektorit lähestyvät O:a, kun t ->  .

Esim 2.4 Satula, reaaliset, erimerkk.

 > A:=<<2,-4>|<-1,-1>>;

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

 > v1:=ov[1..-1,1]; v2:=ov[1..-1,2];

 > C1:='C1': C2:='C2':

 > Y:=C1*exp(lambda[1]*t)*v1+C2*exp(lambda[2]*t)*v2;

 > subs(t=log(s),Y);Ys:=eval(%);

 > # 0 < s < 1, kun t > 0

 > YY:=evalm(Ys);

 > parvi:=seq(seq([YY[1],YY[2],s=0..4],C2=-2..2),C1=-2..2):

 > plot([parvi]);

 > display(%,arrow([v1,v2]),view=[-5..5,-5..5]);

 > kuva:=%:

 > with(DEtools):

 > display(kuva,DEplot([D(y1)(t)=2*y1(t)-y2(t),D(y2)(t)=-4*y1(t)-y2(t)], [y1(t),y2(t)],t=-5..2,y1=-5..5,y2=-5..5,color=grey),title="Satula");

Epästabiili.  Osa trajektoreista --> , lähdettiinpä miten läheltä O:a tahansa (erit. posit. ominaisarvoa vastaavan

ominaisvektorin pisteistä).

Esim. 2.5, epästabiili fokus , kompleksiset, Re( ) > 0

 > A:=<<9,16>|<-8,-7>>;

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

 > v1:=ov[1..-1,1]; v2:=ov[1..-1,2];

 > C1:='C1': C2:='C2':

 > Y:=C1*exp(lambda[1]*t)*v1+C2*exp(lambda[2]*t)*v2;

 > alpha:=Re(lambda[1]);

 > beta:=Im(lambda[1]);

 > u:=map(Re,v1); v:=map(Im,v1);

 > uv:=;

 > C:=<|>;

 > Y:=exp(alpha*t)*(uv.C.); # Sulut laitettava näin, Maple ei täysin hallitse kertolaskujen liitännäisyyttä.

 > YY:=evalm(Y);

 >

 > parvi:=seq(seq([YY[1],YY[2],t=0..2],C2=-2..2),C1=-2..2):

 > plot([parvi]);

 > display(%,arrow([u,v]),view=[-5..5,-5..5]);

 > kuva:=%:

 > with(DEtools):

 > display(kuva,DEplot([D(y1)(t)=9*y1(t)-8*y2(t),D(y2)(t)=16*y1(t)-7*y2(t)], [y1(t),y2(t)],t=-5..2,y1=-5..5,y2=-5..5,color=grey),title="Epästabiili fokus");

Esim. 3.0 Degeneroitununut lähde (KRE s. 168)

 > A:=<<4,-1>|<1,2>>;

 > Eigenvectors(A);

 > expA:=t->exp(3*t)*<<1+t,-t>|>;

 > expA(t);

 > exponential(A,t); # Tarkistus.

 > Y:=exp(3*t)*((expA(t)/exp(3*t)).);

Jälleen piti temppuilla, koska kertolaskun liitännäisyys ei ole ihan hallinnassa.

 > YY:=evalm(Y);

 >

 > parvi:=seq(seq([YY[1],YY[2],t=-2..2],C2=-2..2),C1=-2..2):

 >

 > display(plot([parvi]),arrow(<-1,1>),view=[-5..5,-5..5],title="Degeneroitunut lähde");

Vaihtoehtoisia tapoja piirtoon (Maple-teknikkaa)

Olemme harrastaneet grafiikka-arvoisia funktioita, kuten traje edellä diffrenssiyhtälöissä.

Tässäkin voisimme tehdä vastaavasti:

 > dystraje:=(C1,C2)->plot([-C1*s^2+C2*s,C1*s^2+C2*s,s=0..1]):

 > with(plots):

 > display(dystraje(1,1));

 > display(seq(seq(dystraje(C1,C2),C2=-2..2),C1=-1..1));

Kannattaisi kirjoittaa sellainen, joka laskee ja piirtää trajektorin annetusta alkupisteestä lähtien.

Esim. 1

Tässä laskemme analyyttisen ratkaisun, mutta teemme kaiken piirtämisen DEplot:lla. Se on karkeampi tapa, mutta

hyvä demota sitäkin. Kuvaan voitaisiin lisätä dystraje:lla laskettuja trajektoreita.

 > restart: with(LinearAlgebra):

 > A:=<<-3,1>|<1,-3>>;

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

 > x1:=ov[1..-1,1]; x2:=ov[1..-1,2];

 > Y:=C1*exp(lambda[1]*t)*x1+C2*exp(lambda[2]*t)*x2;

 > Y:=convert(evalm(%),list);

 > A;

 > restart:

 > with(DEtools):

 > alkuarvolista:=[[y1(0)=2,y2(0)=4],[y1(0)=1,y2(0)=1],[y1(0)=-1,y2(0)=1],[y1(0)=1,y2(0)=-1],[y1(0)=-8,y2(0)=4],[y1(0)=-4.04,y2(0)= 7.49]]:

 > DEplot([D(y1)(t)=-3*y1(t)+y2(t),D(y2)(t)=y1(t)-3*y2(t)], [y1(t),y2(t)],t=-5..2,y1=-10..10,y2=-10..10,alkuarvolista,stepsize=.05);

 > #?DEplot

 >

 >