\documentclass[finnish,11pt,twocolumn]{article}
\usepackage{/home/apiola/tex/numsym04}
\usepackage[dvips]{epsfig}
%\include{/home/apiola/tex/defs}
%\usepackage{v23harj}
%\include{defs}
\def\abs#1{\vert #1\vert}      
\def\norm#1{\Vert #1\Vert}      
\def\innerp#1{\langle #1 \rangle}
\def\vnormone#1{\sum_k \vert #1_k \vert}
\def\vnorminf#1{\max_k \vert #1_k \vert}
\def\mnormone#1{\max_j \sum_i \vert#1_{ij}\vert}
\def\mnorminf#1{\max_i \sum_j \vert#1_{ij}\vert}
\def\rv{\hfill\break}
\def\matl{{\sc Matlab}}
\def\maple{{\sc Maple}}
%\def\Arg#1{\textit{Arg} #1}
\def\Arg{{\textit{Arg}}\ \ }
\def\huom {{\bf Huom!}}

\begin{document}

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

\harjoitus{3}{6}{3 -- 5.2.2004}

%\begin{verbatim}
% http://www.math.hut.fi/teaching/v/3/L/fourier.html
% http://www.math.hut.fi/teaching/v/3/L/fourier.mws
% http://www.math.hut.fi/teaching/numsym/02/H/harj1ohje.mws
%\end{verbatim}


Aktivoidaan Matlab:n k�ytt�.

{\bf Harj. 2:n ``r�stiteht�v�t''}: J�tet��n teht�v� 4 omaehtoisen opiskelun
piiriin, luovutaan t�ll� er�� Maplen {\tt dsolve(...,numeric)}-panostuksesta
ja keskityt��n numeerisessa ratkaisussa Matlabiin.

Sensijaan H2t4 on hyv� malli sille, miten lineaarisia Maplella k�sitell��n
analyyttisten ratkaisujen suhteen. Sen mukaan on sopivaa tehd� teht�v� 5,
jota t�ydennet��n Matlabin pplane5-hommalla ja t�m�n harjoituksen
teht�v�ll� 4. 


\begin{enumerate}

\item %1
Tutustu tai palauta mieleen pinta- ja korkeusk�yr�piirrokset
Matlabissa. Lue viitteiss� olevista oppaista ja katso \\
{\tt help mesh, meshgrid, surf, contour}. Kokeile my�s
vuorovaikutteista, hiirell� valittavaa korkeusk�yr��.

Sopiva testikuvattava on harj2:n teht�v�n 2 funktio.

Piirr� lis�ksi gradienttikentt� {\tt quiver}-funktiolla.
K�yt� Maplea gradienttikent�n muodostamiseen. Maplesta saat
kaavat Matlabiin {\tt lprint}:ll�, toisaalta voit halutessasi kokeilla
my�s symbolic toolboxin funktioita Matlab-komennolla {\tt maple}.

\item %2
Perehdy 
\verb_ http://www.math.hut.fi/~apiola/odense2000/v2/optmlb.html_-sivuun. 
Rakenna sen avulla Matlab-skripti, jolla
voit k�sitell� optimointimenetelmi� ja my�s havainnollistaa
sen vaiheita ja visualisoida monipuolisesti.

Sovella taas harj.2 teht. 2:een.
Voit ottaa SD-menetelm�n lis�ksi my�s Newtonin menetelm� mukaan. (Alku SD:ll�,
loppu Newtonilla-tyyliin.) (T�t� lis�yst� voidaan pit�� ``vapaaehtoisena''.)

K�yt� Maplea lausekkeiden generointiin.

Mieluusti saat kokeilla muitakin funktioita. 


\item  %3
Kirjoita Matlab-koodi, jolla saat autonomisen $2\times 2-$ 
diffyht�l�systeemin suuntakent�n piirretyksi.
T�m� on samantapainen teht�v� kuin gradienttikentt� yll�.
T�ss� tapauksessa automaattinen skaalaus voi olla haitallinen, sen
voi poistaa antamalla {\tt quiver}:lle argumentin $0$ sopivaan kohtaan.

