Tentti 3.2.2007

Ratkaisut / HA

                                      Tämä Maple-työarkki täydentää tiedostoa   helmi07ratkaisut.pdf

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

A := Matrix(%id = 20628100)

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

b := Vector(%id = 20624816)

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

Ab := Matrix(%id = 20908016)

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

rAb := Matrix(%id = 21505824)

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

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

ehto := b3+1-1/2*b2 = 0

(b)-kohta:

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

Ab := Matrix(%id = 21656376)

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

b2 := -2

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

Matrix(%id = 21656376)

Takaisinsijoitus:

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

14*x2+10*x3 = 2

x2 := -5/7*x3+1/7

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

x1+13/7*x3+3/7 = 1

x1 := -13/7*x3+4/7

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

matrix([[x[1] = -13/7*t+4/7], [x[2] = -5/7*t+1/7], [x[3] = t]])

>    #latex(%);

 Tehtävä 2.

(c)-kohta

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

P := Matrix(%id = 2741948)

Lasketaan ominaisarvoa 1 vastaava ominaisvektori:

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

A := Matrix(%id = 20755488)

>   

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

A := Matrix(%id = 20730048)

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

Suoritetaan Gauss-askeleet Maplella käsinlaskua matkien:

>    Ab;

Matrix(%id = 20730048)

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

Ab := Matrix(%id = 21060108)

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

Ab := Matrix(%id = 21347756)

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

Ab := Matrix(%id = 21512032)

>    x3:='x3':

>   

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

x2 := .9344262295*x3

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

x1 := .8384074941*x3

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

[8.384074941, 9.344262295, 10]

summa := 27.72833724

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

Vector(%id = 17744280)

(d)-kohtaan:

>    map(Re,Eigenvalues(P));

Vector(%id = 2387680)

>    0.904^40;

.1764969558e-1

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

A := Matrix(%id = 16804484)

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

lambda, ov := Vector(%id = 17515796), Matrix(%id = 17397152)

(a)

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

lambda := 3*I

w := Vector(%id = 17590576)

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

alpha := 0

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

beta := 3

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

u := Vector(%id = 17585736)

v := Vector(%id = 17333440)

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

uv := Matrix(%id = 17705004)

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

Kiertomatr := Matrix(%id = 16282288)

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

Y := Vector(%id = 18063884)

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

Vector(%id = 18031800) = Vector(%id = 17334240)

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

a := 1

b := 0

Vector(%id = 18063884)

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

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

Vector(%id = 16281160)

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

[Maple Plot]

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

Pisteet := [Vector(%id = 3032656), Vector(%id = 2989288), Vector(%id = 2989108), Vector(%id = 2988928)]

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

>   

suunnat := [Vector(%id = 17141288), Vector(%id = 17550860), Vector(%id = 17949792), Vector(%id = 18171944)]

[Vector(%id = 17141288), Vector(%id = 17550860), Vector(%id = 17949792), Vector(%id = 18171944)]

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

suunnat := [Vector(%id = 2468148), Vector(%id = 17336456), Vector(%id = 16310448), Vector(%id = 2625024)]

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

[Maple Plot]

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

[Maple Plot]

 Tehtävä 4.

>    restart:

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

F := proc (t, Y) options operator, arrow; `<,>`(Y[2],-3*Y[1]+t*Y[2]) end proc

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

Y[0] := Vector(%id = 17666156)

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

Y[1] := h*Vector(%id = 18334668)+Vector(%id = 17666156)

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

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

Ya := .2*Vector(%id = 18334668)+Vector(%id = 17666156)

vector([.4000000000, -3.600000000])

>    ya:=0.4;

>   

ya := .4

(b) Kaksi Euler-askelta, h=0.1

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

Yb[1] := .1*Vector(%id = 18334668)+Vector(%id = 17666156)

vector([.7000000000, -3.300000000])

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

Yb[2] := Vector(%id = 17939004)

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

yb := .369999999999999938

>    yb:=0.37;

yb := .37

>    yt:=0.3482;

yt := .3482

>    virhea:=yt-ya;

virhea := -.518e-1

>    virheb:=yt-yb;

virheb := -.218e-1

>    virheb/virhea;

.4208494208

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

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.

f := 50

>    bn;

-100*(cos(n*Pi)-1)/n/Pi

>    bn:=trigsiev(bn,n);

bn := -100*((-1)^n-1)/n/Pi

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

200/Pi, 0, 200/3/Pi, 0, 40/Pi, 0, 200/7/Pi, 0, 200/9/Pi, 0

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

lambda := proc (n) options operator, arrow; 1/10*n*Pi end proc

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

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

>   

Warning, the name changecoords has been redefined

>    u(x,t,3);

200/Pi*(sin(1/10*Pi*x)*exp(-1/100*Pi^2*t)+1/3*sin(3/10*Pi*x)*exp(-9/100*Pi^2*t)+1/5*sin(1/2*Pi*x)*exp(-1/4*Pi^2*t))

Piirretään pintakuva.

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

[Maple Plot]

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

[Maple Plot]

Sarjan 1. termi:

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

termi1 := 200/Pi*sin(1/10*Pi*x)*exp(-1/100*Pi^2*t)

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

maxlampotila := 200/Pi*sin(1/2*Pi)*exp(-1/100*Pi^2*t)

>    ehto:=maxlampotila=10;

ehto := 200/Pi*exp(-1/100*Pi^2*t) = 10

>    t1:=solve(ehto,t);

t1 := -100*ln(1/20*Pi)/Pi^2

Ihminen kirjoittaisi tuon plus-merkkisenä vaihtaen osoittajan muotoon    ln(20/Pi)  .

>    evalf(t1);

18.75457528

Siinäpä se!