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>>;

A := Matrix(%id = 17381784)

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

lambda, ov := Vector(%id = 20234552), Matrix(%id = 16315852)

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

x := [Vector(%id = 20022384), Vector(%id = 20022944)]

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

Yh := C1*exp(-1/4*t)*Vector(%id = 20022384)+C2*exp(-1/20*t)*Vector(%id = 20022944)

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

>    yp:=<v1,v2>;

yp := Vector(%id = 21044072)

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

b := Vector(%id = 21209584)

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

Vector(%id = 21185100) = A.yp+Vector(%id = 21209584)

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

Yp := Vector(%id = 21289740)

>   

>    Y:=Yh+Yp;

Y := C1*exp(-1/4*t)*Vector(%id = 20022384)+C2*exp(-1/20*t)*Vector(%id = 20022944)+Vector(%id = 21289740)

>    diff(Y,t);

-1/4*C1*exp(-1/4*t)*Vector(%id = 20022384)-1/20*C2*exp(-1/20*t)*Vector(%id = 20022944)

>    A.Y+b;

Matrix(%id = 17381784).(C1*exp(-1/4*t)*Vector(%id = 20022384)+C2*exp(-1/20*t)*Vector(%id = 20022944)+Vector(%id = 21289740))+Vector(%id = 21209584)

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,%);

AE := C1*exp(0)*Vector(%id = 20022384)+C2*exp(0)*Vector(%id = 20022944)+Vector(%id = 21289740) = Vector(%id = 21290100)

AE := C1*Vector(%id = 20022384)+C2*Vector(%id = 20022944)+Vector(%id = 21289740) = Vector(%id = 21290100)

>    evalm(lhs(AE));

vector([C1+3/2*C2+42, -2*C1+C2+36])

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

Ceet := Vector(%id = 22332032)

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

Y := 29/8*exp(-1/4*t)*Vector(%id = 20022384)-55/4*exp(-1/20*t)*Vector(%id = 20022944)+Vector(%id = 21289740)

>    y12:=evalm(Y);

y12 := vector([29/8*exp(-1/4*t)-165/8*exp(-1/20*t)+42, -29/4*exp(-1/4*t)-55/4*exp(-1/20*t)+36])

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

y1 := 29/8*exp(-1/4*t)-165/8*exp(-1/20*t)+42

y2 := -29/4*exp(-1/4*t)-55/4*exp(-1/20*t)+36

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

[Maple Plot]

>    <diff(y1,t),diff(y2,t)>;

Vector(%id = 3081680)

>    A.<y1,y2>+b;

Vector(%id = 19705252)

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(t) = cos(t)  , y(t) = sin(t)  .

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

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

[Maple Plot]

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

[Maple Plot]

>   

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>>;

A := Matrix(%id = 17182076)

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

lambda, ov := Vector(%id = 17896372), Matrix(%id = 21551204)

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

x1 := Vector(%id = 18182272)

x2 := Vector(%id = 17897052)

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

Y := C1*exp(-4*t)*Vector(%id = 18182272)+C2*exp(-2*t)*Vector(%id = 17897052)

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

[Maple Plot]

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

>    s=exp(-2*t);

s = exp(-2*t)

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

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

Ys := C1*s^2*Vector(%id = 18182272)+C2*s*Vector(%id = 17897052)

>    YY:=evalm(Ys);

YY := vector([-C1*s^2+C2*s, C1*s^2+C2*s])

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

plot([-arvo1*s^2+arvo2*s, arvo1*s^2+arvo2*s, s = a .. b])

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

[Maple Plot]

>   

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

>    display(kuva1,kuva2);

[Maple Plot]

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

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

