L5maple.mws, Differenssiyhtälöitä

to 26.9.02 , ti -ke  1-2.10.02   Luennot 9-11

Lähteet

Lay  CH 4.8, s. 277 - 285    (prujut)

Lay Ch 4.9 ss. 288 - 296. (ei prujuja)

Alustukset, funktiomääritykset.

Keräämme funktiomäärityksiä tiedostoon v302.mpl. Toistaiseksi sijoittelemme niitä työarkeille tarpeen mukaan.

Tekstitiedostoston voi lukea komennolla

>    # read("v302.mpl"): # Määrittelyt (ja muut komennot) suorittuvat.

>    with(plots): 

Warning, the name changecoords has been redefined

>    jana:=(p1,p2)->plot([p1,p2],thickness=2,color=green):  # Hyödyllinen funktio diskreeteissä visualisoinneissa. Jos haluat muuttaa väriä, paksuutta, tms., modifioi vain tätä.

>    yjana:=(x,y)->jana([x,0],[x,y]); # Jana pisteestä (x,0) pisteeseen (x,y).

yjana := proc (x, y) options operator, arrow; jana([x, 0],[x, y]) end proc

Testataan:

>    f:=x->x^2-2*sqrt(x):

>    h:=.5: display(seq(jana([i*h,0],[i*h,f(i*h)]),i=0..10)):  # Kätevä.

>    h:=.5: display(seq(yjana(i*h,f(i*h)),i=0..10));  # kätevämpi.

[Maple Plot]

>    f:=k->.7^k;

f := proc (k) options operator, arrow; .7^k end proc

>    h:=1: display(seq(yjana(i*h,f(i*h)),i=-3..8));

[Maple Plot]

1. Johdattelua

- Diskreetti, eli digitaalinen data lisääntyvässä määrin

- Myös jatkuvia ilmiöitä lasketaan yleensä differenssiyhtälöillä, numeeriset menetelmät toimivat yleensä siten,

    että esim. differentiaaliyhtälö diskretoidaan, jolloin varsinainen ratkaisu  tapahtuu differenssiyhtälömaailmassa.

2. Diskreettiaikasignaalit

Diskreetti signaali on funktio, jonka määrittelyjoukko on Z . Ts. signaali on jonoavaruuden S  alkio (vektori).

Diskreetti dignaali voi syntyä prosessista luonnostaan tai sillä kuvataan jatkuvaa prosessia ottamalla esim. tasavälisin aikavälein näytteitä.

Esim: 1 : CD-levyllä oleva digitaalinen data edustaa musiikin amplitudiarvoja, kun näytteitä on otettu 44 100 kertaa sekunnissa.

Yllä määritellyt funktiot, erityisesti yjana ovat käteviä diskreettien signaalien esittämiseen.

>    f:=t->sin(t)+2*sin(2*t)-0.5*cos(7*t):

>    h:=0.5: display(seq(yjana(i*h,f(i*h)),i=0..20));

[Maple Plot]

>    display(plot(f(t),t=0..20),seq(yjana(i*h,f(i*h)),i=0..40));

[Maple Plot]

>   

3. Lineaarinen differenssiyhtälö

Lay s. 280

Exa 3. Lineaarinen suodin (filtteri)

>    restart: with(plots): 

Warning, the name changecoords has been redefined

>    jana:=(p1,p2)->plot([p1,p2],color=blue):
yjana:=(x,y)->jana([x,0],[x,y]):

>    F:=(y,k)->a0*y(k+2)+a1*y(k+1)+a2*y(k);

>   

F := proc (y, k) options operator, arrow; a0*y(k+2)+a1*y(k+1)+a2*y(k) end proc

>    a0:=sqrt(2)/4;

a0 := 1/4*2^(1/2)

>    a1:=1/2: a2:=a0:

>    F(y,k)=z(k);

1/4*2^(1/2)*y(k+2)+1/2*y(k+1)+1/4*2^(1/2)*y(k) = z(k)

Tässä ajatellaan y(k) syötteeksi ja z(k) tulokseksi. (Differenssiyhtälöitä ratkaistaessa päinvastoin)

Valitaan ensin signaaliksi y(k)  näytteitä kokonaislukuvälein funktiosta cos(Pi*t/4).

