Tätä vielä kehitellään ... -------------------------- > # harj. 5 loppuviikko > # > # teht. 1 > with(DEtools):with(linalg): Warning, new definition for norm Warning, new definition for trace > dyht:=diff(y(t),t)=sin(t)-y(t)/2; # suositus: kirjoita y(t) joka paikkaan # yhtälössä (tehtäväpaperissa erehdyksessä # vastoin suositusta). d dyht := -- y(t) = sin(t) - 1/2 y(t) dt > DEplot(dyht,[y],t=0..12, {[0,0],[0,1],[0,-1]}, y=-3..3, arrows=THIN); > # versio V.3 todellakin vaatii muodon DEplot(dyht,{t,y], ... > # 1-tehtävän piirroksia: > dyht:= > dyvp:=diff(x(t),t$2)+4*diff(x(t),t)+13*x(t); / 2 \ |d | /d \ dyvp := |--- x(t)| + 4 |-- x(t)| + 13 x(t) | 2 | \dt / \dt / > dy:=dyvp=0; / 2 \ |d | /d \ dy := |--- x(t)| + 4 |-- x(t)| + 13 x(t) = 0 | 2 | \dt / \dt / > dsolve({dyvp=cos(t),x(0)=0,D(x)(0)=1},x(t)); x(t) = 1/30 cos(3 t) cos(4 t) - 1/60 cos(3 t) sin(4 t) + 1/24 cos(3 t) cos(2 t) - 1/24 cos(3 t) sin(2 t) 4 3 + 2/15 sin(3 t) cos(t) + 4/15 sin(3 t) cos(t) sin(t) 2 - 1/20 sin(3 t) cos(t) - 1/20 sin(3 t) cos(t) sin(t) - 1/40 sin(3 t) - 3/40 exp(-2 t) cos(3 t) 11 + -- exp(-2 t) sin(3 t) 40 > # kauhea sotku, siivotaan vähän: > combine("); x(t) = 3/40 cos(t) + 1/40 sin(t) - 3/40 exp(-2 t) cos(3 t) 11 + -- exp(-2 t) sin(3 t) 40 > # saatiin oikea vastaus! > ##################################### > # teht.2 > A4a:=matrix(2,2,[3,-1,-1,3]);A4b:=matrix(2,2,[0,-4,1,-5]); [ 3 -1] A4a := [ ] [-1 3] [0 -4] A4b := [ ] [1 -5] > A5a:=matrix(2,2,[-3,1,-4,2]); A5b:=matrix(2,2,[-2,2,-2,-2]); [-3 1] A5a := [ ] [-4 2] [-2 2] A5b := [ ] [-2 -2] > eigenvals(A4a); 2, 4 > # reaaliset, positiiviset => lähde. Piirto: > ?DEplot > yhtop:=evalm(A4a &* vector([y1(t),y2(t)])); yhtop := [3 y1(t) - y2(t), -y1(t) + 3 y2(t)] > dy4a:={D(y1)(t)=yhtop[1],D(y2)(t)=yhtop[2]}; dy4a := {D(y1)(t) = 3 y1(t) - y2(t), D(y2)(t) = -y1(t) + 3 y2(t)} > type(dy4a,set); true > dsolve(dy4a,{y1(t),y2(t)}); {y1(t) = 1/2 _C1 exp(2 t) + 1/2 _C1 exp(4 t) - 1/2 _C2 exp(4 t) + 1/2 _C2 exp(2 t), y2(t) = - 1/2 _C1 exp(4 t) + 1/2 _C1 exp(2 t) + 1/2 _C2 exp(2 t) + 1/2 _C2 exp(4 t)} > # Versio 4 antaa arrows-komennolle vaihtoehdot SMALL,MEDIUM,LARGE,LINE,NONE > DEplot(dy4a,[y1(t),y2(t)],t=0..3,[[y1(1)=-1,y2(1)=-.5],[y1(1)=0,y2(1)=1],[y1(1)=2,y2(1)=1],[y1(1)=1,y2(1)=1],[y1(1)=-1,y2(1)=-2]],y1=-5..5,y2=-5..5,arrows=MEDIUM); > # yllä laitettiin alkuehtojen evaluointipisteeksi t=1 (siis y1(1)=..., y2(1)=...) jotta käyrä piirtyisi myös taaksepäin alkuehtopisteestä. > eigenvals(A4b); -4, -1 > # reaaliset, negatiiviset => nielu. Piirto: > yhtop:=evalm(A4b &* vector([y1(t),y2(t)])); yhtop := [-4 y2(t), y1(t) - 5 y2(t)] > dy4b:={D(y1)(t)=yhtop[1],D(y2)(t)=yhtop[2]}; dy4b := {D(y2)(t) = y1(t) - 5 y2(t), D(y1)(t) = -4 y2(t)} > DEplot(dy4b,[y1(t),y2(t)],t=0..3,[[y1(1)=2,y2(1)=-2],[y1(1)=0,y2(1)=2],[y1(1)=5,y2(1)=-1],[y1(1)=-3,y2(1)=3],[y1(1)=-1,y2(1)=-2]],y1=-5..5,y2=-5..5,arrows=MEDIUM); > eigenvals(A5a); -2, 1 > # reaaliset, erimerkkiset => satula. Piirto: > yhtop:=evalm(A5a &* vector([y1(t),y2(t)])); yhtop := [-3 y1(t) + y2(t), -4 y1(t) + 2 y2(t)] > dy5a:={D(y1)(t)=yhtop[1],D(y2)(t)=yhtop[2]}; dy5a := {D(y2)(t) = -4 y1(t) + 2 y2(t), D(y1)(t) = -3 y1(t) + y2(t)} > DEplot(dy5a,[y1(t),y2(t)],t=0..3,[[y1(1)=2,y2(1)=-2],[y1(1)=0,y2(1)=2],[y1(1)=1,y2(1)=1.7],[y1(1)=-1,y2(1)=-2],[y1(1)=1,y2(1)=1]],y1=-5..5,y2=-5..5,arrows=MEDIUM); > eigenvals(A5b); -2 + 2 I, -2 - 2 I > # kompleksiset, reaaliosa negatiivinen => stabiili fokus, eli spiraali joka suppenee kohti origoa. Piirto: > yhtop:=evalm(A5b &* vector([y1(t),y2(t)])); yhtop := [-2 y1(t) + 2 y2(t), -2 y1(t) - 2 y2(t)] > dy5b:={D(y1)(t)=yhtop[1],D(y2)(t)=yhtop[2]}; dy5b := {D(y1)(t) = -2 y1(t) + 2 y2(t), D(y2)(t) = -2 y1(t) - 2 y2(t)} > DEplot(dy5b,[y1(t),y2(t)],t=0..3,[[y1(1)=5,y2(1)=2],[y1(1)=0,y2(1)=2],[y1(1)=5,y2(1)=-1],[y1(1)=3,y2(1)=3],[y1(1)=-2,y2(1)=1]],y1=-5..5,y2=-5..5,arrows=MEDIUM); > ################################ > # teht. 3 > dy:={D(y1)(t)=-2*y1(t)+2*y2(t)+2*E,D(y2)(t)=1/5*(-2*y1(t)+y2(t)+2*E)}; dy := {D(y1)(t) = -2 y1(t) + 2 y2(t) + 2 E, D(y2)(t) = - 2/5 y1(t) + 1/5 y2(t) + 2/5 E} > dy1:=subs(E=220*sin(t),dy); dy1 := {D(y1)(t) = -2 y1(t) + 2 y2(t) + 440 sin(t), D(y2)(t) = - 2/5 y1(t) + 1/5 y2(t) + 88 sin(t)} > dsolve(dy1,{y1(t),y2(t)}); 11 1/2 11 1/2 {y1(t) = 1/2 _C1 %2 + -- _C1 41 %1 - -- _C1 41 %2 + 1/2 _C1 %1 82 82 10 1/2 10 1/2 - -- _C2 41 %1 + -- _C2 41 %2 + 616/3 sin(t) 41 41 1/2 1/2 - 352/3 cos(t), y2(t) = 2/41 _C1 41 %1 - 2/41 _C1 41 %2 11 1/2 11 1/2 + 1/2 _C2 %2 - -- _C2 41 %1 + -- _C2 41 %2 + 1/2 _C2 %1 82 82 + 44 sin(t) - 44/3 cos(t)} 1/2 %1 := exp(- 1/10 (9 + 41 ) t) 1/2 %2 := exp(1/10 (-9 + 41 ) t) > # kokeillaan luonnollisilla alkuehdoilla y1(0)=0=y2(0): > ratk:=dsolve({op(dy1),y1(0)=0,y2(0)=0},{y1(t),y2(t)}); 154 1/2 154 1/2 ratk := {y2(t) = 22/3 %1 + 22/3 %2 + --- 41 %2 - --- 41 %1 41 41 + 44 sin(t) - 44/3 cos(t), y1(t) = 176/3 %1 + 176/3 %2 1496 1/2 1496 1/2 + ---- 41 %2 - ---- 41 %1 + 616/3 sin(t) - 352/3 cos(t)} 123 123 1/2 %1 := exp(1/10 (-9 + 41 ) t) 1/2 %2 := exp(- 1/10 (9 + 41 ) t) > ############################ > # teht. 4 > A:=matrix(2,2,[-1,1,-1,-1]); [-1 1] A := [ ] [-1 -1] > yhtop:=evalm(A &* vector([y1(t),y2(t)])); yhtop := [-y1(t) + y2(t), -y1(t) - y2(t)] > dy:={D(y1)(t)=yhtop[1],D(y2)(t)=yhtop[2]}; dy := {D(y2)(t) = -y1(t) - y2(t), D(y1)(t) = -y1(t) + y2(t)} > eigenvals(A); -1 + I, -1 - I > dsolve(dy,{y1(t),y2(t)}); {y2(t) = -exp(-t) sin(t) _C1 + exp(-t) cos(t) _C2, y1(t) = exp(-t) cos(t) _C1 + exp(-t) sin(t) _C2} > [[y1(1)=5,y2(1)=2],[y1(1)=0,y2(1)=2],[y1(1)=5,y2(1)=-1],[y1(1)=3,y2(1)=3],[y1(1)=-2,y2(1)=1]], > DEplot(dy,[y1(t),y2(t)],t=0..5,[[y1(1)=0,y2(1)=-3],[y1(1)=0,y2(1)=3]],y1=-5..5,y2=-5..5,arrows=MEDIUM); > ######################## > # teht. 5 > ?exp > ?linalg[exponential] > A:=matrix(2,2,[1,1,0,-1]); f:= t -> vector([t,1]); [1 1] A := [ ] [0 -1] f := t -> vector([t, 1]) > eksfkt:= t -> exponential(A*t); eksfkt := t -> exponential(A t) > ?int > map(int,f(x),x=0..t); # pakko ajaa map:n läpi. [ 2 ] [1/2 t , t] > x1:=evalm(eksfkt(t) &* x0); [exp(t) 1/2 exp(t) - 1/2 exp(-t)] x1 := [ ] &* x0 [ 0 exp(-t) ] > x2:=map(int,value(evalm(eksfkt(t-u) &* f(u))), u=0..t); x2 := [-t - 2 + (3/2 exp(2 t) + 1/2) exp(-t), 1 - exp(-t)] > x:=x1+x2; /[exp(t) 1/2 exp(t) - 1/2 exp(-t)] \ x := |[ ] &* x0| + x2 \[ 0 exp(-t) ] / > # kokeillaan jollain alkuehdolla: > subs(x0=vector([1,1]),x); /[exp(t) 1/2 exp(t) - 1/2 exp(-t)] \ |[ ] &* [1, 1]| + x2 \[ 0 exp(-t) ] / > value(evalm(")); [3/2 exp(t) - 1/2 exp(-t) - t - 2 + (3/2 exp(2 t) + 1/2) exp(-t), 1] > # onnistui. > ############################## > # teht. 6 >