Muodosta harj. 2:n teht�v�n 5 suuntakentt�. 

\item %4
Laske ja piirr� joitakin harj. 2:n teht. 5:n ratkaisuk�yri� sek� 
aikakuvaan
ett� faasitasoon. K�yt� funktiota {\tt ode45}.

Yhdist� suuntakentt� ja ratkaisuk�yr�piirros. Voit valita ratkaisuk�yri�
hiirell� hy�dynt�m�ll� {\tt ginput}-funktiota.
Er�s idea:
\begin{verbatim}
ya=[0;0]  % "Alkualkupiste"
while length(ya)>0
  [T,Y]=ode45(...,ya);
  plot(...)
  ya=ginput(1);   % Valitaan hiirell� piste.
end; 
\end{verbatim}

\end{enumerate}

T�ydennet��n viel� lis�ohjeilla.


\end{document}



\subsection*{Matriisieksponenttifunktio $e^{At}$}

Tiivis, hyv� esitys aiheesta on {\it Pekka Alestalon} kirjoittama:\\
{http://www.math.hut.fi/teaching/k3/expm.ps}

Jokaiselle neli�matriisille $A$ voidaan m��ritell� eksponenttifunktio
$e^{A}:=\sum_{k=0}^\infty {\frac{A^k}{k!}}$. Diagonalisoituvalle matriisille
se voidaan laskea suoraan soveltamalla diagonalisointiesityksen $D$-matriisin
diagonaalialkioihin tavallista $\exp$-funktiota.
Joskus voi olla k�tev�� soveltaa sopivaa reaalifunktion sarjakehitelm��,
joskus taas yll� mainittu m��ritelm� johtaa katkeavaan sarjaan yll�tt�v�sti
syyst�, jota lukujen maailmassa ei esiinny, nimitt�in matriisin jokin potenssi
voi olla nolla (nilpotentti matriisi). 


Miten sitten toimitaan yleisesti, kun ei kuitenkaan haluta k�ytt�� 
sarjakehitelm��? Silloin etsit��n ns. yleistetyt ominaisvektorit, mik�
vastaa sit� tapaa, joka esitettiin $2\times 2 $-matriiseille.
Yleisess� muodossa esitettyn� se voidaan muotoilla matriisihajotelmaksi,
joka on ``l�hes diagonalisointi'', eli ns. {\it Jordanin muoto}.
Sit� emme t�ll� kurssilla l�hemmin opiskele, kiinnostuneet l�yt�v�t
sen mm. kelpo Laode-kirjasta.

Algoritmi voidaan esitt�� ilman Jordan-hajoitelma-terminologiaa (vaikka
samasta asiasta on kyse).

Itse asiassa p�tee ns. {\em Cayley-Hamiltonin} lause: Jokainen neli�matriisi
toteuttaa oman karakteristisen yht�l�ns�.
(Tarkoittaa, ett� jos lausekkeessa $p(\lambda)=\lambda^n+\ldots $ sijoitetaan
$\lambda$:n paikalle $A$, niin tuloksena on nollamatriisi.)
 T�st� p��dyt��n siihen, ett� matriisi
``k�ytt�ytyy nilpotentisti'' yleistettyjen ominaisvektorien edustamilla
suunnilla. (Kts. tarkemmin alla mainituilta sivuilta).

Helpoimmat tapaukset ovat ��ritapaukset: \newline 
a) Jos ominaisvektoreita on
riitt�v�sti ($n$ kpl) tai \newline
b) jos niit� on vain yksi (jolloin ominaisrvojakin
on vain yksi). \newline
T�ll�inh�n yll� mainitusta {\em Cayley-Hamiltonista} seuraa,
ett� jos $\lambda_0$ on tuo ominaisarvo, niin $(A-\lambda_0)^n=0$.
Koska $e^{At}=e^{\lambda_0 t}e^{(A-\lambda_0)t}$, saadaan $(A-\lambda_0I)t$:n
potenssien mukainen sarja, joka katkeaa n:nness� termiss�.

$2\times 2$-matriisien k�sittely on erityisen helppoa, koska niill� joko
on riitt�v�sti (2) ominaisvektoreita tai $(A-\lambda_0)^2=0$

