Alustukset

>

restart:

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

>	pyorista:=(x,nn)->evalf(round(10^nn*x)/10^nn):

>	%with(plots):

 Tehtävä 1.

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

>	b:=<1,b2,b3>;

>	Ab:=<A | b>;   # Liitännäismatriisi

>	rAb:=ref(Ab);       # Gaussataan ref-muotoon

Jotta ratkaisuja on, jos ja vain jos viimeinen sarake ei ole tukisarake, ts.

>	ehto:=b3+1-b2/2 = 0;  # Kohdan (a) vastaus:

(b)-kohta:

>	Ab:=subs(b3=-2,rAb);   # Sijoitetaan b3:lle arvo -2 Gaussatussa liitännäismatriisin.

>	b3:=-2: b2:=solve(ehto,b2);  # Tietysti nähdään heti suoraan yllä olevastakin.

>	Ab;  # Ab:ssä sioitettu b2:=-2;

Takaisinsijoitus:

>	14*x2+10*x3 = 2; x2:=solve(%,x2);

>	x1+3*x2+4*x3 = 1; x1:=solve(%,x1);

Jos merkitään x3 = t = vapaa,  saadaan yleinen ratkaisu muotoon

>	array([[x[1] = -13*t/7 +1-3/7],[x[2] = -5*t/7 + 1/7], [x[3] = t]]);

>	#latex(%);

 Tehtävä 2.

(c)-kohta

>	P:=<<.93,.03,.04>|<.02,.93,.05>|<.04,.04,.92>>;

Lasketaan ominaisarvoa 1 vastaava ominaisvektori:

>	A:=P-Diag([1,1,1]);

>

>	A:=map(pyorista,A,3);

>	Ab:=A:  # "b"-sarake on 0, joten sitä ei tarvita.

Suoritetaan Gauss-askeleet Maplella käsinlaskua matkien:

>

Ab;

>	kerroin:=Ab[2,1]/Ab[1,1]: Ab[2,1..-1]:=Ab[2,1..-1]-kerroin*Ab[1,1..-1]:Ab:=map(pyorista,Ab,3);

>	kerroin:=Ab[3,1]/Ab[1,1]: Ab[3,1..-1]:=Ab[3,1..-1]-kerroin*Ab[1,1..-1]: Ab:=map(pyorista,Ab,3);

>	kerroin:=Ab[3,2]/Ab[2,2]: Ab[3,1..-1]:=Ab[3,1..-1]-kerroin*Ab[2,1..-1]: Ab:=map(pyorista,Ab,3);

>

x3:='x3':

>

>	x2:='x2': -0.061*x2+0.057*x3=0: x2:=solve(%,x2);latex(%,"huu.tex",append);

>	x1:='x1': -0.07*x1 + 0.02*x2 + 0.04*x3 =0: x1:=solve(%,x1);latex(%,"huu.tex",append);

>	x3:=10: [x1,x2,x3];summa:=sum(%[i],i=1..3);

>	100*<x1,x2,x3>/summa: map(pyorista,%,2);

(d)-kohtaan:

>	map(Re,Eigenvalues(P));

>

0.904^40;

 Tehtävä  3.

>

restart:

>	with(LinearAlgebra):with(plots):#with(plottools):

Warning, the name changecoords has been redefined

>	#arrow(<1,0>,<0,1>);

>	A:=<<0,-9>|<1,0>>;

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

(a)

>	lambda:=3*I; w:=<1,3*I>;

>	#A.w = lambda*w; # Tarkistetaan, että on tosiaan ominaisarvo/-vektori

>

>	# Kompleksinen muoto, jota ei tässä käytetä:

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

Mennään suoraan reaalisen muotoon, kaavahan on annettu.

>	alpha:=Re(lambda);    # Ominaisarvon reaaliosa

>	beta:=Im(lambda);     # Ominaisarvon imag-osa

>	u:=map(Re,w); v:=map(Im,w);  # Omiaisvektorin Re- ja Im-osavektorit

>	uv:=<u|v>;           # Ladotaan vierekkäin

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

