-
mlODE001
mlD001.tex (iv3 harj.1 av 2001)
Tätä tehtävää harjoitellaan Matlab-teknisesti loppuviikon harjoituksissa. Osattava esittää (ke 19.9.) liitutaululla käsin piirtäen periaatteessa.
Tarkastellaan diffyhtälöä \(y'=y-x\) alueessa \(-1 \leq x,y \leq 1\),
Piirrä \(xy\)-koordinaatistoon suuntakenttä ja isokliinejä käsin ja kokeile myös Matlabia laskimen roolissa. Ota hilaväliksi aluksi vaikka \(h=0.5\).
Matlab-laskussa kannattaa muodostaa matriisi, sanokaamme K, jonka alkioina ovat arvot \(y_i - x_j\), \(i=1\ldots 5, j=1\ldots 5\) tähän tapaan:
h=0.5; t=-1:h:1;x=0:h:2
for i=1:5
for j=1:5
K(i,j)=y(i)-x(j)
end
end
Koska matriisissa rivi-indeksi i juoksee alaspäin, on helpompaa sijoittaa arvot koordinaatistoon kääntämällä matriisin sarakkeet ylösalaisin; tämä tapahtuu komennolla flipud. Katso siis matriisista flipud(K) arvot ja merkitse ne kynällä piirrokseen. Tarkista käsin (tai “skalaarilaskimella” (Matlabkin käy)) muutama alkio ainakin. (Hilan tihentäminen käy nyt helposti muuttamalla vain yllä \(h\):ta.)
Myöhemmin opimme, että K-matriisin muodostaminen käy kätevämmin, tehokkaammin ja rutiininomaisemmin näin:
x=...;y=... % Kuten edellä
[X,Y]=meshgrid(x,y);
K=Y-X; % Jos esiintyy kertolaskua, potenssia ym.,
% on varustettava pisteellä, siis
% .*, .^ ./ (taulukko-operaatiot)
Nyt meillä on kaikki data kerättynä suuntakentän piirtämistä varten. Voit katsoa skriptiä suuntak1.m alla aputiedostossa ja miettiä, mitä siinä tapahtuu. (myös: http://www.math.hut.fi/teaching/v/matlab/suuntak1.m)
Lopuksi voit kokeilla LAODE-funktiota dfield8 (versio v. 2012) , jonka käyttö ei vaadi Matlabin tuntemista lainkaan. http://math.rice.edu/dfield/
Vaativuus 2.
Tehtävän Latex-koodi:
../matlabteht/mlODE/mlODE001.tex
Avainsanat:MatlabODE, diffyhtalot, mlDifferentiaaliyhtälöt,suuntakenttä
Viitteet: [LAODE] Golubitzky-Dellnitz: Linear Algebra and Differential Equations using Matlab, Brooks/Cole 1999.
-
mlODE002
mlODE002.tex [vastaava Maple: mplteht/mplODE003.tex] (iv3/2001, harj. 1)
Millä xy-tason käyrillä on ominaisuus: Käyrän tangentin kulmakerroin jokaisessa pisteessä \((x,y)\) on \(-\frac{4 x}{y}\) ?
Ratkaise yhtälö muuttujien erottelulla (“separation of variables”). Piirrä suuntakenttä isokliineja apuna käyttäen käsin vaikkapa alueessa \([-2,2] \times [-2,2]\).
Kokeile myös LAODE-funktiota dfield8. Tässä on käytettävä ahkerasti stop-näppäintä, ratkaisu ajautuu aina ongelma-alueelle, mikäli x-akseli on mukana.
Voit myös täydentää kuvaa alussa laskemillasi ratkaisukäyrillä tyyliin:
x=linspace(-2,2,30);y=x; [X,Y]=meshgrid(x,y);
Z=... % muista pisteitt\"aiset laskutoimitukset.
contour(x,y,Z,1:10); shg
Kokeile ja selitä!
Hae m-tiedosto dfield8 sivulta http://math.rice.edu/~dfield/ ja sijoita se Matlab-polkusi varrelle.
Kirjoita Matlab-istuntoon : dfield8
Vaativuus: 2
Tehtävän Latex-koodi:
../mlteht/mlODE/mlODE002.tex
Viitteet:
[LAODE] Golubitzky-Dellnitz: Linear Algebra and Differential Equations using Matlab, Brooks/Cole 1999.
Avainsanat: mlODE,Matlabdiffyht,differentiaaliyhtaloita Matlab:lla,suuntakenttä, isokliinit
-
mlODE003
mlD003.tex (iv3/2001, harj. 1, teht. 3) (Jatkoa: mlD007.tex)
Lisääntymiskykyinen populaatio, jonka lisääntymistä rajoittavia tekijöitä ei ole, noudattaa yleensä likimain Malthus’n lakia, ts. nykyhetkellä kasvunopeus on verrannollinen populaation nykykokoon.
Muodosta ilmiölle differentiaaliyhtälömalli ja ratkaise.
Ratkaisussa esiintyy vakiot \(y_0\) = populaation koko alkuhetkellä ja \(k\) = "lisääntymiskykyvakio".
Täydennämme vanhaa tuttua USA:n väkilukutaulukkoa arvoilla, jossa toinen rivi ilmaisee väkiluvun miljoonissa ja ensimmäinen vuosiluvun.
1800 1830 1860 1890 1920 1950 1980
5.3 13 31 63 106 150 230
a) Määritä \(k\) kahden ensimmäisen datasarakkeen avulla. Tutki, mihin vuosilukuun saakka malli antaa siedettävän tuloksen ja mistä alkaen sietämättömän.
b) Määritä \(k\) ensimmäisen ja viimeisen datasarakkkeen perusteella (jolloin kyseessä ei enää ole tulevaisuuden ennustaminen), ja tutki tämän mallin käyttäytymistä muissa datapisteissä.
c) Havainnollista kuvilla, käytä tarpeen mukaan logaritmista skaalaa.
Voit valita ajan \(0\)-hetken vuosiluvuksi \(1800\) (miksi?). Toki ei mitenkään välttämätöntä.
Avainsanat:MatlabDy, diffyhtälöt, mlDifferentiaali(yhtälöt), populaationkasvumalli, Malthus’n laki
Vastaavanlaisia tehtäviä: http://matriisi.ee.tut.fi/~piche/ode/ Pichet: TTY 1999, course 73107 (hyvä kokonaisuus):
1) Perusesim tähän kohtaan: In 1980 the population of Altlantia was 254512; in 1995 it was 294726. What is its population in 2000. (Vuosilukuihin voisi lisätä nyt vaikka 15.)
2) Jatkoa: mlD007.tex
-
mlODE004
mlODE004.tex [Maple-versio: ...mplODE004.tex] (iv3/2001, harj. 1, LV teht. 1-2)
Laskuvarjohyppääjän yhtälö. Oletetaan, että hyppääjän + varustuksen massa = \(m\) ja ilmanvastus on verrannollinen nopeuden neliöön, olkoon verrannollisuskerroin \(=b\). Tällöin Newtonin 2. laki antaa liikeyhtälön:
\[mv'=mg - b v^2 .\]
Olkoon yksinkertaisuuden vuoksi \(m=1, b=1\) ja \(g=9.81 m/s^2\).
Piirrä suuntakenttä.
Oletetaan, että laskuvarjo aukeaa, kun \(v=10 m/s\), valitaan tämä alkuhetkeksi \(t=0\). Piirrä tämä ratkaisukäyrä suuntakenttäpiirrokseen. Yritä nähdä suuntakentästä, että kaikki ratkaisut näyttävät lähestyvän rajanopeutta \(v\approx 3.13\) ja että ratkaisut ovat joko kasvavia tai pieneneviä (ja millä alkuarvoilla mitäkin, ja mitä tarkoittaa fysikaalisesti)
Määritä rajanopeus suoraan yhtälöstä.
Käytä Matlab-piirroksiin funktiota dfield8 ja Maplessa DEtools-kirjaston DEplot-funktiota.
Vihje Kts. [HAM] ss. 169-170 tai ?DEplot
> with(DEtools)
> with(plots)
Suuntakenttään:DEplot,
grafiikkojen yhdistämiseen: display.
dfield-ohje:
Hae m-tiedosto dfield8 sivulta http://math.rice.edu/\(\sim\)dfield/ ja sijoita se Matlab-polkusi varrelle.
Kirjoita Matlab-istuntoon : dfield8
Avainsanat: MatlabDy, MapleDy, diffyhtälöt, suuntakenttä, isokliinit, mplDifferentiaali(yhtälöt), mlDifferentiaali(yhtälöt)
Viitteet: [HAM] Heikki Apiola: Symbolista ja numeerista matematiikkaa Maple-ohjelmalla, Otatieto 588, 1998
-
mlODE005
mlODE005.tex (Pichet Course 73107)
Määritä diffyhtälön \[x'=x^2-t\, x + 4t\] suuntakenttä ja pisteen \((x,t)=(-1,-1)\) kautta kulkeva ratkaisukäyrä graafisesti funktion dfieldxx avulla, missä \(xx=8\) (v. 2012).
Mikä on tämän funktion minimipisten likiarvo (vaikka 2:n tai 3:n numeron tarkkuudella).
Vihje:
Hae m-tiedosto dfield8 sivulta http://math.rice.edu/\(\sim\)dfield/ ja sijoita se Matlab-polkusi varrelle.
Kirjoita Matlab-istuntoon : dfield8
Kun sinulla on yhtälö, suuntakenttä ja ratkaisukäyrä, voit valita Edit/Zoom in, jolla pääset tarkentamaan minipisteen hakua (askeleen kerrallaan).
Ennen minimin hakua kannattaa poistaa piirtämäsi ratkaisukäyrät Options/Erase solutions, niin Options/Keyboard input-valinnalla voit antaa tarkan alkuarvon.
Vaativuus: 1+
Tehtävän Latex-koodi:
../mlteht/mlODE/mlODE005.tex
Avainsanat: mlODE,Matlabdiffyht,differentiaaliyhtaloita Matlab:lla,suuntakenttä, dfield8
Matlabfunktioita: ode45,dfield8(Rice University)
-
mlODE006
mlD006.tex
Kun Marsiin laskeutuvan luotain laukaisee laskeutumisrakettinsa, sen nopeutta kuvaa differentiaaliyhtälö \[\frac{dV}{dt} = g_m - \frac{k}{m}V^2,\] missä \(g_m = 3.688\) on Marsin painovoimakerroin, \(k= 1.2\) on aerodynaaminen vastusvoima, ja \(m=150\) on luotaimen massa kilogrammoissa. Käytä vapaan pudotuksen rajanopeutta alkuarvona (Marsissa rajanopeus on 67.056 m/s), ja ratkaise yhtälö välillä \([0:0.05:6]\). Piirrä sekä nopeuden että kiihtyvyyden kuvaajat.
-
mlODE007
mplODE017.tex, mlODE007.tex
Käyrän sovituksen yhteydessä huomasimme, että eksponentiaalinen kasvumalli, ns. Malthus’n laki \(y'=k y\) ei toimi USA:n väestödataan pitkällä aikavälillä. Mallia voidaan tarkentaa lisäämällä sopiva kasvua rajoittava termi, tällöin johdutaan ns. logistiseen kasvulakiin:
\(y'=a y - b y^2\)
USA:n väestödataan liityen Verhulst arvioi v. \(1845\) arvot \(a=0.03\) ja \(b=1.6 10^{-4}\), kun \(t\) mitataan vuosissa ja väkiluku \(y(t)\) miljoonissa.
Opettajalle: Tehtävä voidaan käsitellä ehkä luontevamminkin kokonaan erillisenä numeeristen diffyhtälöratkaisujen opetuksesta. Tällöin otetaan vain alla olevat kohdat (c) ja/tai (d).
(a) Ratkaise tehtävä (\(y(0)=5.3\)) Eulerin menetelmällä käyttämällä askelpituussa \(h=10\)
(b) rk4:llä käyttäen n. nelinkertaista askelta (voit kokeilla pienempiäkin)
(c) Matlabin ode45:llä.
(d) Laske analyyttinen ratkaisu Maplella. (Kyseessähän on Bernoullin yhtälö.)
Piirrä kuvia ja laske kaikissa tapauksessa ratkaisujen arvot annetuissa taulukkopisteissä. (ode45-tapauksessa onnistuu ainakin sovittamalla dataan splini funktiolla spline, joka on maailman helppokäyttöisin.)
kts. http://www.math.hut.fi/teaching/v/matlab/opas.html#splinit
(Nykyään (2012) ei tarvita erillistä splinisovitusta, laskentapisteet voidaan antaa suoraan ode45-funktiolle syötteenä.)
function [T,Y]=eulerS(f,Tspan,ya,n)
% Tämä vain kehittely- ja opettelutarkoituksessa.
% Funktio eulerV hoitaa niin skalaari- kuin vektoriversion.
% (24.2.04, modifioitu 21.8.2010)
% Esim: y'=t+y, y(0)=1
% f=@(t,y)t+y
% [T,Y]=eulerS(f,[0 4],1,6), plot(T,Y,T,Y,'.r');shg
a=Tspan(1);b=Tspan(2);
h=(b-a)/n;
Y=zeros(n+1,1);T=(a:h:b)'; %Pystyvektorit yhdenmukaisesti ode45:n
Y(1)=ya; % kanssa
for j=1:n
Y(j+1)=Y(j)+h*f(T(j),Y(j));
end;
Vaativuus: 2+
Tehtävän Latex-koodi:
../mplteht/mplODE/mplODE017.tex
Ratkaisu: *** Hukassa, etsi/tee ***
Viitteitä:
http://math.aalto.fi/opetus/kp3-ii/06/L/L14dynumkalvot.pdf
http://www.math.hut.fi/ apiola/matlab/opas/lyhyt/esim/eulerS.m
(Listaus yllä)
Avainsanat: mplODE,Maplediffyht,differentiaaliyhtaloita Maplella, logistinen kasvu, numeeriset menetelmät, Eulerin menetelmä, Heunin menetelmä, Runge Kutta
Maplefunktioita: dsolve
-
mlODE008
mlD008.tex,mplD018.tex
Tarkastellaan yhtälöä \(y'=-2\alpha(t-1)y\). Ratkaise aluksi analyyttisesti (saat käyttää Mapleakin.)
Totea kuvasta ja derivaattaehdosta yhtälön stabiilisuus/epästabiilisuusalueet. Ota kuvassa ja aina tarvittaessa vaikkapa \(\alpha =5\).
Ratkaise yhtälö sekä Eulerilla että BE:llä. Sopivia arvoja voisivat olla vaikkapa \(h=0.2\), väli: \([1,4.5]\), \(y(1)=1\).
Vertaa kokeellisesti stabiilisuukäyttäytymistä teorian ennustamaan ja pane merkille, miten epästabiilisuus käytännössä ilmenee.
Tämä tehtävä soveltuu erityisen hyvin Maplella tehtäväksi, se on pitkälle ideoitu [HAM] sivulla 124, myös Euler ja BE ovat valmiina. (Koodit saa kurssin maple-hakemistosta.) ** Tulee aputiedostoon **
** apu puuttuu, editoi viitteet! **
Viitteitä:
http://math.aalto.fi/opetus/kp3-ii/06/L/L14dynumkalvot.pdf
http://www.math.hut.fi/~apiola/matlab/opas/lyhyt/esim/eulerS.m
(Listaus yllä)
-
mlODE016
Maple,Matlab (H2T10)
Ratkaise alkuarvotehtävä \[y' -y = \cos x,\ \ y(0)=1\]
analyyttisesti Maplella ja numeerisesti Matlabilla. Piirrä ratkaisukäyrä.
Anna alkuarvoksi symboli \(c\) ja piirrä ratkaisukäyräparvi sopivalla välillä, kun \(c=-0.9,-0.8,\ldots ,0.\)
Miltä parvi näyttää suurilla \(x:\)n arvoilla. Tässä pitäisi erottua kolmenlaista käytöstä.
Maple: dsolve, Matlab: ode45
Avainsanat: Differentiaaliyhtälö, alkuarvotehtävä, analyyttinen ratkaisu, numeerinen ratkaisu.
-
mlODE017
mplD007.tex (iv3/2001, harj. 2, AV teht. 1)
Ratkaise (AA)-tehtävä \(y' - 2 x y = 1,\ \ y(0)=-0.5\)
Tässä näyttää siltä, että (EHY):n erikoinen olisi helppo löytää, mutta huomaat pian, että luonnolliset yritteet eivät toimi. (Kyseessähän on lineaarinen, mutta ei-vakiokertoiminen yhtälö.)
Ratkaise vaan sitten kiltisti integroivan tekijän menettelyllä.
Integrointi johtaa erf-funktioon, Maple antaa sen suoraan, voit myös konsultoida KRE-kirjaa hakusanalla erf. Lausu siis ratkaisu erf:n avulla.
Piirrä suuntakenttäpiirros Maplen DEtools-pakkauksen DEplot-funktion avulla (kts [HAM] s. 169), voit toki käyttää myös Matlab:n dfield8-funktiota (ohje alla).
Valitse alkuarvoja \(y_0\) väliltä \((-1,-0.5)\) yrittäen löytää kriittistä arvoa \(y_0\), joka jakaa ratrkaisukäyrät plus tai miinus ääretöntä lähestyviin. (Tuo kriittinen ratkaisukäyrä on rajoitettu.) Käytä hyväksesi erf-funktion ominaisuutta \(\lim_{x\to \infty} erf(x) = 1\) laskeaksesi tarkan arvon \(y_0\):lle.
dfield-ohje: Hae m-tiedosto dfield8 sivulta http://math.rice.edu/~dfield/ ja sijoita se Matlab-polkusi varrelle.
Kirjoita Matlab-istuntoon : dfield8
Avainsanat: MapleDy, diffyhtälöt,erf, mplDifferentiaali(yhtälöt)
Viitteet: [KRE] E. Kreyszig: Advanced Engineering Mathematics, Wiley
[HAM] Heikki Apiola: Symbolista ja numeerista matematiikkaa Maple-ohjelmalla, Otatieto 588, 1998.
-
mlODE020
Piirrä tasokäyrä \((x,y)=(e^{-t/20}\cos\, t, e^{-t/10 \sin\, t})\) (faasitaso).
Piirrä edelliset koordinaattikäyrät \(x(t)\) ja \(y(t)\) (aikataso).
Muista (.*), kun muodostat vektorilausekkeita. kokeile eri parametrivälejä. Kokeile linspace(a,b,n):ssä muitakin kuin oletusarvoa \(n=100.\) Käytä axis-komentoa. Ja kysele help:llä.
Uusi grafiikkaikkuna: >> figure, Useampia samaan kuvaan: >> hold on
Avainsanat: Differentiaaliyhtälöryhmä, aikakuva,faasikuva,grafiikka,plot
-
mlODE100
mlODE100.tex
Mallinnetaan heiluria kolmessa ulottuvuudessa. Pisteen \(\mathbf{x}(t) = [x_1(t),x_2(t),x_3(t)]^T\), massa on \(m\), ja siihen kohdistuu voima \(\mathbf{F} = [0,0,-gm]^T\) (\(g = 9.81m/s^2\)). Pistettä rajoittaa pallon pinnan yhtälö \(\Phi(x)=x_1^2+x_2^2+x_3^2-1=0\). Pisteen liikettä mallintaa seuraava differentiaaliyhtälö: \[\mathbf{x}''(t) = \frac{1}{m}\bigg( \mathbf{F}- \frac{m\mathbf{x}'(t)^T\mathbf{H}\mathbf{x}'(t)+ \nabla \Phi\mathbf{F} }{|\nabla \Phi|^2} \nabla \Phi \bigg)\]
Yhtälössä \(\nabla \Phi\) on funktion \(\Phi\) tilagradientti (\(x_i\) komponenttien suhteen) ja \(\mathbf{H}\) on Hessin matriisi funktiolle \(\Phi\) – tässä tapauksessa diagonaalimatriisi, jonka alkiot ovat kaikki \(2\).
Lisäksi yhtälö tarvitsee alkuarvot \(\mathbf{x}(0)= \mathbf{x}_0\) ja \(\mathbf{x}'(0) = \mathbf{v}_0\).
Ratkaise yhtälö kertalukua tiputtamalla, ja käyttämällä MATLAB-ratkaisinta ode45. Tarkastele tulosta graafisesti: ovatko trajektorit sitä mitä odotit? Ratkaise ongelma sitten käyttämällä metodia ode23. Vastaavatko tulokset nyt sitä mitä odotit?
Vaativuus: 3+
Tehtävän Latex-koodi:
../mlteht/mlODE/mlODE100.tex
Avainsanat: mlODE,Matlabdiffyht,differentiaaliyhtaloita Matlab:lla, 3d-heiluri
Matlabfunktioita: ode45,ode23, odexx
-
mlODE102
mlODE102
Ratkaise toisen kertaluvun alkuarvotehtävä \[y''-0.05y'+0.15y = 2t; y'(0)=0,y(0)=0.\] Matlabin ratkaisijoita (ODE45, ODExx) varten 2. kertaluvun yhtälö tulee muuntaa 1. kertaluvun systeemiksi määrittelemällä \(y_1 = y, y_2= y'\), jolloin saadaan systeemi: \[\begin{cases}
y_1' = y_2 \\
y_2' = ...
\end{cases}\]
Piirrä ratkaisukäyrä ja derivaatta sekä faasitasokuva.
Vaativuus: 1+
Tehtävän Latex-koodi:
../mlteht/mlODE/mlODE102.tex
Avainsanat: mlODE,Matlabdiffyht,differentiaaliyhtaloita Matlab:lla, numeeriset menetelmät, faasitaso
Matlabfunktioita: ode45,odexx
-
mlODE103
Eulerin menetelmä on klassinen menetelmä differentiaaliyhtälöiden ratkaisemiseksi. Menetelmän lähtökohtana on yhtälö \[y'(t) = f(t,y(t)); y(t_0)= y_0,\] josta koetetaan ratkaista \(y(t)\) kun \(t= [t_0, t_n]\). Menetelmä perustuu \(t\):n diskretoinnille: \(t\) muunnetaan joksikin vektoriksi \(T = (t_0, t_0+h, \ldots, t_n)\), jonka jälkeen koetetaa etsiä vastaavia arvoja \(y(t_0+ih)\). Merkintöjen selvyyden vuoksi kirjoitetaan \(t_0+jh = t_j\). Eulerin menetelmässä approksimoidaan arvoa \(y(t_i)\) seuraavasti: \[y(t_i) = y(t_{i-1})+ h\, f(t_{i-1},y(t_{i-1})).\]
Ratkaise Eulerin menetelmällä differentiaaliyhtälö \[y'(t) = y(t) \sin(t).\] Käytä eri \(h\):n arvoja, ja vertaa oikeaan ratkaisuun \(e^{1-\cos(t)}\).
Eulerin menetelmä MATLABissa on helpointa toteuttaa yksinkertaisena silmukkana. Olettaen, että yhtälön oikea puoli on kirjoitettu funktioon f(t,y), Euler-iteraatio voidaan tehdä seuraavasti:
for n = 2:length(t)
y(n) = y(n-1)+h*f(t(n-1),y(n-1));
end
missä \(t\) on vektori joka on ratkaisuvälin diskretointi, ja \(h\) on \(t\):n pisteiden välinen etäisyys.
-
mlODE104
Rungen-Kutan menetelmä differentiaaliyhtälöiden ratkaisemiseksi on huomattavasti kehittynyt versio Eulerin menetelmästä. Menetelmässä määritellään seuraavasti: \[\begin{array}{l }
k_1 = f(x_n, y_n),\\
k_2= f(x_n+\frac{h}{2},y_n+\frac{1}{2}hk_1),\\
k_3= f(x_n+\frac{h}{2},y_n+\frac{1}{2}hk_2),\\
k_4= f(x_n+h,y_n+hk_3).\\
\end{array}\] Diskretointiaskeleeksi tulee nyt \[y_{n+1} = y_n+hk = y_n+ \frac{h}{6}(k_1 + 2k_2+2k_3+k_4).\]
Ratkaise edellisen tehtävän yhtälö \[y'(t) = y(t) \sin(t); y(0)=1.\] käyttämällä Rungen-Kutan menetelmää. Vertaa tulosta sekä todelliseen ratkaisuun \(y(t)=e^{1-cos(t)} \), että Eulerin menetelmällä saavutettuun. Mitä huomaat tarvittavasta askelkoosta?
Ratkaise alkuarvotehtävä \[y'(t) = \sqrt{1-y^2}; y(0) = 0\] käyttämällä Rungen-Kutan menetelmää.
-
mlODE105
Ratkaise differentiaaliyhtälöryhmä \[\begin{cases}
\frac{dx}{dt} = -\sigma x+ \rho y \\
\frac{dy}{dt} = \sigma x -y-xz \\
\frac{dz}{dt} = -\beta z + xy \\
\end{cases}\] numeerisesti välillä \([0,20]\), kun \(\sigma = 10, \rho = 28 \text{ ja } \beta = 8/3\). Piirrä ratkaisukäyrät samaan kuvaan, ja piirrä käyrät \(x(t)\) ja \(z(t)\) parametrisesti. Tämän jälkeen piirrä 3-ulotteinen parametrisoitu käyrä kaikista koordinaateista.
Onko ratkaisu rajoitettu? Suppeneeko se kohti jotain arvoa?
Kokeile muuttaa alkuarvoja, sekä parametrien arvoja. Vallitsevan teorian mukaan systeemi on kaoottinen dynaaminen systeemi, jonka käyttäytyminen voi muuttua merkitsevästi jo pienistä muutoksista lähtötilanteessa; itse asiassa termi perhosvaikutus keksittiin kuvaamaan juuri tämän systeemin käytöstä.
Kolmiulotteinen parametrisoitu käyrä (tai pistejoukko) piirretään MATLABissa funktiolla plot3.
-
mlODE106
Kirjoita heiluriyhtälö \(\Theta'' + \frac{g}{L} \sin(\Theta) = 0\) ensimmäisen kertaluvun systeemiksi ja samantien Matlab-funktioksi (joko inline tai m-tiedosto). Voit ottaa \(g/L = 1.\)
Laske ratkaisu sopivalla aikavälillä (esim. \(\left[0,10\right]\)) ja kolmella erilaisella alkuarvolla, joilla saat erityyppiset ratkaisut. Käytä ode45-funktiota.
Piirrä ratkaisukäyrät aikatasoon ja trajektorit faasitasoon.
-
mlODE107
Reuna-arvotehtävä tähtäysmenetelmällä
(Moler odes teht. \(7.19\) ss. 45–47.) Olkoon ratkaistavana reuna-arvotehtävä \(y''=y^2-1,\ \ y(0)=0,\ \ y(1)=1.\) Tee Matlab-funktio, joka ottaa argumentikseen "alkunopeuden" \(0\):ssa ja palauttaa vastaavan AA-tehtävän ratkaisun arvon \(1\):ssä, vähennä siitä vielä 1, niin sinulla on funktio, jota sopii tarjota fzero:lle.
-
mlODE109
Ratkaise differentiaaliyhtälö \[\frac{\partial u(t,x)}{\partial t} = -u-5e^{-t}\sin(5t)\] nelikulmiossa \([0,5]\times[-2,2]\), alkuarvokäyrällä \(u(0,x)= e^{-x^2}\).
Piirrä tuloksista kuva:
Määrittele ensin sopivan harva diskretaatio välille \([-2,2]\), ja sen jälkeen tähän diskretaatioon sopiva alkuarvokäyrä, sen jälkeen ratkaise differentiaaliyhtälö numeerisesti (ode45) vuoroittain kullekin alkuarvolle.
-
mlODE110
Ratkaise nk. Rösslerin systeemi \[\begin{cases}
y_1'(t) = -y_2(t)-y_3(t), \\
y_2'(t) = y_1(t)+ay_2(t), \\
y_3'(t) = b+y_3(t)(y_1(t)-c),
\end{cases}\] missä \(a,b,c\) ovat parametreja, välillä \(t \in [0,100]\). Käytä alkuarvoja \(y(0)=[1,1,1]^T\) ja \((a,b,c) = (0.2,0.2,2.5).\)
Kyseessä on kaoottinen systeemi, ts. ratkaisurata riippuu suuresti joko alkuarvoista ja parametreista. Saatuasi ratkaisun, piirrä ratkaisun kuvia alkuarvon muuttujana, ja tutki kuinka pienet muutokset aiheuttavat havaittavia eroja.
-
mlODE200
Ratkaise differentiaaliyhtälö \[y'(t) = -\frac{1}{t^2}+10 \bigg( y-\frac{1}{t} \bigg); y(1) = 1.\] käyttäen Backward-Euler menetelmää: \[y_{n+1} = y_{n} + hf(t_{n+1,y_{n+1}}).\] Jotta tämä onnistuisi, tulee jokaisella iteraation askeleella ratkaista yhtälö \[y_{n+1} - hf(t_{n+1,y_{n+1}}) = y_n\] \(y_{n+1}\):n suhteen. Tämä tarkoittaa, että yleisen ratkaisimen kirjoittaminen olisi vaikeaa, mutta tässä eksplisiittisessä tapauksessa se onnistuu.
Vertaa sitten saamaasi tulosta MATLABin ode45 funktiolla saamaasi: kuinka selität eron?
