

\section{Differentiaaliyht�l�iden numeriikkaa}

\today
%18.11.05

Kirjallisuutta:

[KRE] CH 19.\\
[NR] Press-Flannery-Teukolski-Vetterling: Numerical Recipes \\
[KNM] Kahaner-Moler-Nash: Numerical Methods and Software\\
[CvL] C. van Loan\\

(T�h�n tulee t�smennyst�)
%{\tiny Luentoaiheet:

%L29 to 18.11.04  : euleresim1virhe

%}


\subsection*{Johdanto}

Differentiaaliyht�l�iden k�sittelyss� on kolme menetelm�tyyppi�:

\begin{itemize}
\item[] Analyyttiset
\item[] Kvalitatiiviset
\item[] Numeeriset
\end{itemize}

N�ill� on vahva keskin�inen vuorovaikutus. Joskus ne my�s esiintyv�t 
``valepuvuissa'', kuten kvalitatiivisen analyysin suuntakentt�piirros, johon
sijoitetaan ratkaisutrajektoreita. Oikeasti ne lasketaan numeerista menetelm��
k�ytt�en, vaikka k�ytt�j� kuvittelee olevansa keskell� puhdasta ``kvalitatiikkaa''.

Kun k�sitell��n differentiaaliyht�l�iden (systeemien) numeriikkaa, tarkasteltavana
on aivan yleinen tilanne, siisp� emme oleta autonomisuutta, emmek� muutakaan 
erityisominaisuutta.

Tarkastellaan differeniaaliyht�l�n alkuarvoteht�v��:

%\begin{equation}\label{ode1}
$$\y'=\f(t,\y),\ \ \y(0)=\y_0.$$
%\end{equation}

P��paino t�ss� on yleisell� 1. kertaluvun diffyht�l�systeemill�.
Kuten tunnettua, jokainen korkeamman kertaluvun yht�l� tai systeemi
voidaan palauttaa 1. kertaluvun systeemiksi ottamalla derivaatat uusiksi
muuttujiksi.
Yll� oleva yht�l� edustaa siten yht� hyvin yht� skalaariyht�l�� kuin yht�l�systeemi�.
J�lkimm�isess� tapauksessa
ymm�rr�mme $\y$:n vektoriksi ja $\f$:n vektoriarvoiseksi 
funktioksi.

%saamme diffyht�l�systeemin aivan samaan muotoon \ref{ode1}


Mit� tapahtuu, kun {\tt pplane}-kuvassa painnalletaan hiirt� jossain 
suuntakent�n
pisteess�. Numeerinen ratkaisija laskee ja piirt�� ratkaisuapproksimaation.

%Numeeriset menetelm�t diffyht�l�ille ovat ��rimm�isen t�rkeit� ...

Periaatteessa edet��n alkuarvopisteest� pienin askelin, joko tangentin
suuntaan (minne nen� osoittaa) tai tehd��n "maaston tunnustelua" erilaisilla
koeaskelilla ja edet��n sopivaan painotettuun keskiarvosuuntaan.

Edellisen puhdasoppisin edustaja on \emph{Eulerin menetelm�}, j�lkimm�isen
tyyppiesimerkki on \emph{Runge--Kutta}-menetelm�, josta hetken kuluttua hiukan
lis��.

\begin{center}
\begin{tabular}{|l|l|}
\hline
\textbf{Menetelm�tyyppej�} & \textbf{Yht�l�iden/menetelmien ominaisuuksia}\\ \hline
Yksiaskel/moniaskel  & Stabiilisuus/h�ri�alttius \\
Ekspilisiittiset/implisiitiset & kankeus.\\
\hline
\end{tabular}
\end{center}




%\begin{itemize}
%\item
% \emph{Yksiaskel/moniaskel}menetelm�t 
%\item
% \emph{Eksplisiittiset/implisiittiset}
%\end{itemize}

Toteustustekninen t�rke� periaate on \textbf{adaptiivisuus}, jolla tarkoitetaan
automaattista askelpituuden s��t�� paikallisen virhearvion avulla siten,
ett� tavoitteeksi astettu virhetoleranssi saavutetaan niin, ett� laskentaty�t�
s��stet��n "�lykk��ll�" askelstrategialla.

\subsection{Ratkaisun olemassaolo ja yksik�sitteisyys}



\subsection{Yksiaskelmenetelm�t}

Kuten aina, aloitetaan \emph{Eulerin menetelm�st�}.

Jos $y(t)$ on ratkaisufunktio, muodostetaan \emph{Taylorin} kehitelm�
pisteess� $t$:

$y(t + h)=y(t)+hy'(t) + O(h^2) = y(t)+hf(t,y(t))+O(h^2).$

Jos pudotetaan j��nn�stermi, saadaan siten approksimaatio ratkaisufunktiolle
pisteess� $t+h,$ kun tunnetaan ratkaisu(approksimaatio) pisteess� $t.$