>    Y:=k->cos(Pi*k/4);

Y := proc (k) options operator, arrow; cos(1/4*Pi*k) end proc

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

1, 1/2*2^(1/2), 0, -1/2*2^(1/2), -1, -1/2*2^(1/2), 0, 1/2*2^(1/2), 1, 1/2*2^(1/2), 0

>    seq(F(Y,k),k=0..10);

1/2*2^(1/2), 0, -1/2*2^(1/2), -1, -1/2*2^(1/2), 0, 1/2*2^(1/2), 1, 1/2*2^(1/2), 0, -1/2*2^(1/2)

>    display(plot(Y(t),t=0..20),seq(yjana(k,Y(k)),k=0..40));

[Maple Plot]

>    display(plot(Y(t),t=0..20),seq(yjana(k,F(Y,k)),k=0..40));

[Maple Plot]

Signaali (Y(k)) menee filtterin läpi sellaisenaan muuten, paitsi siirtyy askeleen taaksepäin..

Pidetään sama näytteidenottoväli, kasvatetaan alkup. signaalin taajuutta:

>    Y:=k->cos(3*Pi*k/4);

Y := proc (k) options operator, arrow; cos(3/4*Pi*k) end proc

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

1, -1/2*2^(1/2), 0, 1/2*2^(1/2), -1, 1/2*2^(1/2), 0, -1/2*2^(1/2), 1, -1/2*2^(1/2), 0

>    seq(F(Y,k),k=0..10);

0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0

Kolminkertainen taajuus suodattui pois (näytteenottoväli pysyi samana).

Alipäästösuodin .

>    display(plot(Y(t),t=0..10),seq(yjana(k,Y(k)),k=0..20));

[Maple Plot]

>    display(plot(Y(t),t=0..10,axes=boxed),seq(yjana(k,F(Y,k)),k=0..20));

[Maple Plot]

>   

Exa 4. Lineaarinen homogeeninen differenssiyhtälö

>    restart:

>    differy:=y(k+3)-2*y(k+2)-5*y(k+1)+6*y(k)=0;

differy := y(k+3)-2*y(k+2)-5*y(k+1)+6*y(k) = 0

Yrite:

>    Y:=k->r^k:

>    subs(y=Y,differy);

Y(k+3)-2*Y(k+2)-5*Y(k+1)+6*Y(k) = 0

>    vdiffery:=value(%);

vdiffery := r^(k+3)-2*r^(k+2)-5*r^(k+1)+6*r^k = 0

>    vdiffery/r^k;simplify(%,symbolic);

1/(r^k)*(r^(k+3)-2*r^(k+2)-5*r^(k+1)+6*r^k) = 0

r^3-2*r^2-5*r+6 = 0

>    factor(%);

(r-1)*(r-3)*(r+2) = 0

Siis jonot u,v,w toteuttavat differy:n, kun

>    u:=k->1^k; v:=k->3^k; w:=k->(-2)^k;

u := 1

v := proc (k) options operator, arrow; 3^k end proc

w := proc (k) options operator, arrow; (-2)^k end proc

>    subs(y=u,differy);

1(k+3)-2*1(k+2)-5*1(k+1)+6*1(k) = 0

>    value(%);

0 = 0

>    subs(y=v,differy);

v(k+3)-2*v(k+2)-5*v(k+1)+6*v(k) = 0

>    value(%);

3^(k+3)-2*3^(k+2)-5*3^(k+1)+6*3^k = 0

>    simplify(%,symbolic);

0 = 0

>    subs(y=w,differy); value(%);simplify(%,symbolic);

w(k+3)-2*w(k+2)-5*w(k+1)+6*w(k) = 0

(-2)^(k+3)-2*(-2)^(k+2)-5*(-2)^(k+1)+6*(-2)^k = 0

0 = 0

>   

Siis saatiin 3 ratkaisua u,v,w.

Ovatko LRT? Ovatko ratkaisuavaruuden kanta?

Kyllä, lauseen 16 mukaan!

Exa 5. Lineaarinen EHY

>    restart:
HY:=y(k+2)-4*y(k+1)+3*y(k)=0;

HY := y(k+2)-4*y(k+1)+3*y(k) = 0

>    karpoly:=r^2-4*r+3: factor(%);

