Harj. 8 LV

to 9.11.00

1

> restart:with(linalg):with(inttrans):alias(L=laplace,IL=invlaplace):

Warning, the name changecoords has been redefined

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

Warning, the name hilbert has been redefined

> A:=matrix([[3,-2],[4,1]]);

A := matrix([[3, -2], [4, 1]])

> eAt:=exponential(A,t);

eAt := matrix([[exp(2*t)*cos(t*sqrt(7))+1/7*sqrt(7)...

> id:=diag(1,1):karmat:=inverse(s*id-A);

karmat := matrix([[(s-1)/(s^2-4*s+11), -2*1/(s^2-4*...

> map(L,eAt,t,s): map(simplify,%);

matrix([[(s-1)/(s^2-4*s+11), -2*1/(s^2-4*s+11)], [4...

Tuo henkeäsalpaava kaava siis toimii käytännössä. Tästä innostuneina käänteismuunnetaan vielä edellinen,

tulos ei enää yllätä.

> map(IL,%,s,t);

matrix([[exp(2*t)*cos(t*sqrt(7))+1/7*sqrt(7)*exp(2*...

>

> g:=t->vector([sin(t),-cos(t)]);

g := proc (t) options operator, arrow; vector([sin(...

> G:=map(L,g(t),t,s);

G := vector([1/(s^2+1), -s/(s^2+1)])

> x0:=vector([0,0]); X:=inverse(s*id-A)&*(x0+G);X:=evalm(%);

>

x0 := vector([0, 0])

X := `&*`(matrix([[(s-1)/(s^2-4*s+11), -2*1/(s^2-4*...

X := vector([(s-1)/((s^2-4*s+11)*(s^2+1))+2*s/((s^2...

> map(convert,X,parfrac,s);

vector([1/58*(-11+13*s)/(s^2+1)-1/58*(-63+13*s)/(s^...

> map(IL,%,s,t);

vector([13/58*cos(t)-11/58*sin(t)-13/58*exp(2*t)*co...

Näin saatiin ratkaisu matriisi-Laplace-muunnostekniikalla.

Verrataan vielä vakioiden varioimiskaavaan:

> map(int,evalm(exponential(A,t-tau)&*g(tau)),tau=0..t);

vector([13/58*cos(t)-11/58*sin(t)-13/58*exp(2*t)*co...

Kaikki tiet johtavat samaan, matematiikassa on taikaa ...

2

> restart: alias(Z=ztrans,IZ=invztrans):

Warning, the name changecoords has been redefined

> sum(1*z^(-k),k=0..infinity); # z-muunnoksen määritelmän mukaan.

z/(-1+z)

> series(%,z=infinity,10); # Käänteismuunnos on kerroinjono, kun kehitetään 1/z:n potenssien mukaan, ts. "äärettömyydessä".

1+1/z+1/(z^2)+1/(z^3)+1/(z^4)+1/(z^5)+1/(z^6)+1/(z^...

Siltäpä tosiaan näyttää.

Tehtäväpaperissa on sulut huonosti.

> mdz:=F->(-z*diff(F,z));

mdz := proc (F) options operator, arrow; -z*diff(F,...

> F:=z->z/(-1+z); # tässä siis 1:n z-muunnos

F := proc (z) options operator, arrow; z/(z-1) end ...

> a:=1/2: F(z/a):normal(%); # Kertolasku a^k:lla toimii:

2*z/(-1+2*z)

> Z(k*(1/2)^k,k,z)=simplify(mdz(F(z/a),z));

2*z/((-1+2*z)^2) = 2*z/((-1+2*z)^2)

> Z(k^2*(1/2)^k,k,z)=simplify((mdz@@2)(F(z/a),z));

2*z*(1+2*z)/((-1+2*z)^3) = 2*z*(1+2*z)/((-1+2*z)^3)...

> seq(Z(k^n*(1/2)^k,k,z)=simplify((mdz@@n)(F(z/a),z)),n=1..6);

2*z/((-1+2*z)^2) = 2*z/((-1+2*z)^2), 2*z*(1+2*z)/((...
2*z/((-1+2*z)^2) = 2*z/((-1+2*z)^2), 2*z*(1+2*z)/((...
2*z/((-1+2*z)^2) = 2*z/((-1+2*z)^2), 2*z*(1+2*z)/((...

> plot(2*z*(1+114*z+1208*z^2+2416*z^3+912*z^4+32*z^5)/((-1+2*z)^7),z=2..10);

[Maple Plot]

3

> restart: alias(Z=ztrans, IZ=invztrans):

Warning, the name changecoords has been redefined

a)

> X:=(z+2)/(z+1);X:=1+1/(z+1);normal(%): # tarkistus # Kirjoitettiin itse ajatellen.

X := (z+2)/(z+1)

X := 1+1/(z+1)

> x:=unapply(IZ(X,z,k),k); seq(x(k),k=0..10);

x := proc (k) options operator, arrow; 2*charfcn[0]...

1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1

> series(X,z=infinity,10);

1+1/z-1/(z^2)+1/(z^3)-1/(z^4)+1/(z^5)-1/(z^6)+1/(z^...

b)

> restart: alias(Z=ztrans, IZ=invztrans):

Warning, the name changecoords has been redefined

> X:=z/(z-1)/(z+2);

X := z/((z-1)*(z+2))

> X:=expand(z*convert(X/z,parfrac,z));

X := 1/3*z/(z-1)-1/3*z/(z+2)

> x:=unapply(IZ(X,z,k),k);

x := proc (k) options operator, arrow; 1/3-1/3*(-2)...

> seq(x(k),k=0..10);

0, 1, -1, 3, -5, 11, -21, 43, -85, 171, -341

> series(X,z=infinity,10);

1/z-1/(z^2)+3/(z^3)-5/(z^4)+11/(z^5)-21/(z^6)+43/(z...

c)

> restart: alias(Z=ztrans, IZ=invztrans):

Warning, the name changecoords has been redefined

> X:=z/(z^2+1);X:=expand(z*convert(X/z,parfrac,z,I));

X := z/(z^2+1)

X := -1/2*I*z/(z-I)+1/2*I*z/(z+I)

> IZ(X,z,k);evalc(%);x:=unapply(%,k);

-1/2*I*I^k+1/2*I*(-I)^k

sin(1/2*k*Pi)

x := proc (k) options operator, arrow; sin(1/2*k*Pi...

> seq(x(k),k=0..15); series(X,z=infinity,15);

0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1

1/z-1/(z^3)+1/(z^5)-1/(z^7)+1/(z^9)-1/(z^11)+1/(z^1...

d)

Tässä ei ole mitään laskemista:

tulos: (1, 0, 3, 0,0,0,0,0,0,-2,0,0, ...) eli

> charfcn[0](k)+3*charfcn[2](k)-2*charfcn[9](k);

charfcn[0](k)+3*charfcn[2](k)-2*charfcn[9](k)

> seq(%,k=0..10);

1, 0, 3, 0, 0, 0, 0, 0, 0, -2, 0

e)

> X:=(2*z^3+6*z^2+5*z+1)/((z^2)*(2*z+1));#convert(X,parfrax,z);P:=numer(X):Q:=denom(X);

X := (2*z^3+6*z^2+5*z+1)/(z^2*(2*z+1))

> X:=convert(X,parfrac,z);

X := 1+1/(z^2)+3/z-1/(2*z+1)

> invztrans(1/(2*z+1),z,k);

charfcn[0](k)-(-1/2)^k

Tästä eteenpäin on asian kannalta epärelevanttia.

Viimeksi ei osamurtoilu jostain syystä onnistunut, muuten kuin polynomijaon jälkeen. Syynä täytyi olla se, että jollekin muuttujalle oli jäänyt outo arvo. No, tulipa harjoitelluksi "pitkää" jakoa. Sinänsä on hyödyllistä tietää, että se sadaan aikaan Maplen quo ja rem -funktioilla:

> P:=numer(X);Q:=denom(X);osam:=quo(P,Q,z);jakoj:=rem(P,Q,z);

P := 2*z^3+6*z^2+5*z+1

Q := z^2*(2*z+1)

osam := 1

jakoj := 1+5*z^2+5*z

> expand(osam*Q+jakoj);

2*z^3+6*z^2+5*z+1

> X:=osam+jakoj/Q;normal(%);

X := 1+(1+5*z^2+5*z)/(z^2*(2*z+1))

(2*z^3+6*z^2+5*z+1)/(z^2*(2*z+1))

> convert(X,parfrac,z);series(%,z=infinity,10);

1+1/(z^2)+3/z-1/(2*z+1)

1+5/2/z+5/4/(z^2)-1/8*1/(z^3)+1/16/(z^4)-1/32*1/(z^...

> IZ(1/(1+2*z),z,k);

charfcn[0](k)-(-1/2)^k

> convert(X,parfrac,z);IZ(%,z,k);

1+1/(z^2)+3/z-1/(1+2*z)

charfcn[2](k)+3*charfcn[1](k)+(-1/2)^k

4. Käsitellään viikolla 46 (ajattele, että asuntolainasi kohtalo riippuu tästä)

5. Kompleksitapaus b) viimeistellään viikon 46 tehtävänä.

a)

> restart:alias(Z=ztrans,IZ=invztrans):

Warning, the name changecoords has been redefined

> dy:=6*y(k+2)+5*y(k+1)-y(k)=5;zdy:=Z(dy,k,z);

dy := 6*y(k+2)+5*y(k+1)-y(k) = 5

zdy := 6*z^2*Z(y(k),k,z)-6*y(0)*z^2-6*y(1)*z+5*z*Z(...

>

> Y:=solve(zdy,Z(y(k),k,z));y(0):=0:y(1):=0:

Y := z*(6*y(0)*z^2-y(0)*z+6*y(1)*z-6*y(1)-5*y(0)+5)...

> Y;

5*z/(6*z^3-z^2-6*z+1)

> Yperz:=convert(Y/z,parfrac,z);Y:=expand(z*Yperz);

Yperz := 1/2*1/(z-1)-36/7*1/(6*z-1)+5/14/(z+1)

Y := 1/2*z/(z-1)-36/7*z/(6*z-1)+5/14*z/(z+1)

> invztrans(Y,z,k);yy:=unapply(%,k);

1/2-6/7*(1/6)^k+5/14*(-1)^k

yy := proc (k) options operator, arrow; 1/2-6/7*(1/...

> plot([seq([k,yy(k)],k=0..8)]);

[Maple Plot]

>

> rsolve({6*y(k+2)+5*y(k+1)-y(k)=5,y(0)=0,y(1)=0},y(k));

5/14*(-1)^k-6/7*(1/6)^k+1/2

b)

> restart:alias(Z=ztrans,IZ=invztrans);

Z, IZ

> dy:=y(k+2)-3*y(k+1)+3*y(k)=1;

>

dy := y(k+2)-3*y(k+1)+3*y(k) = 1

> Zdy:=Z(dy,k,z);y(0):=1:y(1):=0:Zdy;

Zdy := z^2*Z(y(k),k,z)-y(0)*z^2-y(1)*z-3*z*Z(y(k),k...

z^2*Z(y(k),k,z)-z^2-3*z*Z(y(k),k,z)+3*z+3*Z(y(k),k,...

> Y:=solve(Zdy,Z(y(k),k,z));solve(z^2-3*z+3=0,z);Y:=expand(z*convert(Y/z,parfrac,z,{I,sqrt(3)}));

Y := z*(z^2-4*z+4)/(z^3-4*z^2+6*z-3)

3/2+1/2*I*sqrt(3), 3/2-1/2*I*sqrt(3)

Y := -2/3*I*z*sqrt(3)/(2*z-3+I*sqrt(3))+-2/3*I*z*sq...

> yy:=IZ(Y,z,k);

yy := -1/3*I*6^k*(1/(I*sqrt(3)+3))^k*sqrt(3)+1/3*I*...

> seq(evalf(yy),k=0..10);

1., -.4e-9+0.*I, -2.000000000+0.*I, -5.000000002+0....
1., -.4e-9+0.*I, -2.000000000+0.*I, -5.000000002+0....

> evalc(%);

-1/3*6^k*exp(-1/2*k*ln(12))*sin(1/6*k*Pi)*sqrt(3)-1...

Tässä ei tuo mainio evalc auta. On syytä laskea käsin (evalc ei esim. tunne "koulupojan kolmiota"). Edellä suoritttu seq näyttää, että reaalilukujono siitä sentään saadaan (ainakin laskentatarkkuuden puitteissa). On syytä varmistua asiasta kynän ja paperin avustuksella. Parhaiten se onnistunee yllä olevasta yy -muodosta, kun kirjoitetaan potenssien kanataluvut polaarimuotoon (ja käytetään tuota ennen mainittua - sanotaan nyt tasa-arvon nimissä - "koulutytön" kolmiota. Käsitellään tämän loppuunvienti ensi viikon (46) harj:ssa.

c)

6 Viimeistellään viikolla 46

> restart:alias(Z=ztrans,IZ=invztrans):

> dy1:=c(k+1)=1.5*c(k)-e(k);dy2:=e(k+1)=0.21*c(k)+0.5*e(k);

dy1 := c(k+1) = 1.5*c(k)-e(k)

dy2 := e(k+1) = .21*c(k)+.5*e(k)

> dys:={dy1,dy2};

dys := {c(k+1) = 1.5*c(k)-e(k), e(k+1) = .21*c(k)+....

> Zdys:=map(Z,dys,k,z);

Zdys := {z*Z(e(k),k,z)-e(0)*z = .21*Z(c(k),k,z)+.5*...

> Ratk:=solve(Zdys,{Z(c(k),k,z),Z(e(k),k,z)});

Ratk := {Z(c(k),k,z) = 12.50000000*z*(2.*c(0)*z-2.*...

> ratk:=map(IZ,Ratk,z,k);

ratk := {c(k) = 0, e(k) = 0}

> C:='C';E:='E':Sys:={(z-1.5)*C+E=z*c(0),-0.21*C+(z-0.5)*E=z*e(0)};

C := 'C'

Sys := {(z-1.5)*C+E = z*c(0), -.21*C+(z-.5)*E = z*e...

> solve(Sys,{C,E});

{E = .2500000000*z*(100.*z*e(0)+21.*c(0)-150.*e(0))...

> assign(%);

> Eperz:=E/z;

Eperz := .2500000000*(100.*z*e(0)+21.*c(0)-150.*e(0...

> expand(z*convert(Eperz,parfrac,z));

1.750000000*z*e(0)/(z-.8000000000)-.5250000000*z*c(...

> ee:=unapply(collect(IZ(%,z,k),{c(0),e(0)}),k);

ee := proc (k) options operator, arrow; (-.52500000...

> Cperz:=C/z:expand(z*convert(Cperz,parfrac,z)):cc:=unapply(collect(IZ(%,z,k),{c(0),e(0)}),k);

cc := proc (k) options operator, arrow; (-.75000000...

> C;

12.50000000*z*(2.*c(0)*z-2.*e(0)-1.*c(0))/(25.*z^2-...

> c(0):=c0:e(0):=e0:

> C;

12.50000000*z*(2.*c0*z-2.*e0-1.*c0)/(25.*z^2-50.*z+...

> 1.75*c0-2.5*e0;

2400.00

> c0:=12000: e0:=7440:

>