N�in johdumme \\
\textbf{Eulerin menetelm��n}:\\
Annettu alkupiste: $(t_0, \y_0).$ Lasketaan 
$\y_{k+1}=\y_k + h \f(t_k,\y_k), k=0,\ldots, n$



%\section*{Lis�ohjeita harjoituksiin 9}

%14.11.2004


%\textbf{Eulerin menetelm�}

%\begin{verbatim}
%\end{verbatim}

%Ratkaistavana on alkuarvoteht�v�:
%\begin{equation}\label{ode1}
%\vec{y}'=\vec{f}(t,\vec{y}),\ \ \vec{y}(0)=\vec{y}_0.
%\end{equation}



%\[
% \ \  \ \ \vec{y}'=\vec{f}(t,\vec{y}),\ \ \vec{y}(0)=\vec{y}_0.
%\]


%Voidaan ajatella joko skaariyht�l�n� tai vektoriyht�l�n�, 
%eli 1. kertaluvun yht�l�systeemin�.

%\textbf{Algoritmi} \\
%Alkuarvo $(t_0,\y_0)$\\
%\y_{k+1}=\y_k + h\,\f(t_k,\y_k), k=1,2,\ldots 


\begin{esim} Oikein kunnolla ep�lineaarinen ja ep�autonominen:

$$y'=\sin(t\, y), \ \ y(0)=3.$$

T�ss� siis $f(t,y)=\sin(t\,y).$
Olkoon $h=0.1$. J�rjestet��n laskut taulukon muotoon:
\begin{center}
\begin{tabular}{|c|c|l|}
\hline      %1tabular1taulukko1latex
$t_n$ & $y_n$ & $y_{n+1}=y_n + h\,f(t_n,y_n)$ \\
\hline        
$t_0=0.0$ & $y_0=3.000$ & $y_1=3+0.1\,\sin(0\cdot 3)=3.000$\\ 
$t_1=0.1$ & $y_1=3.000$ & $y_2=3+0.1\,\sin(0.1\cdot 3)=3.030$\\ 
$t_2=0.2$ & $y_2=3.030$ & $y_3=3.030+0.1\,\sin(0.2\cdot 3.03)=3.087$\\ 
$t_3=0.3$ & $y_3=3.087$ & $y_4=3.087+0.1\,\sin(0.3\cdot 3.087)=3.167$\\ 
$t_4=0.4$ & $y_4=3.167$ & \ldots \\
\hline

\end{tabular}
\end{center}
\end{esim}


Esimerkki voitaisiin tehd� \matl:lla vaikka n�in:

\begin{verbatim}
clear; close all
h=0.1; n=10;
t=0:h:1;
y(1)=3;   % Indeksointi alkaa 1:st� Matlabissa.
for k=1:n
    y(k+1)=y(k)+h*sin(t(k)*y(k));
end;
[t;y]
\end{verbatim}
N�ytet��n tuloksesta 7 ensimm�ist� arvoa
\begin{verbatim}
ans =
  Columns 1 through 7
         0    0.1000    0.2000    0.3000    0.4000    0.5000    0.6000 
    3.0000    3.0000    3.0296    3.0865    3.1664    3.2618    3.3616 
\end{verbatim}
Tulosten pienet eroavuudet johtuvat siit�, ett� "k�sinlasku" on tehty $4:$n numeron tarkkuudella, kun taas \matl\, laskee n. 
$16:$lla numerolla, joista n�ytet��n oletusarvona $5$ ja systeemikomennon \verb_ >> format long_ j�lkeen kaikki $16.$
 
Voidaan tietysti jatkaa piirt�en, edustakoon \matl-koodi kuvaa.
\begin{verbatim}
plot(t,y)
ylim([2.8 3.6])
\end{verbatim}


\subsection{Diffyht�l�systeemit}

Kuten todettu, tulkitaan vain $\y$ ja $\f$ yll� vektoriarvoisiksi.
Otetaan aluksi esimerkki lineaarisesta systeemist�.
(Kiirehdit��n senkin takia, ett� AV-harjoituksissa on t�llainen ja luennolla
tuli v�h�n pikaisesti.)

\begin{esim} Lineaarinen diffyht�l�systeemi.

$$\y'=A \,y, \ \ \y(0)=\lbrack 1,0\rbrack.$$

T�ss� siis
$\f(t,\y)= A\, \y$ (systeemi on autonominen, joten $t$ ei oikeasti esiinny
oikealla.)
Siten iteraatioaskel on


$$\y_{n+1}=\y_n + h\, \f(t_n,\y_n) = \y_n + h\, A\,\y_n = (I+hA)\y_n.$$

Toisin sanoen iteraatio on yksinkertaisesti  matriisilla $B=I+hA$ kertomista.
\matl:lla vaikka t�h�n tapaan:

\begin{verbatim}
clear
A=[1 -4;-1 1];
h=0.2; n=5;                % Askelpituus ja askelten lukum��r�. 
B=eye(size(A))+h*A;        % Iteraatiomatriisi.

Y=zeros(2,n);   % Alustetaan 2-rivinen matriisi, jonka sarakkeiksi sijoitetaan 
                % Eulerin iteraation tuottamat pisteet. (Alustus ei ole v�ltt�m�t�n,
                % mutta lienee selkeytt�v�, lis�ksi parantaa tehokkuutta.)
Y(:,1)=[1;0];   % Alkuarvo sijoitetaan 1. sarakkeeksi. (Matlabin indeksointi alkaa 1:st�)

for k=1:n
    Y(:,k+1)=B*Y(:,k);   % Iteraatioaskel
end;
Y              % Katsotaan, mit� Y-matriisin sarakkeissa on.
Y =
    1.0000    1.2000    1.6000    2.3040    3.4816    5.4067
         0   -0.2000   -0.4800   -0.8960   -1.5360   -2.5395
\end{verbatim}
Sarake 1: alkuarvo, sarake 2: ensimm�inen iteraatio, $\ldots$
\end{esim}



Seuraavaksi katsotaan tuttua, esimerkki� reippaasti ep�lineaarisesta
tapauksesta, nimitt�in heiluria.
\begin{esim}
Heiluriyht�l� $y''=-k\sin(y)$ muunnetaan 1. kl:n systeemiksi:

$$
\begin{cases}
y_1'=y2\\
y_2'=-k\, \sin\, y_1
\end{cases}
$$

Lasketaanpa vaikkapa alkupisteest� $(0,1)$ l�htien muutama askel.

$\y_0=\left[ \begin {array}{c} 0\\\noalign{\medskip}1\end {array} \right] $,\quad $h= 0.1$\\
$\f(\y_0)=\left[ \begin {array}{c} 1\\\noalign{\medskip}0\end {array} \right] $

$\y_1= {\y}_{{0}}+ 0.1\,f \left( {\y}_{{0}} \right)= \left[ \begin {array}{c}  0.100\\\noalign{\medskip}
 1.0\end {array} \right] $

$\y_2=
 \left[ \begin {array}{c}  0.200\\\noalign{\medskip}
 0.990\end {array} \right] 
$


\end{esim}



\end{document}

%\subsubsection{Eulerin menetelm�n Matlab-toteutus}


%Ensin 

%heisys=inline('[y(2);-sin(y(1))]','t','y')




\subsection*{Stabiilisuus}

Helpointa on suorittaa analyysi Eulerin menetelm�n avulla. 
Muut eksplisiittiset yksiaskelmenetelm�t ovat periaatteessa samantyylisi�, 
samat k�sitteet esiintyv�t, kaavat hiukan mutkikkaampia.

Eulerin menetelm�lle saadaan suoraan m��ritelm�st� ja soveltamalla v�liarvolausetta:

Globaali virhe 
$t_{k+1}$:ss� = $(1+ h J)$(Globaali virhe $t_k$:ssa) + lokaali virhe 
$t_{k+1}$:ss�.

T�ss� $J$ on Jacobin matriisi muuttujien $y_k$ suhteen, yhden yht�l�n 
tapauksessa $J=f_y$. $J=J(\xi)$, siis $J$ lasketaan jossain pisteess�,
joka yhden muuttujan tapauksessa on $y_k$:n ja $y_{k+1}$:n v�liss�.
(kts [KNM] s. 228) 
Yht�l�systeemin tapauksessa luku 
$1$ on sopivissa kohdin tulkittava yksikk�matriisiksi.

Usein n�it� kysymyksi� tarkastellaan prototyyppiyht�l�iden

$y' = Ay$ avulla, vrt.. esim. [CvL],  [NR]






\subsection*{Kankeus}


$${\frac {d}{dt}}y \left( t \right) =-1000\,y \left( t \right) +1000\,
\sin \left( t \right) +\cos \left( t \right) $$


$$
\left[ \begin {array}{cc} 998&1998\\\noalign{\medskip}-999&-1999\end {array} \right] 
$$

$$
\left[ \begin {array}{cc} -{e^{-1000\,t}}+2\,{e^{-t}}&2\,{e^{-t}}-2\,{e^{-1000\,t}}\\\noalign{\medskip}-{e^{-t}}+{e^{-1000\,t}}&2\,{e^{-
1000\,t}}-{e^{-t}}\end {array} \right] 
$$

$$
\left[ \begin {array}{c} -{e^{-1000\,t}}+2\,{e^{-t}}\\\noalign{\medskip}-{e^{-t}}+{e^{-1000\,t}}\end {array} \right] 
$$