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 )