Amaalis07.html

Tentti 31.3.2007

Alustukset

> restart:

> with(LinearAlgebra):with(linalg):alias(Inv=MatrixInverse,Diag=DiagonalMatrix,Det=Determinant,Id=IdentityMatrix,ref=GaussianElimination,rref=ReducedRowEchelonForm):

Warning, the name GramSchmidt has been rebound

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

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

Katsotaan ensin tulos ja periaate valmiista Maple-funktiota k�ytt�en.

> <A | Id(3)>,rref(<A|Id(3)>); # A:n viereen 3 x 3-yksikkömatriisi ja rref koko h�sk�lle k�ytt�en Maplen valmista rref-funktiota

Koska vasemmalle saadaan yksikk�matriisi, ts. kaikki sarakkeet ovat tukisarakkeita, on k��nteismatriisi olemassa ja se koostuu yll� olevista

kolmesta viimeisest� rivist�.

Tehd��n laskut nyt vaiheittain (Maplella k�sinlaskua simuloiden):

> AE:=<A | Id(3)>;

> AE:=Matrix(swaprow(AE,1,2)); # Vaihdetaan 1. ja 2. rivi

> AE[3,1..-1]:=AE[3,1..-1]-4*AE[1,1..-1]: AE; # rivi3:=-4*rivi1+rivi3

> AE[3,1..-1]:=AE[3,1..-1]+3*AE[2,1..-1]:AE; # rivi3:=3*rivi2+rivi3

T�st� n�hd��n, ett� A:n kaikki sarakkeet ovat tukisarakkeita ja niin ollen k��nteismatriisi on olemassa.

> AE[3,1..-1]:=AE[3,1..-1]/2:AE; # Skaalataan alin tukialkio 1:ksi.

Alhaalta yl�s:

> AE[2,1..-1]:=AE[2,1..-1]-2*AE[3,1..-1]:AE; # rivi2:=-3*rivi3+rivi2

> AE[1,1..-1]:=AE[1,1..-1]-3*AE[3,1..-1]:AE; # rivi1:=-3*rivi3+rivi1

P��stiin siten samaan tulokseen kuin yll� valmiilla rref-funktiolla ja k��nteismatriisi on siten

Huomioita: T�ss� ei todellakaan tarvitse laskea yht��n determinanttia, johtop��t�s tulee suoraan laskenta-algoritmista.

Teoriassa voisi esiinty� ratkaisu, jossa johtop��t�s olemassaolosta tehd��n determinantista. Kyll� sekin hyv�ksyt��n, vaikka

onkin hiukan "huonoa makua" osoittava.

> A:=<<5,0,0,0>|<-2,3,0,0>|<6,c,5,0>|<-1,0,4,1>>;

Koska kyseess� on yl�kolmiomatriisi, ovat ominaisarvot matriisin diagonaalialkiot. Jos et sattuisi t�t� muistamaan, saat sen suoraan kehitt�m�ll�

n (alideterminantteineen) ensimm�isen sarakkeen mukaan.

Matriisi (n x n) on diagonalisoituva, jos sill� on n LRT ominaisvektoria, mik� on yht�pit�v�� sen kanssa, ett� kunkin ominaisrvon algebrallinen ja

geometrinen kertaluku yhtyv�t. T�ss� tapauksessa n�in on jos ja vain jos ominaisarvon geometrinen kertaluku yhtyy algebralliseen = 2

Toisin sanoen, laskettava ominaisarvoa vastaavat ominaisvektorit ja j�rjestett�v� parametrin c avulla tilanne sellaiseksi, ett� saadaan

2 LRT ominaisvektoria.

> (lam,V):=Eigenvectors(A);

> eigenvectors(A);

Kumpikaan ei ota huomioon mahdollisuutta: c=6, poikkeusarvot on aina ja kaikkialla k�sitelt�v� erikseen.

> A-5*Id(4); # A - 5*I

Ominaisvektoriyht�l� on siis . Merkit��n B :ll� M5:n kanssa riviekvivalenttia matriisia, joka syntyy rivioperaatioilla.

(Olkoon B yleisnimi t�llaisille.)

Seuraava rivioperaatio riippuu siit�, onko vai ei. Selvyyden vuoksi nollataan ensin alin rivi:

Jos c ei ole 6, on systeemin matriisilla 3 tukialkiota (-2,c-6,4) ja siis yksi vapaa parametri, ja siis 1-ulotteinen ominaisavaruus.

Jos c=6, on tukialkioita vain 2 (-2 ja 1), joten vapaita parametreja on 2 ja siten ominaisavaruus on 2-ulotteinen.

Siis matriisi on diagonalisoituva, jos ja vain jos c=6.

Voidaan vied� selvyyden vuoksi loppuun saakka : c=6, lis�t��n -4*rivi2 kolmanteen riviin:

Tukialkiot sarakkeissa 2 ja 4, vapaat muuttujat sarakkeissa 1 ja 3.

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

Warning, the name GramSchmidt has been rebound

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

Warning, the name changecoords has been redefined

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

Aloitetaan laskemalla ominaisarvot ja -vektorit. Sit� laskua en en�� n�yt�, kaikkihan sen osaavat!

Antaa Maplen laskea valmiilla funktiolla:

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

Siis ominaisarvot ovat -1 ja 3 ja vastaavat ominaisvektorit ovat ov-matriisin sarakkeet:

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

Yleinen ratkaisu:

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

Alkuarvoteht�v�:

> y(0)=<1,3>;

> LinearSolve(ov,<1,3>); # Ratkaistaan vakiot C1 ja C2:

> ya:=unapply(subs(C1=1/2,C2=5/2,y(t)),t): # Maple-tekniikkaa, ei tarvitse lukea.

> ya(t); # Alkuarvoteht�v�n ratkaisu koordinaattimuodossa:

> Ya:=ya(t):

> plot([Ya[1],Ya[2],t=-2..2],view=[0..5,0..5]);

> AAratkuva:=%:

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

> YY:=evalm(Ys);

> #YYa:=subs(C1=1/2,C2=1/2,op(YY));

> #YYa;

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

> parvikuva:=plot([parvi]):

> display(parvikuva,AAratkuva,plot([[1,3]],style=point,symbol=circle,symbolsize=17),arrow([v1,v2]),view=[-5..5,-5..5]);

> kuva:=%:

> with(DEtools):

Warning, the previous bindings of the names RationalCanonicalForm and adjoint have been removed and they now have an assigned value

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

> display(kuva,%,title="Satula");

Selityksi�:

Yleinen ratkaisu:

miss�

> 'v1'=v1, 'v2'=v2;

Ominaisvektorilla (C2=0) kuljetaan O:a kohti () ja ominaisvektorilla v[2] (C1=0) O:sta poisp�in ().

Kyseess� on satulapiste, joka on ep�stabiili (jo pelk�st��n ominaisvektorik�yt�s kertoo sen).

Suuntakentt�nuolet annetuissa pisteiss� lasketaan yksinkertaisesti n�in:

> p1:=<1,0>; s1:=A.p1: s1:=s1/10;#arrow([p1],[s1/10]);

Piirrettiin kuhunkin annettuun pisteeseen p derivaatan suuntainen nuoli s, joka saadaan kertomalla , kun kerran y'=A p. Normeerataan t�ss� jakamalla 10:ll�.

Kriittiset pisteet: Ratkaistaan yht�l�ryhm�: "oikeat puolet = 0", ts.

> yhtaloryhma:={200*x-4*x*y=0,-150*y+2*x*y=0};

> solve(yhtaloryhma,{x,y});

K�sin laskien: n�hd��n v�litt�m�sti ratkaisuksi (sukupuuttoratkaisu). Jos x ei ole 0, saadaan 1. yht�l�st�: ja sitten toisesta .

Muita mahdollisuuksia ei ole.

> krp:=<75,50>; # Hengiss� s�ilymisen kriittinen piste

Linearisointi:

> f:=200*x-4*x*y;g:=-150*y+2*x*y;

> J:=Matrix(jacobian([f,g],[x,y])); # Jacobin matriisi

> A:=subs(x=75,y=50,J); # Sijoitetaan kriittinen piste

> Eigenvalues(A);

Puhtaasti imaginaariset ominaisarvot, KRP:n tyyppi on siten keskus ja siis stabiili.

Huom: Keskuksen tapauksessa linearisoidun systeemin (O:n) tyyppi ei v�ltt�m�tt� kerro alkuper�isen systeemin (KRP:n) tyyppi�.

T�ss� tapauksessa se on sama, mutta sit� ei teht�v�ss� pyydetty pohtimaan (keinoja siihen ei kurssilla ole edes opetettu).

Teht�v� on periaatteessa samanlainen kuin kahdessa edellisess� tentiss�. Nyt on alkuarvofunktio m��ritelty paloittain kahdessa osassa.

Se on tehty edelleenkin mahdollisimman humaanisti integrointivaivoja s��st�en.

�l� v�lit� ylim��r�isist� Maple-komennoista. (Vaikka ovat punaisella, niin �l� sin� ala n�hd� punaista!)

> with(plots):setoptions3d(orientation=[-30,50],axes=box):

Warning, the name changecoords has been redefined

> read("c:/usr/heikki/numsym05/maple/ns05.mpl");

Ratkaisun muoto on annettu.

Jotta alkuehto toteutuisi, on - kertoimet valittava alkuehtofunktion f sinisarjan kertoimiksi: .

Tarvitsemme siis sinisarjaa.

> bn:=(2/L)*int(f(x)*sin(n*Pi*x/L),x=0..L);

> bn := (2/10)*(Int(100*sin(n*Pi*x/10),x=0..5)+Int(O*sin(n*Pi*x/10),x=5..10));

> seq(-1+cos((n*Pi)/2),n=1..10);

Kertoimet ovat siis muotoa:

Muodostetaan yleinen n:s osasumma:

> lambda:=n-> n*Pi/10; summaf:=(x,t,N)-> sum(B[n]*sin((n*Pi*x)/10)*exp(-lambda(n)^2*t),n=1..N);

Kysytty 3:n termin summa on siis:

> plot3d(summaf(x,t,3),x=0..10,t=0..20,title="n=3");plot3d(summaf(x,t,20),x=0..10,t=0..20,title="n=20");plot3d(summaf(x,t,40),x=0..10,t=0..20,title="n=40");

Huomataan, ett� n=3 on varsin karkea alkuapproksimaatio, n=20 n�ytt�� Gibbsin ilmi�n hyvin selv�sti. Toki Gibbsin ilmi� esiintyy isommillakin,

mutta kokonaisuus n�kyy sile�mp�n� isommalla n:n arvolla. (Gibbs-piikit ovat hyvin l�hell� ep�jatkuvuuskohtia, eik� 3D-grafiikka v�lttm�tt� n�yt�

niit�.) Kun aikaa kuluu, niin ratkaisu sili�� varsin nopeasti, johtuen exp-termeist�.