parvi := [2*s^2-2*s, -2*s^2-2*s, s = 0 .. 1], [2*s^2-s, -2*s^2-s, s = 0 .. 1], [2*s^2, -2*s^2, s = 0 .. 1], [2*s^2+s, -2*s^2+s, s = 0 .. 1], [2*s^2+2*s, -2*s^2+2*s, s = 0 .. 1], [s^2-2*s, -s^2-2*s, s =...
parvi := [2*s^2-2*s, -2*s^2-2*s, s = 0 .. 1], [2*s^2-s, -2*s^2-s, s = 0 .. 1], [2*s^2, -2*s^2, s = 0 .. 1], [2*s^2+s, -2*s^2+s, s = 0 .. 1], [2*s^2+2*s, -2*s^2+2*s, s = 0 .. 1], [s^2-2*s, -s^2-2*s, s =...
parvi := [2*s^2-2*s, -2*s^2-2*s, s = 0 .. 1], [2*s^2-s, -2*s^2-s, s = 0 .. 1], [2*s^2, -2*s^2, s = 0 .. 1], [2*s^2+s, -2*s^2+s, s = 0 .. 1], [2*s^2+2*s, -2*s^2+2*s, s = 0 .. 1], [s^2-2*s, -s^2-2*s, s =...
parvi := [2*s^2-2*s, -2*s^2-2*s, s = 0 .. 1], [2*s^2-s, -2*s^2-s, s = 0 .. 1], [2*s^2, -2*s^2, s = 0 .. 1], [2*s^2+s, -2*s^2+s, s = 0 .. 1], [2*s^2+2*s, -2*s^2+2*s, s = 0 .. 1], [s^2-2*s, -s^2-2*s, s =...
parvi := [2*s^2-2*s, -2*s^2-2*s, s = 0 .. 1], [2*s^2-s, -2*s^2-s, s = 0 .. 1], [2*s^2, -2*s^2, s = 0 .. 1], [2*s^2+s, -2*s^2+s, s = 0 .. 1], [2*s^2+2*s, -2*s^2+2*s, s = 0 .. 1], [s^2-2*s, -s^2-2*s, s =...
parvi := [2*s^2-2*s, -2*s^2-2*s, s = 0 .. 1], [2*s^2-s, -2*s^2-s, s = 0 .. 1], [2*s^2, -2*s^2, s = 0 .. 1], [2*s^2+s, -2*s^2+s, s = 0 .. 1], [2*s^2+2*s, -2*s^2+2*s, s = 0 .. 1], [s^2-2*s, -s^2-2*s, s =...
parvi := [2*s^2-2*s, -2*s^2-2*s, s = 0 .. 1], [2*s^2-s, -2*s^2-s, s = 0 .. 1], [2*s^2, -2*s^2, s = 0 .. 1], [2*s^2+s, -2*s^2+s, s = 0 .. 1], [2*s^2+2*s, -2*s^2+2*s, s = 0 .. 1], [s^2-2*s, -s^2-2*s, s =...

>    plot([parvi]);

[Maple Plot]

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

[Maple Plot]

>    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);

[Maple Plot]

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

[Maple Plot]

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>>;

A := Matrix(%id = 16910140)

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

lambda, ov := Vector(%id = 3182904), Matrix(%id = 16709764)

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

v1 := Vector(%id = 17508248)

v2 := Vector(%id = 17552412)

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

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

y := proc (t) options operator, arrow; C1*exp(lambda[1]*t)*v1+C2*exp(lambda[2]*t)*v2 end proc

>    Y:=y(t);

Y := C1*exp(t)*Vector(%id = 17508248)+C2*exp(2*t)*Vector(%id = 17552412)

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

[Maple Plot]

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

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

Ys := C1*s*Vector(%id = 17508248)+C2*s^2*Vector(%id = 17552412)

>    YY:=evalm(Ys);

YY := vector([C1*s+C2*s^2, C2*s^2])

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

>    plot([parvi]);

[Maple Plot]

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

[Maple Plot]

>    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);

[Maple Plot]

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

[Maple Plot]

Esim. 2.3, Nielu, ominaisarvot < 0

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

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

A := Matrix(%id = 2646876)

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

lambda, ov := Vector(%id = 2491888), Matrix(%id = 3180340)

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

v1 := Vector(%id = 19087640)

v2 := Vector(%id = 16352512)

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

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

y := proc (t) options operator, arrow; C1*exp(lambda[1]*t)*v1+C2*exp(lambda[2]*t)*v2 end proc

>    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.

Ys := C1/s^3*Vector(%id = 19087640)+C2/s*Vector(%id = 16352512)

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

Lasketaan koordinaattimuodossa:

>    YY:=evalm(Ys);

YY := vector([C1/s^3-C2/s, C1/s^3+C2/s])

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

>    plot([parvi]);

[Maple Plot]

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

[Maple Plot]

>    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");

[Maple Plot]

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

Esim 2.4 Satula, reaaliset, erimerkk.

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

A := Matrix(%id = 3120688)

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

lambda, ov := Vector(%id = 3086900), Matrix(%id = 19189308)

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

v1 := Vector(%id = 3021452)

v2 := Vector(%id = 2842128)

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

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

Y := C1*exp(3*t)*Vector(%id = 3021452)+C2*exp(-2*t)*Vector(%id = 2842128)

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

C1*exp(3*ln(s))*Vector(%id = 3021452)+C2*exp(-2*ln(s))*Vector(%id = 2842128)

Ys := C1*s^3*Vector(%id = 3021452)+C2/s^2*Vector(%id = 2842128)

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

>    YY:=evalm(Ys);

YY := vector([C1*s^3+C2/s^2, -C1*s^3+4*C2/s^2])

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

>    plot([parvi]);

[Maple Plot]

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

[Maple Plot]

>    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");

[Maple Plot]

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

ominaisvektorin pisteistä).

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

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

A := Matrix(%id = 16942800)

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

lambda, ov := Vector(%id = 17309028), Matrix(%id = 16481584)

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

v1 := Vector(%id = 17519860)

v2 := Vector(%id = 17619196)

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

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

Y := C1*exp((1+8*I)*t)*Vector(%id = 17519860)+C2*exp((1-8*I)*t)*Vector(%id = 17619196)

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

alpha := 1

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

beta := 8

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

u := Vector(%id = 17685188)

v := Vector(%id = 17680892)

>    uv:=<u|v>;

uv := Matrix(%id = 2703092)

>    C:=<<cos(beta*t),-sin(beta*t)>|<sin(beta*t),cos(beta*t)>>;

C := Matrix(%id = 2929688)

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

Y := exp(t)*Vector(%id = 18275120)

>    YY:=evalm(Y);

YY := vector([exp(t)*(cos(8*t)*C1+sin(8*t)*C2), exp(t)*((cos(8*t)+sin(8*t))*C1+(sin(8*t)-cos(8*t))*C2)])

>   

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

>    plot([parvi]);

[Maple Plot]

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

[Maple Plot]

>    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");

[Maple Plot]

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

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

A := Matrix(%id = 2780164)

>    Eigenvectors(A);

Vector(%id = 16319216), Matrix(%id = 17216456)

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

expA := proc (t) options operator, arrow; exp(3*t)*`<|>`(`<,>`(1+t,-t),`<,>`(t,1-t)) end proc

>    expA(t);

exp(3*t)*Matrix(%id = 17998988)

>    exponential(A,t); # Tarkistus.

matrix([[exp(3*t)+t*exp(3*t), t*exp(3*t)], [-t*exp(3*t), exp(3*t)-t*exp(3*t)]])

>    Y:=exp(3*t)*((expA(t)/exp(3*t)).<C1,C2>);

Y := exp(3*t)*Vector(%id = 21355320)

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

>    YY:=evalm(Y);

YY := vector([exp(3*t)*((1+t)*C1+t*C2), exp(3*t)*(-t*C1+(1-t)*C2)])

>   

>    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");

[Maple Plot]

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));

[Maple Plot]

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

[Maple Plot]

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>>;

A := Matrix(%id = 2976816)

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

lambda, ov := Vector(%id = 16497644), Matrix(%id = 16488400)

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

x1 := Vector(%id = 16568092)

x2 := Vector(%id = 16537988)

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

Y := C1*exp(-4*t)*Vector(%id = 16568092)+C2*exp(-2*t)*Vector(%id = 16537988)

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

Y := [-C1*exp(-4*t)+C2*exp(-2*t), C1*exp(-4*t)+C2*exp(-2*t)]

>    A;

Matrix(%id = 2976816)

>    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);

[Maple Plot]

>    #?DEplot

>   

>