Algoritmi on selitetty tarkemmin luentosivuilla: (Emme sit� kuitenkaan kokeissa vaadi osattavaksi.)


\verb+www.math.hut.fi/teaching/v/3/L/L_ldys.html+ ja t�m�n viitteiss�:\hfill\\
\verb_www.math.hut.fi/teaching/k3/luentomateriaali/expA.html_ ja erit. \hfill\\
\verb_www.math.hut.fi/teaching/k3/luentomateriaali/eat8.gif_ \\
\verb_www.math.hut.fi/teaching/k3/luentomateriaali/eat9.gif_

K�yt�nn�ss� voidaan turvautua ohjelmistoihin, joihin tuo Jordan-muoto
on k�tketty. Maple:ssa on {\tt linalg[exponential]} ja Matlabissa {\tt expm}. 

%\matl:n tapauksessa saamme
%tietysti ``numeerisen tuloksen nopeasti'' ja sanoisimmeko, ett� \maple:n
%tapauksessa ``symbolisen tuloksen hitaasti''


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

%\input{/home/apiola/tex/AmsLTXfiMacro.tex}

\harjoitus{9}{46}{10 --14.11.2003}

{\bf Kirjallisuusviitteit�: }
%KRE8 Ch 3, KRE7 Ch 4, Lop Ch 12. 

%KRE8: 5.3, 5.5, 5.7, KRE7 6.3,6.4,6.6,6.7,
%Lop 6.6, 6.9, 6.10\\
%Luentosivuja: {\tt www.math.hut.fi/teaching/k3/03/L/}. 
%Matlab-oppaat:
%\verb_ http://www.math.hut.fi/~apiola/matlab/opas/lyhyt/_ \\
%\verb_ http://www.math.hut.fi/~apiola/matlab/opas/ei-niin-lyhyt/ _

\subsubsection*{Alkuviikko (AV)}

\begin{enumerate}

\item %2
M��rit� systeemin
$$\begin{cases}
        % \dot x=x-x^2+2y+3 \\
         %\dot y=-x+2y 
    x'=x-x^2+2y+3 \\
    y'=-x+2y 
\end{cases}$$
kriittiset pisteet.

Linearisoi systeemi kunkin kriittisen pisteen suhteen ja m��rit� niiden luonne
sek� hahmottele faasikuvaa.

%Ratk. Jeffrey s. 579--580


\item
Klassinen saalis-saalistaja-ongelma ("ketut ja j\"anikset", Volterran-Lotkan
malli) on seuraavanlainen:  Alueella olevien j�nisten ja kettujen
lukum��ri� merkit\"a\"an $x(t)$:ll\"a ja $y(t)$:ll\"a vastaavasti.

Ajatellaan, ett� j�nikset sy�v�t ruohoa ja ketut j�niksi� (ja vain niit�).
Kasvua rajoittavia tekij�it� ei ole, muuta kuin j�niksille ketut
ja ketuille j�nisten puute (eli kilpailevat ketut). Tilannetta mallintavat 
yht�l�t voisivat olla
vaikkapa:
$\begin{cases}
x'= 2 x -\alpha x y \\
y'=-y  + \alpha x y.
\end{cases} (\alpha > 0).$

M��rit� systeemin kriittiset pisteet ja linearisoi niiss�.
M��rit� linearisoitujen systeemien luonne.

\item
Tarkastellaan vaimentamatonta heiluria, jonka yht�l� on 
$\Theta''+k\sin(\Theta)=0.$
Kun se kirjoitetaan 1. kertaluvun systeemiksi, voidaan ratkaista
kriittiset pisteet ja suorittaa linearisointi, piirt�� faasitasokuva jne.

T�m� on tehty mm. KRE8 Exa 1 ss. 176--177 (Linearisointi voitaisiin
hyvin tehd� alla neuvotulla Jakobiaani-tyylill�, joka on yleisp�tev�mpi.)

Katso KRE-kirjan kuvaa tai piirr� {\tt pplane5}:ll� (valitsemalla 
{\em Gallery} ja sielt� {\em pendulum} (anna $D$:n olla oletusarvossaan $0$)).

