Matlab-tehtäviä 4

Ratkaisuja tehtäviin voit kurkistella täältä. Kannattaa toki miettiä ensin itse, toisaalta ei ole syytä jokaiseen tehtävään käyttää liikaa aikaa. Hae itsellesi sopiva etenemistyyli, kunhan otat opiksesi.

Käyrän sovitus, epälineaariset yhtälöt

Teht. 1.
Eräs kemiallinen koe tuotti seuraavat datapisteet:
tdata=[-1 -0.960 -0.86 -0.79 0.22 0.5 0.93];
ydata= [ -1 -0.151 0.894 0.986 0.895 0.5 -0.306];
Tarkoitus on estimoida y(t):n arvoja välillä [-1,1]
  1. Piirrä datapisteet.
  2. Muodosta interpolaatiopolynomi ja piirrä samaan kuvaan.
  3. Muodosta asteita 2,3,4 olevat PNS-polynomit ja piirrä samaan kuvaan
  4. Sovita vielä splini ja piirrä samaan.
  5. Asettele grafiikkaikkuna paremmaksi vaikka tyyliin:
      a=min(xdata)-0.1;b=max(xdata)+0.1;
      c=min(ydata)-0.1;d=max(ydata)+0.1;
      axis([a b c d])
    
    Käytä myös legend-komentoa ja kokeile grid on
Teht. 2.
Edellä kirjoittamasi omanewton soveltuu luonnollisesti polynomiyhtälöön, mutta vaatii yleensä hyvän alkuarvauksen. Historiallisesti mielenkiintoinen yhtälö on
             x3 - 2x - 5 = 0,
jota Wallis-niminen matemaatikko käsitteli, kun hän ensi kertaa esitteli Newtonin menetelmää Ranskan akatemialle. [Lähde: Moler NCM]
  1. Piirrä kuvaaja saadaksesi alkuarvon Newtonin menetelmälle reaalijuurta varten.
  2. Määritä reaalijuuri omanewton:lla
    Neuvo: Polynomifunktion voit määritellä tähän tapaan, kun p on kerroinvektori:
                  pf=@(x) polyval(p,x)
    
    Derivaatan saat polyder- (tai omapolyder)-funktiolla.
  3. Yritä löytää kompleksijuuri antamalla kompleksisia alkuarvoja. (Jos löydät yhden, niin toinen on sen liittoluku.)
  4. Määritä juuret roots-funktion avulla.

  
Teht. 3. Bisektio
Tutustu [CvL]-skriptiin ShowBisect (Löydät sen myös täältä). Kopioi se omaan työhakemistoosi (vaikka avaamalla Matlab-esitori ja cut/paste.)
Suorita joitakin esimerkkiajoja, ilmeinen eka kokeilu voisi olla 'sin' ja pi/2, 3*pi/2. (Kirja on senverran vanha, että funkio-osoittimia ei vielä ollut. Siitä johtuen koodissa on feval-komentoja arvon laskentaan. Nykyisin niitä ei tarvita, kun käytetään @-tyyliä.)

Tässä on bisektio-algoritmin hahmotelma. Täydennä se funktioksi (vastaavaan tapaan kuin Newton edellä). Ota tol valinnaiseksi parametriksi, jonka oletusarvona on eps. Tässä testataan nargin-arvoa, jos se on 2, asetetaan tol=eps.

maxiter=100;
k=0
 while(abs(b-a)) > tol*abs(b)
     if (k>maxiter) 
         return
     end;    
     x=(a+b)/2;
     if sign(f(x)) == sign(f(b))
       b=x;
     else
       a=x;
     end
     %disp([k a b])
     k=k+1;
 end
 

Teht. 4.
Matlabin fzero käyttää yhdistelmää bisektiosta, "sekanttimenetelmästä" (karkeasti ottaen Newton, jossa derivaatallem käytetään differenssiapproksimaatiota), sekä käänteistä interpolaatiota.

Tutstu ja kokeile oppaan fzero-esimerkkiin Kopioi se m-tiedostoksi ja suorita.
Suorita help fzero . Määritä 10 pienintä positiivista ratkaisua yhtälölle tan(x) = x