>	Y:=exp(alpha*t)*(uv.Kiertomatr.<a,b>);

>	subs(t=0,Y)=<1,0>; # Alkuarvo y(0)=[1,0]

>	a:=1; b:=0; Y; # Sijoitetaan ratkaisut a:=1 ja b:=0 ja katsotaan Y.

Trajektori on ellipsi, jonka puoliakselit ovat 1 ja 3, kierretään myötäpäivään, kun t kasvaa.

>	subs(t=Pi/6,Y);

>	plot([Y[1],Y[2],t=0..3*Pi/3],scaling=constrained):kuva:=%:

>	display(kuva,plots[arrow](<1,0>,-0.5*<0,1>,width=1/16));

>	Pisteet:=[<0,-3>,<-1,0>,<0,3>,<1,0>];

>	suunnat:=map(p->A.p,Pisteet);suunnat;

>

>	suunnat:=map(v->1/5*v,suunnat);

>	display(kuva,seq(plots[arrow](Pisteet[i],suunnat[i]),i=1..4));

>	parvi:=seq(seq([Y[1],Y[2],t=0..2],a=-2..2),b=-2..2):

>	plot([parvi],scaling=constrained):

>

kuva:=%:

>	with(DEtools):

>	display(kuva,DEplot([D(y1)(t)=y2(t),D(y2)(t)=-9*y1(t)], [y1(t),y2(t)],t=0..2,y1=-3..3,y2=-9..9,color=grey,stepsize=.02),title="",thickness=1,labels=["y1","y2"]);

 Tehtävä 4.

>

restart:

>	F:=(t,Y)-><Y[2],-3*Y[1]+t*Y[2]>;

>	Y[0]:=<1,-3>;

>	Y[1]:=Y[0]+h*F(0,Y[0]);

(a) Yksi Euler-askel: h=0.2.

>	Ya:=subs(h=0.2,Y[1]); evalm(%);

>

ya:=0.4;

>

(b) Kaksi Euler-askelta, h=0.1

>	Yb[1]:=subs(h=0.1,Y[1]);(evalm(%));

>	Yb[2]:=Yb[1]+0.1*F(0.1,Yb[1]);

>	yb:=Yb[2][1];

>

yb:=0.37;

>	yt:=0.3482;

>	virhea:=yt-ya;

>	virheb:=yt-yb;

>	virheb/virhea;

Tehtävä 5

Tehtävä on (aivan tarkoituksella) aivan samanlainen kuin joulukuun 2006 tentissä teht. 5. Pari numeroarvoa vain on muutettu.

Tässä tyydytään suorittamaan lasku Maplella.

Ensin alustus ja latauskomentoja ("ignore").

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

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

Tarvitsemme sinisarjaa.

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

>	L:=10: f:=x->50; # Kehitettävä vakiofunktio f=50 välillä [0,10], jakson puolikas L=10.

>

bn;

>	bn:=trigsiev(bn,n);

>	seq(bn,n=1..10);

>	lambda:=n->n*Pi/10;

>	Sinisarjan osasummat:

>	u:=(x,t,N)->200/Pi*sum(1/(2*k-1)*sin((2*k-1)*Pi*x/10)*exp((-lambda(2*k-1)^2)*t),k=1..N);

>

Warning, the name changecoords has been redefined

>

u(x,t,3);

Piirretään pintakuva.

>	plot3d(u(x,t,20),x=0..10,t=0..20);

Lisäksi lämpötilaprofiilit 0.5 s välein. Alin on ajanhetkellä t=20, nähdään, että silloin maxlämpötila

menee jo hiukan alle 10 asteen.

>	plot([seq(u(x,t,20),t=[seq(j*0.5,j=0..40)])],x=0..10);

Sarjan 1. termi:

>	termi1:=u(x,t,1);

>	maxlampotila:=subs(x=5,termi1); # maximi on keskipisteessä.

>	ehto:=maxlampotila=10;

>	t1:=solve(ehto,t);

Ihminen kirjoittaisi tuon plus-merkkisenä vaihtaen osoittajan muotoon     .

>	evalf(t1);

Siinäpä se!