Selvit� heilurin tila (paikka, nopeus, liikkeen suunta) \\
(a) Umpinaisen trajektorin ja vaaka- ja pystyakselin leikkauspisteess�
(ota mukaan my�s ``rajatapaustrajektori, joka kulkee pisteen $(\pi,0)$
kautta).

(b) Aaltoilevaa (ei umpinaista) trajektoria seuraillen v�lill� 
$\left[0,\pi\right]$.


\item
Sovella Eulerin menetelm�� lineaariseen systeemiin
$\begin{cases} 
y_1'=-y_1+y_2\\
y_2'=-y_1-y_2
\end{cases} y_1(0)=0,\ \ y_2(0)=4.$

K�yt� askelpituutta $h=0.2$ ja laske $5$ askelta.
Piirr� ratkaisupolku $y_1y_2-$ tasoon (faasitasoon) ja vertaa tarkkaan
ratkaisuun.
 
\item
Sovella {\em Eulerin} menetelm�� yht�l��n $y''=t y' -3 y, \ \ y(0)=0,
\ \ y'(0)=-3.$
Laske $5$ askelta, $h=0.05$. Totea, ett� tarkka ratkaisu on
$y(t)=t^3 - 3t,$ ja laske virhe.

Ohje: Kirjoita yht�l� 1. kertaluvun systeemiksi.

\item
Laske teht�v� 4 Runge-Kutta4-menetelm�ll� ottamalla vain yksi $h=1$:n
pituinen askel. Vertaa teht�v�n 4 tulokseen.

\end{enumerate}


\subsection*{Loppuviikko}


\begin{enumerate}

\item %1
Osoita, ett� systeemin
$
\begin{cases}
x'=y + x(x^2+y^2)\\ 
y'=-x +y(x^2+y^2) 
\end{cases}$
ainoa KRP on $(0,0).$ Linearisoi systeemi $O$:ssa ja selvit� sen
tyyppi. 

Piirr� faasikuva {\tt pplane5}:ll�. Piirtele muutama trajektori.
Valitse ``PPLANE5 Display''-ruudun
``solutions''-valikosta ``Find an equilibrium point''. Kuljeta hiusristikkoa
l�helle origoa ja klikkaa. Vertaa saamasi ``Equilibrium point''-ikkunan
antamia tuloksia omiin laskuihisi ja valitse ``Display the linearization''.
Piirr� linearisoituun joitakin trajektoreita.

Antaako $O$:ssa linearisoitu systeemi oikean vai v��r�n tyypin 
alkuper�iselle systeemille siin�?


\item %2
Eulerin menetelm� 1-ulotteisessa tapauksessa.
Voit operoida joko {\tt L/dynum1.m} tyyliin tai n�in:

Ota tiedosto {\tt H/euleri.m} ja sijoita se Matlab-polun varteen.
(Koodi annettu my�s alla.)

Jos ratkaistavana olisi yht�l� $y'=t-y^2$ ja alkuarvona vaikkapa $y(0)=1$,
voisit kokeilla vaikka tyyliin:
\begin{verbatim}
>> clear; clf
>> f1=inline('t-y.^2','t','y')  % euleri-koodista n�kyy, ett� y^2 toimii
>> n=2;hold on; grid on            % yht� hyvin (mutta hyv� totutella(.^))
>> [T,Y]=euleri(f1,0,1,1,n);n=2*n,plot(T,Y,T,Y,'or');shg
\end{verbatim}
(Pelkk� $f$ voi olla vaarallinen funktionnimi, niit� voi siunaantua eri paikkoihin. Kannattaa tarkistaa: {\tt which f} ja {\tt type f}.)

Iteroi alimmaista rivi� nuolin�pp�imell� ja tee havaintoja.
Kun $n=32$ (sit� luokkaa), niin lakkaa piirt�m�st� rinkuloita (editoi
piirtokomentoa).
Huomannet my�s prosessin hidastumisen.

Tarkkaa ratkaisua ei ole helppo laskea, eik� onnistukaan alkeisfunktioiden
avulla. Virhek�yt�st� voidaan tutkia vertaamalla Matlabin valmiiseen
ode-ratkaisijaan {\tt ode45} vaikkapa.