Teht. 5.
Keplerin mallin mukaisessa planeettaradassa esiintyy luku E, "eksentrinen anomalia", joka toteuttaa yhtälön
           M = E - e sin(E)
missä M on keskianomalia ja e radan eksentrisyys. Tässä tehtävässä voit käyttää arvoja: M=24.851090, e=0.1.
Tee m-tiedosto (skripti,pääohjelma), jossa määrittelet kohdefunktion funktio-osoittimena tyyliin:
       M=.... ; e=....;
       f=@(E) .....
Tämän f:n annat (funktio)argumenttina funktiolle fzero. (Voit jopa antaa sen "anonyymina funktiona" suoraan fzero:n kutsussa.)
Kts. help fzero, jossa tehdään (vain) hiukan mutkikkaammin. Molerin kirjassa käytetään tällä kohdalla vanhentunutta tapaa(!). Paras esitys on [Hig-Hig]-kirjassa (siinä "parhaassa").
Teht. 6.
Sopiva lisätehtävä: [Moler NCM] Ch4 s. 138 teht. 4.16 (vesijohdon jäätyminen)

ODE - differentiaaliyhtälö(systeemit)

help ode45 antaa tietoa kaikista muistakin Matlabin ODE-ratkaisijoista. Perussyntaksi on kaikilla sama.

Esimerkki: yksi yhtälö

              y'=f(t,y)=-y-5e-tsin(5t), y(0)=y0=1 
Ratkaistaan numeerisesti aikavälillä 0 < t < 3 Matlab:n ratkaisijalla ode45.
Diffyhtälön määrittelevä funktio f(t,y) voitaisiin tässä tapauksessa kirjoittaa suoraan funkito-osoittimena (function handle), mutta tehdään ensin perustyylillä, joka sopii myös mutkikkaampiin yhtälöryhmiin.
  1. Aseta polku työhakemistoosi, avaa editori ja kirjoita tiedostoon omafun.m
      function yp=omafun(t,y)
      % omafun.m   ode-esimerkkifunktio
      yp=-y-5*exp(-t).*sin(5*t);
    
    Huom: Tässä toimisi myös "pisteetön" kertolasku, sillä Matlab suorittaa laskennan skalaaritoimituksin. Samahan pätee mm. numeerisessa integroinnissa. Pisteestä ei toisaalta ole haittaa, tämä huomautus on siksi, ettet ihmettelisi, miksi monissa dokuissa piste on jätetty pois.

    Testataan ensin, että funktiomme on oikein kirjoitettu:

    >> omafun(0,1)
    ans =
        -1
    >> omafun(0,2)
    ans =
        -2
    
    Syntaksi on kunnossa ja ainakin nämä lasketaan oikein.
  2. Komentoikkunassa ryhdytään ratkaisupuuhiin:
    >> tvali=[0 3];
    >> [T,Y]=ode45(@omafun,tvali,y0);
    >> plot(T,Y,'*--')
    >> xlabel t, ylabel y(t)
    
    Kuvassa * osoittaa adaptiivisen ratkaisijan askelpituuksien valinnat. Jos halutaan ratkaisun arvo muissa pisteissä, voitaisiin [T,Y]-dataan sovittaa splini. Mutta ode45 tekee senkin puolestamme. Jos halutaan arvot vaikka tasavälisessä t-pisteistössä, jossa on 40 pistettä, annetaan t-välin sijasta tuo pisteistö siis näin:
    >> tpisteet=linspace(0,3,40);
    >> [T,Y]=ode45(@omafun,tpisteet,y0);
    >> figure    % uusi kuvaikkuna
    >> plot(T,Y,'*--')
    
  3. Vaihtoehtoisesti voitaisiin välttää koko m-tiedostomäärittely tällaisessa yksinkertaisdessa tilanteessa ja hoitaa diffyhtälön määrittelykin "lennossa", jopa nimettömänä funktiona: Koko homma voitaisiin tiivistää tähän riviin:
    >> [T,Y]=ode45(@(t,y) -y-5*exp(-t).*sin(5*t) ,[0 3],1);
    