(r-1)*(r-3)

HY:n yleinen:

>    yh:=k->c1+c2*3^k;

yh := proc (k) options operator, arrow; c1+c2*3^k end proc

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

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

>   

Etsitään EHY:n erityistä. Luonnollinen lähtökohta on k:n polynomi, kun oikealla on sellainen (1. asteen).

a+bk - muodossa a on turha, koska pelkkä vakio toteuttaa HY:n.  

>    yrite:=k->b*k:

>    EHYyr:=subs(y=yrite,EHY);

EHYyr := yrite(k+2)-4*yrite(k+1)+3*yrite(k) = -4*k

>    eval(EHYyr);

b*(k+2)-4*b*(k+1)+3*b*k = -4*k

>    map(collect,%,b);

-2*b = -4*k

1. asteen polynomi ei siis toteuta. Yritetään 2. asteen (vakiotermi on edelleenkin turha).

>    yrite:=k->a*k^2+b*k:

>    EHYyr:=subs(y=yrite,EHY);

EHYyr := yrite(k+2)-4*yrite(k+1)+3*yrite(k) = -4*k

>    eval(%);

a*(k+2)^2+b*(k+2)-4*a*(k+1)^2-4*b*(k+1)+3*a*k^2+3*b*k = -4*k

>    expand(%);

-4*a*k-2*b = -4*k

Nähdään, että b=0, a=1, joten olemme löytäneet EHY:n erikoisen

>    yp:=k->k^2;

yp := proc (k) options operator, arrow; k^2 end proc

EHY:n yleinen ratkaisu on siis:

>    Y:=yh+yp;

Y := yh+yp

>    Y(k);

c1+c2*3^k+k^2

>   

Yleisiä lauseita ja periaatteita

>    restart:

(EHY): Epähomogeeniyhtälö, (HY) Homogeeniyhtälö (oikea puoli = 0).

Lineaarikuvaus T:S--> S,

>    T(({y[k]}))={y[k+n]+sum(a[j]*y[k+n-j],j=1..n)}[k];

T({y[k]}) = {y[k+n]+sum(a[j]*y[k+n-j],j = 1 .. n)}[k]

>    H=N(T); # T:n nolla-avaruus (ydin) on (EHY):n ratkaisujoukko.

H = N(T)

Lause 16  Lineaarisella n:nnen kertaluvun differyllä (a[n] # 0)

>    y[k+n]+sum(a[j]*y[k+n-j],j=1..n)=z[k];

y[k+n]+sum(a[j]*y[k+n-j],j = 1 .. n) = z[k]

on jokaisella alkuarvojevektorin   [ y[0] , ..., y[n-1] ]  valinnalla yksikäsitteinen ratkaisu ( y[k] ).

Tod: Rekursiokaavasta voidaan yksikäsitteisesti laskea sekä eteenpäin että taaksepäin (jos tarpeen) miten pitkälle tahansa.

(Tarkemmin luennolla tai kirjassa(prujussa).)

Lause 17  

Lineaarisen n:nnen kl:n  (HY):n ratkaisuavaruuden dimensio = n.

  Tod : H=N(T) on aliavaruus. Osoitetaan, että dim(H)=n. Sitä varten muodostetaan lineaarikuvaus F: H -> R^n  , joka osoitetaan

 isomorfismiksi. Isomorfismi säilyttää kaikki VA-ominaisuudet, erityisesti dimension.

No, määritellään "projektiokuvaus"   F: ( y[k] )  -->( y[0] , ... , y[n-1] )

Osoitetaan, että F: H --> R^n   on isomorfismi. ts. lineaarinen bijektio.

1) Lineaarisuus on hyvin selvä.

2) Surjektio & injektio: Olkoon ( y[0] , . . ., y[n-1] )  mielivaltainen R^n :n vektori. Lauseen 16 mukaan on olemassa 1-käs. ratkaisu ( y[k] ),

jonka koordinaatit 0,...,n-1 ovat nämä annetut. Olemassaolo tarkoittaa F:n surjektiivisuutta, yksikäsitteisyys injektiivisyyttä.

No niin, kuten sanottu, dim H = dim R^n  = n.   [QED]

HY:n ratkaisumenetelmä:  Etsitään n LRT ratkaisua.