Kutsu voisi olla:

\verb_>> [T1,Y1]=ode45(f1,[0,1],1);_

Tutki virhetta pisteess� $t=1$ ottamalla ratkaisun arvoksi {\tt Y1(end)}.
Muodosta siis erotuksia ja katso, onko virhek�yt�s tyyppi� $O(h)$, eli
askeleen puolittuessa virhe puolittuu (kutakuinkin).


%0.8333833913


%[T,Y]=euleri(f1,0,1,1,20); plot(T,Y)

"Vapaaehtoinen" lis�teht�v�: Tutki Heunin ja Runge-Kutta4-menetelmien virhett� vastaavaan
tapaan. 

\item %3
%s97 harj9AV teht. 6
Totea, ett� $O$ on systeemin 
$\begin{cases} 
x'=-\frac{3}{2}x -10(1+x)y \\
y'=5x - \frac{1}{2}y
\end{cases}$
kriittinen piste ja linearisoi systeemi siin�.
Mik� on KRP:n tyyppi ja stabiilisuus?

Piirr� asianmukainen {\tt pplane5}-piirros, tarkista k�sinlaskusi vertaamalla ohjelman 
laskemaan linearisointiin.

Ratkaise sitten AA-teht�v� $x(0)=1, y(0)=1$ Eulerin menetelm�ll� ({\tt eulers.m})
k�ytt�en askelpituutta
$h=0.1.$ Piirr� faasimurtoviiva. Voit piirt�� sen samaan kuvaan pplane-piirroksen kanssa
(voi olla viisasta ensin harjoitella eri ikkunassa, muista {\tt figure(n)}).
Huomannet, ett� menetelm� k�ytt�ytyy jopa kvalitatiivisesti v��rin.
Voit toki kokeilla pienemm�ll� askeleella.

P��asia on verrata parempiin menetelmiin, {\em Heun, RK4} ja "adaptiivinen" RK45.
Viimemainitun tuloksen n�et klikkaamalla alkupistett� hiirell� tai soveltamalla
Matlab-funktiota {\tt ode45}.

Voit tehd� tyylikk��n kokonaisuuden kuvineen ja johtop��t�ksineen (omaksi opiksesi
ja iloksesi).


\item %4
Sovella RK4-menetelm�� vaimentamattoman heilurin tapaukseen: 
$y''+\sin y = 0, \ \ y(\pi)=0, y'(\pi)=1.$

Ota sopiva aikav�li ja kokeile askelpituutta. Piirr� trajektori faasitasoon, voit 
hyvin sijoittaa sen pplane5-kuvaan, kuten ennenkin.


%\item
%Vaimennettu heiluri

%\item
%Laske teht. AV 5 sek� Eulerilla ett� RK4:ll�.
%Virhek�yt�s?


\end{enumerate}

\subsection*{Ohjeita}

{\bf Linearisointi} tapahtuu laskemalla kriittiset pisteet.
Jos $ {\bf p_0}$ on KRP, niin linearisoidun systeemin matriisi on 
diff.yht�l�n m��rittelev�n (vektoriarvoisen) funktion ${\bf f}$
Jacobin matriisi 
$({\frac {\partial}{\partial y_{{j}}}}f_{{i}}(y_1,...y_n)))$

Jos linearisoidun systeemin O:n luonne on keskus (ominaisarvot puhtaasti
imaginaariset), se ei v�ltt�m�tt� kerro alkuper�isen ep�lineaarisen
systeemin KRP:n luonnetta. Muissa tapauksissa ($\Re(\lambda)\neq 0$)
sensijaan kertoo. 

Numeerisia menetelmi� laskiessasi voit hyvin k�ytt�� tavallisen laskimen 
sijasta Matlabia laskimena. Tarkoituksenmukaisin tapa lienee m��ritell�
systeemin oikea puoli (vektori)funktioksi, jolloin laskeminen on pelkk��
riemua!
Esim:
\begin{verbatim}
>> F=inline('[-y(1)+y(2);-y(1)-y(2)]','y')
>> h=0.2;
>> y0=[0;4]
>> y1=y0+h*F(y0)
...
\end{verbatim}