Teht. 7.
Suorita Matlabilla edellinen lasku, vaikkapa vain lyhintä kutsutapaa käyttäen. (Mutta sisäistä myös pitempi.)
Analyytinen ratkaisu saadaan tässä tapauksessa helposti (käsin tai symbolic toolboxilla tai jollain symbolilaskentaohjelmalla (Maple, Mathematica).)
                y(t) = e-tcos(5t) 

Määritä maksimivirhe välillä [0,3]. Kokeile eri pisteistöjä, jatka hienontamista niin kauan, että arvo ei muutu. Vast: (View source, mutta malta ...)

Suuntakenttä

Diffyhtälö määrää kussakin ty-tason pisteessä ratkaisukäyrän tangentin suunnan (y'=f(t,y)). Jos tasoalueen pistehilan pisteisiin asetetaan vastaava suuntanuoli, saadaan suuntakenttä. Matlab tarjoaa hyvät välineet myös tähän: Kokeilepa seuraavaa:
close all
n=16;
tpisteet=linspace(0,3,n); ypisteet=linspace(-1,1,n);
[t,y]=meshgrid(tpisteet,ypisteet); % Pistehila, muistele 3d-grafiikkaa.
pt=ones(size(y));
py=-y-5*exp(-t).*sin(5*t);
quiver(t,y,pt,py,1.5);
title('y''=-y-5e^{-t}sin(5t)','FontSize',14)
xlabel('t');ylabel('y','Rotation',0)
xlim([0 3.2]);ylim([-1.3 1.15]); % Viritellään akselirajoja.

Komento quiver(x,y,u,v,skaalaus) piirtää (meshgrid)hilapisteisiin [x,y] nuolet, joiden suuntavektorit ovat [u,v]. Nuolien pituus = skaalaus*[u(i),v(i)]-vektorin pituus.

Teht. 8.
Kirjoita edellä olevat komennot skriptiksi tiedostoon suuntak1.m ja lisää siihen diffyhtälön ratkaisu. Valitse joukko alkuarvoja ja piirrä samaan kuvaan suuntakenttä ja ao. alkuarvoja vastaavat ratkaisukäyrät. Voit täydentää skriptin vuorovaikutteiseksi käyttämällä alkuarvoille hiirisyöttöä funktiota ginput hyödyntäen.

Teht. 9. Esimerkki systeemistä, heiluri
Vaimentamaton heiluriyhtälö:
               Θ'' + sin(Θ)=0
Muutetaan 1. kertaluvun systeemiksi merkitsemällä: y1 = Θ, y2=Θ' , jolloin yhtälö saa muodon:
                   y1' = y2 
                   y2' = -sin(y1)
Systeemi koodataan ODE-ratkaisijalle muodossa:
  function yp = heiluri(t,y)
  % heiluri.m
  yp=[y(2);-sin(y(1))];
Siis kyseessä on funktio, jonka argumentteina on skalaari t ja vektori y (tässä tapauksessa y:n pituus on 2). Funktio palauttaa 2-pituisen derivaattavektorin yp. (yp:n on oltava pystyvektori.)
Huom: Tässä tapauksessa kysymyksessä on autonominen systeemi, ts. oikea puoli ei riipu t:stä. Siten t ei osallistu heiluri-funktion arvon laskentaan ja on siten "turha". Sen on oltava mukana siksi, että ode45 ja muut ode-funktiot toimivat niin autonomisille kuin yleisillekin, joten parametrilistan rakenteen on aina oltava sama.

  1. Muodosta heiluriyhtälön ratkaisut aikavälillä [0 10] ja alkuarvoilla y0a=[1;1], y0b=[-5;2], y0c=[5;2]. Piirrä sekä aikakuva, että erityisesti faasikuva.
  2. Piirrä faasikuvaan myös suuntakenttä Tätä ohjetta noudatellen.
  3. Kehittele vuorovaikuteinen hiirikäyttö ginput-funktion avulla, jos tuntuu mielekkäältä. Tarkempia ohjeita on saatavilla.

pplane7

Matlabin gui-ohjelmoinnilla tehty hieno väline (J. Pohlking).
pplane7 - Rice university, download
(pplane7 - Youtube )