Alustukset

Tässäpä sattui hankala paikka. Rita  oli miltei epätoivon partaalla, kun GramSchmid t ei toiminut, vaikka syöte

aivan varmasti oli oikea. Onnistuin nyt itsekin kaatumaan samaan seikkaan, vaikka aiemmin en ollut (ja ymmärsin omakohtaisesti

Ritan epätoivon).

Syy selvisi tarkkailemalla sinisiä latausvaroituksia;

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

?GramSchmidt

paljastaa, että myös linalg:ssa on aivan samanniminen funktio. Nyt leikkisästi riippuu kirjastojen latausjärjesyksestä,

kumpi on voimassa.

Lupaan kirjoittaa tästä palautetta MUG-listalle ensi tilassa (Maple Users' Group).

Kohtuullisen robusti ratkaisu on alla oleva alias-määrittely, jossa viitataan koko nimeen kirjasto mukaanlukien.

Huomaa! Aina voidaan käyttää pitkää nimeä, siis LinearAlgebra[GramSchmidt],  jolloin kirjastoa ei edes tarvitse ladata. (Voidaan ajatella, että LinearAlgebra  on ison perheen

sukunimi ja esim. GramSchmid t etunimi.)

>

restart:

>	alias(GS=LinearAlgebra[GramSchmidt]):

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

>	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

>	# ?GramSchmidt

>

>	alias(ref=GaussianElimination):

>	alias(TR=Transpose):

>	sarakkeet:=A->[seq(A[1..-1,j],j=1..LinearAlgebra[Dimension](A)[2])]; # Tämä voi olla kätevä.

>	randV:=n->RandomVector(n,generator=rand(-1000..1000)/1000.);

>	randM:=(m,n)->RandomMatrix(m,n,generator=rand(-1000..1000)/1000.);

1.

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

>

TR(A).A;

Tämä on käyttäjän kannalta tehokkain tapa. Siinähän muodostetaan juuri kaikki sisätulot , missä  on

A:n sarakevektori.

Koska saatiin diagonaalimatriisi, on ortogonaalisuus osoitettu.

Nähdään myös, että jokaisen vektorin normi = 10, joten normaaraus saadaan yksinkertaisesti:

>

U:=A/10;

>

TR(U).U;

Tietenkään ortogonaalisen systeemin vektoreilla ei yleensä ole sama normi, joten normeeraus on tehtävä vektoreittain.

>	GS(sarakkeet(A)); # Käytetään nyt riemumielin GS:ää, vaikka ei tässä olekaan tarpeen.

>	GS(sarakkeet(A),normalized);

b)

>	TR(U).U,U.TR(U);

Ainakin ne ovat eri kokoiset neliömatriisit. Kokeillaan nyt vielä tällaista:

>	(U.TR(U))^2;

Ihmeellistä,  on sama kuin sen neliö. (Idempotentti.)

c)

Kun lasketaan , niin vektorimuodossa kirjoitettuna saadaan :

         + ... +  ,

  missä =< ,y>, ja  on U:n k:s sarakevektori,

Mieti tämä läpi kynän ja paperin kanssa.

U:n sarakkeet ovat A:n sarakkeet ortonormeerattuina, joten ne virittävät saman aliavaruuden, joka on (määritelmän mukaan)

A:n sarakeavaruus.

Toinen tapa olisi käyttää ref:iä:

>	y:=randV(8);

>	p:=U.TR(U).y:

>	ref(<A|p>);

Luokkaa 10^-16 oleva luku on ilman muuta nolla. Liitännäismatriisin rangi pysyy samana kuin A:n, joten :

Tässä siis ref ei ollut tarpeen, koska ortogonaalisuus antoi tiedon suoraan.

p on y:n ortog. projektio aliavaruudella col(A), joten se on yleisen projektiolauseen mukaan sellainen, että

          (y - p) _|_ col(A)

Lasketaan nyt tässä tapauksessa vielä:

>

z:=y-p;

>

TR(z).A;

Vektorin z sisätulo jokaisen sarakkeen kanssa = 0, joten z on todellakin kohtisuorassa koko aliavaruutta col(A) vastaan

(vrt. teht. 1 AV, "lapsellinen lause").

d) Tämä esitys on avainasemassa moneen asiaan. Me osaamme ainakinnähdä p:n PNS-approksimaationa y:lle.

2.

>	data:=[[4,1.58],[6,2.08],[8,2.5],[10,2.8],[12,3.1],[14,3.4],[16,3.8],[18,4.32]];

>	xd:=map(p->p[1],data);yd:=map(p->p[2],data);

>	van:=VandermondeMatrix(xd);

>	A:=van[1..-1,2..4];

>	ATA:=TR(A).A;

