\documentclass[finnish,11pt,twocolumn]{article}
\usepackage{/home/apiola/tex/v3harj}
%\include{/home/apiola/tex/defs}
%\usepackage{v23harj}
%\include{defs}

\begin{document}

%\vspace{0.2cm}
%\begin{verbatim}
%\end{verbatim}

\harjoitus{2}{39}{26 -- 28.9.2001}


\subsubsection*{Alkuviikko (AV)}

\begin{enumerate}

\item
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.

Ratkaise vaan sitten kiltisti integroivan tekijän menettelyllä.

Integrointi johtaa {\it erf}-funktioon, Maple antaa sen suoraan,
voit myös konsultoida KRE-kirjaa hakusanalla {\it erf}. Lausu siis
ratkaisu {\it erf}:n avulla.

Piirrä suuntakenttäpiirros vaikkapa {\tt dfield5}:llä (jos haluat,
voit käyttää myös Maplen {\tt DEtools}-pakkauksen {\tt DEplot}-funktiota.
(kts [HAM] s. 169)

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 {\it erf}-funktion
ominaisuutta $\lim_{x\to \infty} erf(x) = 1$ laskeaksesi tarkan arvon $y_0$:lle.


\item %2 
Tarkastellaan (AA)-tehtävää $xy'=4y,\ \ y(0)=1$. 

(a) Osoita, että tehtävällä ei ole ratkaisua. Osoita, että
tämä ei ole ristiriidassa $\exists_1$-lauseen kanssa. (Huom: Lauseen
avulla {\em ei voi todistaa epäeksistenssiä}, koska lauseen ehdot
eivä ole välttämättömät.)

(b) Vaihdetaan alkuehdoksi: $y(0)=0$. Miten nyt on ratkaisujen laita.

(c) Mitä voit sanoa alkuehdon $y(x_0)=y_0$ tapauksessa, jos $x_0 \neq 0$, \\
(1)  suoraan ratkaisukaavan avulla, (2) $\exists_1$-lauseen avulla.


\item
Muodosta {\it Picardin} iteraatiojonon muutama termi (AA)-tehtäville

(a) $y' = x+y,\ \ y(0)=0 $ \quad (b) $y' = x+y, \ \ y(0)=-1$ \\
(c) $y'=y^2,\ \ y(0)=1$.

Määritä myös tarkka ratkaisu.

LV-tehtävässä palataan asiaan Maple-hommana.


\item
Sovella {\it Picardin} iteraatiota (tuttuakin tutumpaan) (AA)-tehtävään

$y'=y,\ \ y(0)=1$. Osoita, että iteraatiojono lähestyy ratkaisufunktiota
$y(x)=e^x$ .

%$\forall x \ in \R$.
 

\item
Ratkaise alkuarvoteht\"av\"a $y'=t^2-y,\ \ y(0)=1$ Eulerin menetelm\"all\"a.
Vertaa ratkaisua tarkkaan ratkaisuun.
K\"aytt\"aytyyk\"o virhe odotetulla tavalla ($O(h)$)?
Laske arvot $y_1, ..., y_n$ seuraavasti:



K\"asipelill\"a, v\"ali $[0,0.2]$, 

a) $h=0.2,n=1$, \ \ b) $h=0.1, n=2$, \ \ c)  $h=0.05, n=4$

Matlabilla  V\"ali $[0,2]$, a) $h=0.2$, b) $h=0.1$,
               c) $h=0.05$.


\item
(Olkoon vaihteeksi $x(t)$ .)
Olkoon alkuarvotehtävänä $x'=x,\ \ x(0)=1$, jonka ratkaisu hyvin
tunnetaan: $x(t)=e^t$. 

Osoita, että jos lasketaan likiarvo $x_n=x_h(t_n)$ EM:llä pisteessä $t=t_n$ 
käyttäen askelpituutta $h$, niin $x_h(t_n)=c(h)^{t_n}$, missä
$c(h)=(1+h)^{1/h}$.