{\bf KRP-tyypit.} Linearisointitilanteessa on samassa teht�v�ss�
m��r�tt�v� usein monta KRP-tyyppi�. T�ll�in tr(A),det(A)-ehdot ovat
vaivaa s��st�vi�. Vapautamme ne edellisess� teht�v�paperissa julistetusta
pannasta. Emme edelleenk��n silti oikein sulata ``p'' ja ``q''-nimityksi�,
koska ne eiv�t ole ``luonnollisia'' etenkin, kun p sis�lt�� miinusmerkin.

\subsection*{Diff. yht�l�iden numeriikkaa}

\textbf{Eulerin menetelm�}. Yhden skalaariyht�l�n tapaus voidaan toteuttaa
Matlab-funktiolla:

\begin{verbatim}
function [T,Y]=euleri(fnimi,a,b,ya,n)
% fnimi : funktiom��ritys, inline tai m-tiedosto 
%          a: alkupiste, b:loppupiste, ya:alkuarvo (y(a)=ya)
%          n: v�lin (a,b) jakov�lien lukum��r�
% Esim: f1=inline('t-y^2','t','y'); 
%       [T,Y]=euleri(f1,0,1,1,20); plot(T,Y)
h=(b-a)/n;
Y=zeros(1,n+1);T=Y;
T(1)=a;Y(1)=ya;
for j=1:n
   Y(j+1)=Y(j)+h*feval(fnimi,T(j),Y(j));
   T(j+1)=a+j*h;
end;
\end{verbatim}

Vektoriversiot systeemeille ovat {\tt .../matlab/eulers.m, heuns.m, rk4s.m}

Matlabin omat ode-ratkaisijat ovat varsin laadukkaat ja monipuoliset. 
{\tt help ode45}-komennolla p��set alkuun. 

%Viel� nopeammin seuraavien esimerkkien avulla:

Lis�� ohjeita: {\tt H/harj9LVohje.html}

%Perusesimerkki

%systeemille:

%s97 harj9 Numeerisia menetelmi� ...

\subsection*{Diff. yht�l�iden numeeristen menetelmien pikaesittely}


\subsubsection*{1-askelmenetelmi�}

T�ss� merkit��n ratkaistavaa funktiota $x(t)$:ll� ($y(t)$:n sijasta), jos
kyseess� on vektorifunktio, niin notaatio on sama.

Annettu diff. yht�l�(systeemi):
$x'=f(t,x),\ \ x(t_0)=x_0$ %\rv

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$ \rv
$x_{n+1}=x_n+{\frac{1}{2}}h(f(t_n,x_n)+f(t_{n+1},x_{n+1}^*))$. \rv
GKV $=O(h^2)$ 


\paragraph*{Klassinen Runge-Kutta (RK4)}
$$\begin{cases}k_1=h f(t_n,x_n)\cr 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: \rv $-2 < h f_y(t,y) < 0$ 

\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}

clear; clf
f=inline('t-y^2')
h=0.1;
t0=0;y0=1;T=[];Y=[];
t=t0;y=y0;
hold on
%%%%%%%%%%%%%%%%%%%%%%%
t=t+h;y=y+h*f(t,y);T=[T,t];Y=[Y,y];plot(T,Y);plot(T,Y,'or');shg

function [T,Y]=euleri(fnimi,a,b,ya,n)
h=(b-a)/n;
Y=zeros(1,n+1);T=Y;
T(1)=a;Y(1)=ya;
for j=1:n
   Y(j+1)=Y(j)+h*feval(fnimi,T(j),Y(j));
   T(j+1)=a+j*h;
end;

>> [T,Y]=euleri(f,0,1,1,10);
>> plot(T,Y);hold on;plot(T,Y,'or');grid;shg

clf
n=2;hold on
[T,Y]=euleri(f,0,1,1,n);n=2*n;plot(T,Y,T,Y,'or');grid;shg

[T,Y]=ode45(f1,[0,1],1);
 plot(T,Y,'r')




\end{document}







\end{document}