>	ATy:=TR(A).<yd>;

>	c:=LinearSolve(ATA,ATy);  # Normaaliyhtälöt

>	p:=c[1]*x+c[2]*x^2 + c[3]*x^3;

>	display(plot(data,style=point,symbol=cross,symbolsize=15),plot(p,x=4..18,color=grey));

Tehdään vielä QR-hajotelmalla:

>	vvv:=seq(A[1..-1,j],j=1..3);

>	Q:=Matrix(GS([vvv],normalized));

>	R:=TR(Q).A;

>	cqr:=LinearSolve(R,TR(Q).<yd>);

>

c;

Varsin hyvällä tarkkuudella samat.

3. Pariisin lämpötilat

6.

>

restart:

Warning, the name changecoords has been redefined

>	vasen:=m*diff(y(t),t,t)+c*diff(y(t),t)+k*y(t);

>	HY:=vasen=0;

>	EHY:=vasen=r;

>	m:=2: k:=2: c:=1:

Lasketaan muutakin kuin "steady state". Olkoon tämä teoriaa kertaava ja täydentävä esimerkkiratkaisu.

1) HY:n yleinen ratkaisu:

>

HY;

>	karyht:=2*lambda^2+lambda+2=0;

>	lam:=solve(karyht,lambda);

>	alpha:=Re(lam[1]); beta:=abs(Im(lam[1])); # abs vain siksi, että Maple voi saada päähänsä vaihtaa lam-jonon alkioiden järjestystä eri ajokerroilla.

HY:n yleinen on siis:

>	yh:=exp(alpha*t)*(A*cos(beta*t)+B*sin(beta*t));

Siinä se komeilee! Vakiot A ja B määräytyvät alkuehdoista. Mutta, alkuehtoja ei sovelleta HY:öön,  kun yhtälömme  

on EHY . (Jos kyse olisi vapaasta, siis ilman ulkoista voimaa olevasta värähtelystä, niin silloin (ja vain silloin)...)

Joka tapauksessa, olivatpa alkuehdoista saatavat vakiot A ja B mitä tahansa, niin  --> 0, kun t --> .

Siksi  yh:ta sanotaan transientiksi, se edustaa nopeasti nitisyvää alkuvaihetta, joka siis negatiiviseksponenttisen exp-

funktion ansiosta pian kuolee pois.

Katsotaan, miten Maplen dsolve - komentoa käytettäisiin:

>

HY;

>	dsolve(HY,y(t));

2. EHY:n erikoinen:

>	r:=3*cos(3*t)-2*sin(3*t):

Koska ulkoisen voiman (herätteen) taajuus ei yhdy HY:n ratkaisun taajuuteen, niin "teoriasta" tiedämme jo etukäteen,

että samantaajuinen sin/cos-värähtely johtaa tulokseen. Ts. alla oleva yrite toimii.

>	yp:=a*cos(3*t)+b*sin(3*t);

>

vasen;

>	yriteyhtalo:=eval(subs(y(t)=yp,vasen))=r;

>	vas:=lhs(yriteyhtalo): oik:=rhs(yriteyhtalo):

cos-termien kertoimet samat, sin-termien kertoimet samat:

>	kerroinyhtalot:=coeff(vas,cos(3*t))=coeff(oik,cos(3*t)),coeff(vas,sin(3*t))=coeff(oik,sin(3*t));

>	ab:=solve({kerroinyhtalot},{a,b});

>	yp:=subs(ab,yp);

Tätä kutsutaan tasapainotilaratkaisuksi ("steady state") , pitkällä aikavälillä yhtälön ratkaisu on hyvin tarkkaan tämä

(mitä suurempi t, sen tarkempi), alkuarvoista riippumatta.

Kirjoitetaan muotoon: Vaihekulma ja amplitudi.

>	C:=sqrt(a^2+b^2);

>	assign(ab): C;

>	delta:=arctan(b,a);

>	C:=evalf(C); delta:=evalf(delta);

>	yss:=C*cos(3*t-delta);

Alkuehdot:

Muutetaan ratkaisulauseke funktioiksi, jotta alkuehdot olisi kätevä muodostaa.

>	Y:=unapply(yh+yp,t);

>	solve({Y(0)=1,D(Y)(0)=0},{A,B});

>	assign(%);

>

Y(t);

>

>	Y(0);D(Y)(0);

>	plot([r,Y(t),yp],t=0..15,color=[red,green,blue],title="Vaimennettu pakkovärähtely",legend=["Heräte (r(t))","Vaste (Y=yh+yp)","Tasapainotila (yp)"]);

>	plot([r,Y(t),yp],t=0..15,color=[red,green,blue],legend=[p,g,oo]);

>