Osoita tämän nojalla, että kiinteällä $t=t_n$ pätee $\lim_{h\to 0}x_h(t)=e^t$.
(Vrt. Picard-tehtävään 4)

%Osoita lisäksi, että (globaali) katkaisuvirhe on muotoa $O(h)$.


%Sano myös, minkätyyppisiä diffyhtälöt ovat ominaisuuksiin
%{\em lineaarinen ja homogeeninen/epähomog., kertaluku, epälineaarinen}
%nähden.



\end{enumerate}


\subsubsection*{Loppuviikko (LV)}


\begin{enumerate}


\item
Muodosta {\it Picardin} iteraatiojonoa pitemmälle kuin AV-tehtävässä
samoille (AA)-tehtäville (a), (b), (c) ja lisäksi vielä (d): lle.

(a) $y' = x+y,\ \ y(0)=0 $ \quad (b) $y' = x+y, \ \ y(0)=-1$ \\
(c) $y'=y^2,\ \ y(0)=1$. \quad (d) $y'=3\frac{y}{x}$

Laske myös tarkka ratkaisu Maplella ja piirrä se ja iteraatiojonon funktioita.
(Jos tuntuu liian pitkältä, niin jätä yksi pois, hyvä olis saada kaikki
yhteisesti katetuksi (vaikka parityöskentelyssä sopimalla).

Malli: {\tt L/L4Picard1.html}, kts. myös [HAM] ss. 162--165 ({\tt dsolve})ja 
s. 126 {\it Picard--Lindelöf}


\item
Eulerin menetelmän etenemistä on helppo seurata rakentamalla
hyvin vähäeleinen skripti tähän tapaan. (Tässä on diffyhtälönä
siis $y'=y$)
\begin{verbatim}
f=inline('y','x','y')
x0=0;X=x0;y0=1;Y=y0;h=0.25;clf; hold on
x1=x0+h;y1=y0+h*f(x0,y0);X=[X,x1];Y=[Y,y1];x0=x1;y0=y1;
[X' Y'], plot(X,Y);shg
\end{verbatim}

Selvitä, mitä tässä tapahtuu ja kokeile. Iterointi sujuu 
Matlab-ikkunassa nuoli-ylös-näppäimellä toistamalla riviä 3.
Välillä voi aina käydä rivillä 4.
Piirrä myös tarkka ratkaisu samaan kuvaan.

Muuttele parametreja, myös diffyhtälöä.

\item
Tarkastellaan (AA)-tehtävää 
$$y'=\frac{3t^2}{(3y^2-4)},\ \ y(1)=0 .$$

(a) Laske EM:llä ratkaisuapproksimaatiot pisteissä $t=1.2,1.4,1.6,1.8$
käyttäen askelta $h=0.1$ .

(b) Tee sama askeleella $h=0.05$ .

(c) Vertaa tuloksia.

(d) Piirrä suuntakenttä ja ratkaisuapproksimaatioita, sekä
EM-ratkaisuja. Osaatko selittää, miksi EM toimii kohtuullisesti
alussa, mutta kelvottomasti lopussa?


\item
Tarkastellaan (AA)-tehtävää
$$y'={\frac{2\sqrt{y-\ln t} }{t}}+{\frac{1}{t}},\ \ y(1)=0$$
välillä $t \in [1,1.8]$
Ratkaise tehtävä 

\kohta{a}{ Eulerin menetelmällä askelpituudella $h=0.1$},

\kohta{b}{ Heunin menetelmällä  askelpituudella $h=0.2$},

\kohta{c}{ RK4- menetelmällä  askelpituudella $h=0.4$.}

Määritä tarkka ratkaisu Maple:n {\tt dsolve}-komennolla ja laske
sen avulla virheet, piirrä ja taulukoi kussakin tapauksessa.

Huomaa, että näillä askelpituuksien valinnoilla funktion arvojen 
laskentamäärät ovat samat.


\item
Huomasimme, että eksponentiaalinen kasvumalli, ns. {\it 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:

\quad $y'=a y - b y^2$

USA:n väestödataan liityen {\it Verhulst} arvioi v. $1845$ arvot
$a=0.03$ ja $b=1.6 10^{-4}$, kun $t$ mitataan vuosissa ja väkiluku
$y(t)$ miljoonissa.

(a) Ratkaise tehtävä ($y(0)=5.3$) Eulerin menetelmällä käyttämllä
askelpituussa $h=10$ 

(b) rk4:llä käyttäen n. nelinkertaista askelta (voit kokeilla pienempiäkin)

(c) Matlabin {\tt ode45}:llä.

(d) Laske analyyttinen ratkaisu Maplella (kyseessähän on
{\it 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. \verb_ http://www.math.hut.fi/teaching/v/matlab/opas.html#splinit_

\item
Tarkastellaan yht\"al\"o\"a $y'=-2\alpha(t-1)y$. 
Ratkaise aluksi analyyttisesti (saat käyttää Mapleakin.)

%$y=C\exp(-\alpha(t-1)^2)$. 

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.)


\end{enumerate}


%\verb_ 17.9.01  home/apiola/opetus/peruskurssi/v2-3/301/L/odeL4.tex _
%/home/apiola/opetus/peruskurssi/v2-3/L/latex/odeL4.tex

\subsubsection*{AA-tehtävän ratkaisun olemassolo ja yksikäsitteisyys}

{\bf Lause} $\exists_1$  (vrt. KRE s. 53 Lauseet 1 ja 2)\\
 (AA) \qquad $ y'=f(x,y), \ \ y(x_0)=y_0 $

Olkoot $f$ ja $\frac{\partial f}{\partial y}$ jatkuvia suorakulmiossa
$$ R=\lbrace(x,y) | |x-x_0|\leq a,  |y-y_0|\leq b\rbrace.$$

Tällöin AA-tehtävällä on yksikäsitteinen ratkaisu $x \mapsto y(x)$, joka
on määritelty välillä $|x-x_0| < \alpha$.

$\alpha=\textit{min}(a, \frac{b}{K})$, missä K on $|f|$:n arvojen yläraja
R:ssä: $|f(x,y)| \leq K,\ \ x\in R$

{\bf Huom!} Jatkuvuusoletuksesta seuraa, että $f$ ja $\frac{\partial f}
{\partial y}$ ovat rajoitettuja, koska suorakulmio R on suljettu.


\subsection*{Diff. yhtälöiden numeerisia menetelmiä}

\subsubsection*{1-askelmenetelmiä}
Annettu diff. yhtälö(systeemi):
$x'=f(t,x),\ \ x(t_0)=x_0$ \\
%Jos kyseessä on systeemi, niin $f$ ja $x$ ovat vektoriarvoisia.

\paragraph*{Eulerin menetelm\"a}
$x_{n+1}=x_n+h f(t_n,x_n),\ \ x(t_0)=x_0$
GKV $=O(h)$ (GKV={\bf G}lobaali {\bf k}atkaisu{\bf v}irhe) 

\paragraph*{Heunin menetelm\"a}
$x_{n+1}^{*}=x_n+h f(t_n,x_n),\ \ x(t_0)=x_0$ \\
$x_{n+1}=x_n+{\frac{1}{2}}h(f(t_n,x_n)+f(t_{n+1},x_{n+1}^*))$. \\
GKV $=O(h^2)$ 


\paragraph*{Klassinen Runge-Kutta (RK4)}


$$
\begin{cases}
k_1=h f(t_n,x_n)\\ k_2=h f(t_n+{\frac{h}{2}},x_n+{\frac{k_1}{2}}\\
        k_3=h f(t_n+{\frac{h}{2}},x_n+{\frac{k_2}{2}})\\
        k_4=h f(t_n+h,x_n+k_3) \\
       x_{n+1}=x_n+ {\frac{1}{6}}(k_1+2 k_2 + 2 k_3 + k_4) \end{cases}
$$
GKV$=O(h^4)$


\paragraph*{Taylorin menetelmä}
$y(t+h)=y(t)+hy'(t)+...+{\frac{h^p y^{(p)}}{p!}}+O(h^{p+1})$, missä
$y'(t)=f(t,y(t)), y''(t)={\frac{d}{dt}}f(t,y(t))$ (eli derivoidaan 
ketjusäännön avulla, sijoitetaan y'(t):n paikalle $f(t,y(t))$ ja jatketaan
näin.

\paragraph*{``Backward Euler''} (BE) on 
{\em implisiittisten menetelmien} kantaäiti
$x_{n+1}=x_n+h f(t_{n+1},x_{n+1})$.

\subsubsection*{Stabiilisuus}

Yhtälöä sanotaan {\tt stabiiliksi} pisteess\"a $(t,y)$, 
jos $f_y(t,y)<0$, systeemin tapauksessa $J_y$:n spektri on vasemmassa
puolitasossa (om. arvojen reaaliosat $<0$), missä $J_y$ on Jacobin matriisi
derivoituna $y$-muuttujien suhteen.

Numeerinen menetelm\"a on stabiili, jos kokonaisvirhe pysyy rajoitettuna.
T\"am\"a merkitsee Eulerin menetelmän  
%\footnote{EM=Eulerin menetelm\"a}
kohdalla ehtoa: \\ $-2 < h f_y(t,y) < 0$ 

\subsection*{Matlab-koodeja}

\verb_http://www.math.hut.fi/teaching/v/matlab/mathews/_

Listataan tämä, muita ovat mm.
 {\tt meuler.m, heun.m, milne.m, taylor.m}

\begin{verbatim}
function [T,Y] = rk4(f,a,b,ya,m)
% [T,Y] = rk4(f,a,b,ya,m)
% Runge-Kutta solution for y' = f(t,y) with y(a) = ya.
% f  is the function, input.
% a  is the left endpoint, input.
% b  is the right endpoint, input.
% ya is the initial condition, input.
% m  is the number of steps, input.
% T  is the vector of abscissas, output.
% Y  is the vector of ordinates, output.
h = (b - a)/m;
T = zeros(1,m+1);
Y = zeros(1,m+1);
T(1) = a;
Y(1) = ya;
for j=1:m,
  tj = T(j);
  yj = Y(j);
  k1 = h*feval(f,tj,yj);
  k2 = h*feval(f,tj+h/2,yj+k1/2);
  k3 = h*feval(f,tj+h/2,yj+k2/2);
  k4 = h*feval(f,tj+h,yj+k3);
  Y(j+1) = yj + (k1 + 2*k2 + 2*k3 + k4)/6;
  T(j+1) = a + h*j;
end
\end{verbatim}



\end{document}

\subsubsection*{Kankeus}
Tehtävää, jonka (Jacobin) matriisin ominaisarvot ovat negatiiviset
ja itseisarvoltaan suurimman ja pienimmän suhde on suuren suuri, sanotaan
{\em kankeaksi (``stiff'')}. Tämä on tyypillistä kemiallisille 
reaktioyhtälöille, joissa esiintyy eri suuruusluokkaa olevia aikavakioita.
Stabiilisuusvaatimus pakottaa askelpituuden kohtuuttoman pieneksi, vaikka
oltaisiin kiinnostuneita vain pitkän aikavälin ``steady state'' ratkaisusta
ja ``syyllinen'' on pienessä hetkessä pois kuoleva ``transientti'' komponentti.
Näiden käsittely vaatii aivan omat metkunsa. (\matl:n versiossa 5 on sellainen,
4-versiossa ei).

\end{document}


\subsubsection*{Matlab-avustusta}


\end{document}



%Tarkastellaan diffyhtälöitä a) $x'=x-t$ alueessa $-1 \leq t \leq 1$,
%$0 \leq x \leq 2$ ja b) $x'=x/t$ alueessa $0 < t \leq 2$, $0 < x \leq 2$.
%(Huomaa, että vaihdoimme {\it KRE}-tyylistä {\it LAODE}-tyyliseen notaatioon,
%ts $x\to t$ ja $y\to x$).

%Piirrä kummassakin tapauksessa $tx$-koordinaatistoon suuntakenttä
%ja isokliinejä ensin käsin ja sitten Matlab-skriptiä {\tt suuntak1.m}
%sopivasti editoiden ja Matlabissa ajaen. (Askel $h=0.25$ on aika sopiva,
%voit kokeilla muitakin.)
%Kommentoi jälkimmäisessä tapauksessa isokliinipiirroskomento 
%pois skriptistä.
%Miksi suuntakenttä ei erotu ennen tätä? (Perustelu kynällä ja paperilla;
%ei siis riitä: "koska siltä näyttää".)

Tämän voit (soveltuvin osin) esittää puhtaasti käsinlaskutehtävänäkin, 
jos Matlab tuntuu vielä vieraalta.
Toisaalta sitäkään en pane pahakseni, jos kokeilet LAODE-funktiota 
{\tt dfield5}.



%\item % LAODE s. 117, teht. 12: Newton's law of cooling
      % 2000 Harj. 2 teht. 1

%Newtonin jäähtymislaki sanoo, että kappaleen lämpötilan muutosnopeus
%ajanhetkellä $t$ on verrannollinen kappaleen ja ympäröivän väliaineen
%lämpötilaeroon hetkellä $t$. Differentiaaliyhtälöllä ilmaistuna:

%$$\frac{dT}{dt}=\alpha (T-T_v),$$

%missä $T_v$ on ympäröivän väliaineen vakiolämpötila ja $\alpha$ on
%verrannollisuuskerroin.

%Oletetaan, että kappaletta ympäröivän ilman lämpötila on $50^\circ$
%ja kappale jäähtyy $100^\circ$:sta $75^\circ$:een $20$ min:ssa.
%Mikä on kappaleen lämpötila tunnin kuluttua?





\item %2
Radioaktiivisessa hajoamisessa ainemäärän muutosnopeus on verrannollinen
nykyhetkellä jäljellä olevaan ainemäärään.
Muodosta tämän perusteella prosessia kuvaava differentiaaliyhtälö.
 
Radium-isotoopin $_{88}$Ra$^{224}$ puoliintumisaika on $3.6$ vrk. 
Kuinka paljon yhdestä grammasta on jäljellä a) yhden vrk:n, b)
yhden vuoden kuluttua?

{\bf Ohjetta:}
{\em Puoliintumisaika} tarkoittaa aikaa, jonka kuluttua puolet aineesta on
hajonnut. Kts. myös KRE$^8$ example 4 ss. 6--7. 




\item %5  HUB-W  

Tarkastellaan diffyhtälöitä

(a) $x'=-\frac{x}{t}$ \quad (b) $x'=\frac{x}{2t}$  \\
(c) $x'=x/t$ \quad %(d)  $x'=\frac{1-t}{1+x}$  \\

Hahmottele käsin isokliinejä ja muodosta niiden avulla suuntakentät sekä
hahmottele ratkaisukäyriä.

Selvitä, missä yhtälöt on määritelty.

Ratkaise lopuksi yhtälöt muuttujien erottelulla.

Vast: (a) $x=\frac{C}{t}$ \quad (b) $x=C\sqrt{|t|}$, 
\quad (c) $x=Ct$ \quad % (d) $(t-1)^2+(x+1)^2=C^2$



\item
Mieti (hetki), miten piirtäisit suuntakentän, kun käytössäsi on 
{\tt suuntak1.m} skriptissä saatava data. Katso sitten skriptistä
ja selvitä itsellesi (ja taululla muille), miten for-silmukassa
tapahtuva piirto tapahtuu.



\item %2 KRE 18, Verhulstin logistinen populaatiomalli.
Sopiva myös tähän kohtaan.
Harkitse !


\item KRE t. 20 CAS-project, sopiva myös, mutta ...

\item
Viime syksyn harj2 LV-tehtävissä on varteenotettavia myös, valmiiksi
ratkaistuja. Kysellään lineaarisuutta, autonomisuutta,...