Tämä on ratkaisuavaruuden kanta , joten kaikki ratkaisut ovat näiden lineaarikombinaatioia.

>   

yh[k] = sum(c[j]*y[j][k],j = 1 .. n)

Miten löydetään? Yrite

>   

y[k] = r^k

johtaa karakteristiseen polynomiyhtälöön ("auxiliary equation"):

>   

r^n+sum(a[j]*r^(n-j),j = 1 .. n) = 0

  • Tapaukset:
  • Yksinkertaiset reaalijuuret (eli n erillstä juurta). Helppo osoittaa, että tällöin ratkaisut (r[i]^k) ovat LRT ja siis ratkaisuavaruuden kanta. (Harj 4 AV)
  • Yksinkertaiset kompleksijuuret. Tällöin saadaan myös kanta aivan suoraan. Eulerin kaavan avulla siirrytään reaaliseen kantaan, jossa on sin/cos-termejä.
  • Moninkertaisia juuria. Tällöin tarvitaan a) temppuja tai b) teoriaa. Palataan ominaisarvojen yhteydessä.

EHY:n ratkaisumenetelmä: HY:n yleinen + EHY:n erikoinen

Lause. (EHY):n kaikki ratkaisut ovat muotoa

>   

y = y[h]+y[p]

missä y[h] on HY:n yleinen ratkaisu (siis lineaarikombinaatio n:stä kantaratkaisusta) ja y[p] on jokin EHY:n (vaikkapa hatusta tempaistu) "erityis"ratkaisu.

(p niinkuin "particular")

Tod: Perustuu pelkästään lineaarisuuteen. (Aivan sama lause pätee vaikkapa tavalliselle lineaariselle yhtälösysteemile.)

1) Olkoon H (EHY):n ratkaisuavaruus ja

>   

y[h]*epsilon*H

>   

T(y[h]+y[p]) = T(y[h])+T(y[p])

Tämä = z, koska T(y[h]) = 0  ja T(y[p]) = z .

Siis mielivaltainen väitettyä muotoa oleva summa on EHY:n ratkaisu.

2) Olkoon nyt y mielivaltainen EHY:n ratkaisu.

Tällöin

>   

y-y[p]*epsilon*H

sillä

>   

T(y-y[p]) = 0

koska

>   

T(y-y[p]) = T(y)-T(y[p])

ja molemmat viimemainitut ovat = z.

Siispä

>   

y-y[p] = y[h]

missä y[h]  on  jokin HY:n ratkaisu.

Niinpä mielivaltaisesti annettu EHY:n ratkaisu y  voitiin kirjoittaa vaaditussa muodossa.

Lineaarinen differyhtälösysteemi

Aivan samoin, kuten näemme differentiaaliyhtälöiden kohdalla, myös differenssiyhtälöiden  teoria ja käytäntö kannattaa opetella

palauttamaan systeemien käsittelyyn. (Tämä on ns. "moderni" käsittelytapa.)

  • Jokainen n:nnen kertaluvun yhtälö palautuu 1. kertaluvun nxn-systeemiin.
  • Systeemien teoria on kiintoisaa sinänsä (kaikki systeemit eivät ole peräisin jostain n:nnen kertaluvun yhtälöstä).

Lineaarisen homogeenisen systeemin yleinen muoto:

>   

x[k+1] = A*x[k]

missä A on nxn-matriisi ( x[k]   on   R^n :n  vektori.)

Esim. 7. Muuntaminen systeemiksi.

>    y[k+3]-2*y[k+2]-5*y[k+1]+6*y[k] = 0

y[k+3]-2*y[k+2]-5*y[k+1]+6*y[k] = 0

Merkitään

>   

x[k] = Vector(%id = 12583436)

>   

x[k+1] = Vector(%id = 12758760)

y[k+3] lasketaan yhtälöstä:

>   

x[k+1] = Vector(%id = 12583596)

Siis

>    x[k+1] = A*x[k]

x[k+1] = A*x[k]

missä

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

A := Matrix(%id = 12847508)

Jos onnistumme ratkaisemaan tämän 1. kertaluvun systeemin, niin ratkaisuvektorin x[k]  1. koordinaatti antaa yhtälön ratkaisukoordinaatin y[k] .