\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{2}{5}{27 -- 29.1.2004}

Versio to 22.1., pari ode-tehtävää lisätty.\\
Tarkoitus: valmiiksi to 29.1. mennessä. 
(Täsmennetään aikataulu ti).

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


\begin{enumerate}

\item %1
Suorita pari "Steepest descent"-askelta kohdefunktiolle 
$f(\V{x})=2(x_1^2+x_2^2)+x_1x_2 -5(x_1+x_2)$, kun alkupisteenä on 
$\V{x_0}=[1, -2]^T$.


Kts. {\tt maple/harj2ohjeita.mws} ja alla olevat ohjeet.
(vrt. harj. 1, teht. 3). Tällä kertaa käytetään funktioita lausekkeiden sijasta.

Suorita ohjetyöarkilla ehdotetut visualisoinnit.

\item %2
Piirrä funktion 
$f(x,y)=e^x(4x^2+2y^2+4xy+2y+1)$
korkeuskäyrä- ja pintapiirros. Piirrä myös gradienttikenttä plots-pakkauksen 
\texttt{fieldplot}:lla.

Määritä $f$:n kriittiset pisteet (gradientin 0-kohdat) sekä niiden luonne.
Voit käyttää {\tt linalg}-pakkauksen funktioita {\tt grad} ja {\tt hessian}.
Tässä tapauksessa saattaa jopa {\tt solve} osata, yleensä joudutaan turvautumaan
{\tt fsolve}:een.
Laske Hessen matriisin ominaisarvot kriittisissä pisteissä ja päättele "luonne"
niiden avulla. Käytä {\tt LinearAlgebra}-pakkauksen {\tt Eigenvalues} tai 
{\tt Eigenvectors}-komentoa. (Opiskele ..maple/LA.mws:ää)

Sovella SD-menetelmää funktion
$f(x,y)=e^x(4x^2+2y^2+4xy+2y+1)$
jonkin minimin (tai maksimin) etsimiseksi lähtemällä sopivan kaukaa edellä saaduista.

Suorita samoja visualisointeja kuin edellisessä.

\item %3
"Hiihtäjän reittiä" voidaan katsoa myös approksimaationa reitille, jossa kuljetaan joka
hetki jyrkimmän laskeuman suuntaan. Jos reitti on $(x(t),y(t)),$ niin on oltava
$(x'(t),y'(t))=-\nabla f(x(t),y(t)) \ \ \forall t,$ ts. funktiot $x(t)$ ja $y(t)$ ovat
diffyhtälösysteemin
$$\begin{cases}
 x'=f_x(x,y) \\
 y'=f_y(x,y) 
\end{cases}$$
ratkaisuja, ja reitti on systeemin trajektori.\\
Suorita latauskomento {\tt with(DEtools)} ja piirrä {\tt DEplot}-komennolla suuntakenttä
ja trajektoreita, erityisesti nyt aluksi se, joka alkaa harj. 1:n hiihtäjätehtävän
alkupisteestä. Mikä on alkuperäinen hiihtäjän reitti diffyhtälömenetelmänä ajateltuna?

Tee kuva myös Matlab-funktion {\tt pplane5} avulla. (Unix:ssa varo potenssia, 
grafiikkavalikot eivät salli väkästä näppäimistöltä, se on tehtävä leikkaus/liimaus-
tyyliin tai sitten luettava  ".pps-tekstitiedostosta".)

Huom! Matlab-työkalu antaa paljon tarkemman tuloksen, koska Maplen {\tt DEplot} laskee
varsin karkealla, kiinteän askelpituuden numeerisella menetelmällä. Maplessa on paitsi 
hyvin kehittynyt analyyttinen ratkaisija {\tt dsolve}, myös hyvä numeerinen 
{\tt dsolve(...,numeric, method=..)} (useita menetelmiä valittavissa).
(Toki Matlab on tehokkaampi ja monipuolisempi numeerisella puolella.)

\item %4
Opiskele {\tt L/odemaple.mws}-työarkin kohdat 1 ja 2 alakohtineen. \footnote{Viittaamme L:llä kurssin luentohakemistoon
ja H:lla harjoitushakemistoon} (myös odemaple.html olemassa).

Muodosta kohdissa 1.1--1.5 esitetyt $2\times 2$ lineaaristen systeemien faasitasopiirrokset
näillä välineillä.

Tee sitten sama Matlabin {\tt pplane5}:lla. (On suoraan TKK:lla polun varrella, muita ympärostöjä
varten haettavissa {\tt matlab}-hakemistostamme. )


\item %5
Määritä systeemin 
$$
\begin{cases}
x'=x(1-x/2-y)\\
y'=y(x-1-y/2)
\end{cases}
$$

kriittiset pisteet (KRP) ja linearisoi systeemi niissä.
Piirrä faasikuva systeemistä ja erikseen kustakin linearisoidusta.
Voit laittaa joukkoon myös isokliinejä.

Käytä Jacobin matriisin muodostukseen {\tt linalg[jacobian]}-funktiota.

Tässä riittää tehdä faasikuvat joko Maplella tai Matlabin pplane5:llä.


\end{enumerate}


\subsubsection*{Ohjeita}

\textbf{Jyrkimmän laskeuman menetelmä, "method of steepest descent":}

Etsitään kohdefunktion $f(x,y)$ (yleisemmin $f(x_1,\ldots,x_n)$ ) minimiä etenemällä
negatiivisen gradientin suuntaan. Harj. 1 teht. 3:ssa harjoiteltiin "hiihtäjän reittiä",
joka on yksi versio asiasta.
Siinä otettiin kiinteä askelpituus $h$ (prujussa "timestep"). \\
Askelpituuden määrääminen suoritetaan perusSD-menetelmässä niin, että kohdefunktio
minimoidaan etsintäsuunnalla.

Työarkilla {\tt H/harj2ohjeita.mws} on ehdotelmia alkuuunpääsemiseksi.

Määritellään funktio $\phi(t)=f(\V{x}_n+t\V{u}),$ missä $\V{u}$ on 
%$-\nabla f(\V{x_n})$ normeerattuna yksikkövektoriksi, siis kohdefunktio laskettuna 
etsintäsuunnalla. (Nykypisteestä lähtevällä, $\V{u}$:n suuntaisella (puoli)suoralla)
Määrätään $\phi$:n minimi tällä suunnalla, ohjetyöarkin tapauksessa se on helppoa.


%Kts. myös ohjeet: v2


\subsubsection*{Kriittisten pisteiden (gradientin 0-kohtien) luokittelu}

Tässä peruskurssi 2:n kertausta:

Kriittinen piste (KRP) $p$: $\nabla f(p)=0$.\\
Kriittisen pisteen laatu selviää (jos selviää) Hessen matriisin $H_f(p)$ definiittisyydestä seuraavasti:

{\bf Lause} Olkoon $f:\R^n \to \R$ 2 kertaa jatkuvasti derivoituva KRP:ssä $p$.
Olkoon $H=H_f(p)$ Hessen matriisi pisteessä $p$, ja olkoot $\lambda_1,\ldots \lambda_n$ $H$:n ominaisarvot.
 %$q(h)=\frac{1}{2} h^T H_f(p) h$.

(a) Jos $\lambda_i > 0, i=1\ldots n$, niin KRP $p$ on minimi. \\
(b) Jos $\lambda_i < 0, i=1\ldots n$, niin KRP $p$ on maksimi. \\
(c) Jos jokin $\lambda_i >0 $ ja jokin $\lambda_j < 0$, 
niin KRP $p$ on satulapiste (siis ei ääriarvo). \\
(d) Jos kaikki ominaisarvot ovat samanmerkkiset ja jokin on $=0$, niin
"testi pettää".

%\textbf{Lineaaristen faasitasonimityksiä}

%- Lähdenoodi (epästabiili) ``source node'', nielunoodi (vahvasti stabiili)
%\qquad \qquad \quad ``sink node''  \\
%- Satulapiste ``saddle point'' (epästabiili)\\
%- Keskus (stabiili) ``center''\\
%- Lähdespiraalipiste (epästabiili spiraalip.), nieluspiraali (vahvasti stabiili)
 

\subsubsection*{Dys:n linearisointi ja (KRP):t = (TP):t}. 


Diffyhtälösysteemin linearisointi tapahtuu 
laskemalla ensin kriittiset pisteet (KRP), 
eli oikean puolen 0-kohdat.
(Käytetään myös nimitystä tasapainopiste (TP).)
Jos $ {\bf p_0}$ on KRP, niin linearisoidun systeemin (KRP siirretään O:oon) 
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}))_{i,j}$
pisteessä $ {\bf p_0}$.

\end{document}


\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


V3/02/H8
\item %2
Kirjoita vaimenevaa värähdysliikettä kuvaava diff. yhtälö 

$$y'' + 0.125 y' + y = 0$$

ensimmäisen kertaluvun
systeemksi, ja määritä sen yleinen ratkaisu. 
Riittää, kun jätät tuloksen kompleksieksponenttifunktiomuotoon.
(Älä säiky, vaikka saisit tyyppiä $\sqrt{255}$ olevia lukuja.)

Tarkista tulos alkuperäisen yhtälön avulla (myös kompleksimuodossa).


\item %4
Olkoon
$A=
\left[ \begin {array}{cc} 2&3\\\noalign{\medskip}-1&-2\end {array} \right]
$
Ratkaise alkuarvotehtävä ("AA-tehtävä") 
$\V{y}' = A \V{y}$, kun $\V{y}(0)=(3,2)$.

Selvitä kriittisen pisteen $\V{0}$ luonne (lähde, nielu, satula tms.) ja piirrä
kuvaan ominaisvektorit, edellä laskettu ratkaisutrajektori, sekä muutama muu 
trajektori suuntanuolineen (käsin hahmotellen).

\item %5
 

Olkoon
$A=
\left[ \begin {array}{cc} 7&-1\\\noalign{\medskip}3&3\end {array} \right] 
$
Ratkaise taas AA-tehtävä 
$\V{y}' = A \V{y}$ samalla alkuehdolla $\V{y}(0)=(3,2)$.

Tee vastaavat asiat kuin edellisessä tehtävässä.



\item %6
Osoita, että alla ohjeessa sanottu
matriisieksponenttifunktion sarjan suppenee kaikilla neliömatriiseilla $A$, 

Tarvitset kahta asiaa:

1) $e^x=1+x+\frac{x^2}{2!}+ \ldots + \frac{x^n}{n!} + \ldots$
suppenee kaikilla $x \in \R$.

2) Yläraja-arviota matriisin $A^k$ alkioille,jotta pääset käyttämään edellä
olevaa sarjaa vertailusarjana kullekin alkiosarjalle.  

Voit rajoittua matriisiin, jonka alkiot ovat samoja, koska sillekin
asian tulee päteä, ja toisaalta mielivaltaisen matriisin tapauksessa
saadaan yläraja-arvio korvaamalla kaikki matriisin alkiot itseisrvoltaan
suurimmalla. 
Jos tuntuu mukavalta, voit hyvin rajoittua $2 \times 2-$ matriisiin.
Kaikki olennaiset seikat saat näin ollen esille tarkastelemalla matriisia

$A=\left[ \begin {array}{cc} 1&1\\\noalign{\medskip}1&1\end {array} \right]$.

\end{enumerate}





\subsubsection*{Loppuviikko (LV)}


%\subsubsection*{Ohjeita}

Katso nyt ainakin tätä: {\tt L/L12maple.mws ja .html}


\begin{enumerate}


\item % 1
Piirrä tehtävien 4 ja 5 faasitasokuvat. Piirrä myös ominaisvektorit {\tt arrow}-komennolla ja 
suuntakenttä {\tt DEplot}:lla. (Malli: {\tt L/L12maple.mws}.)
Piirrä lisäksi kysyt alkuarvoratkaisut aikatasossa.


\item %2
Ratkaise systeemi $\V{y}'=A \V{x}$, missä

$A=\left [\begin {array}{cc}           
                            2&10\\      
                            -1&4    
          \end {array}\right ] $



alkuehdolla $\V{y(0)}=[1,0]^T$

Piirrä faasitasokäyriä, suuntakenttä ja selvitä KRP:n $0$ luonne.

%\kohta{b}{alkuehdolla $\V{x(0)}=[0,1]'$}
%\kohta{c}{Piirrä pplane5-piirroksia, katsele niin faasitaso
%\footnote{Faasitasokäyrät ovat muotoa $(y_1(t),y_2(t)),$
%kun $t$ on käyräparametrina.} 
%kuin aikakuvaa
%(time series)}



\item %3
Tarkastellaan matriisin  $A=
\left [\begin {array}{cc} \alpha&2\\\noalign{\medskip}-2&0\end {array}
\right ]$
määräämää (HY):öä $y' = A y$.

Kuvaile, miten ratkaisuparven kvalitatiivinen luonne muuttuu 
$\alpha$:n muuttuessa.
Selvitä trajektorien käyttäytyminen sekä origon luonne ja stabiilisuus
eri tapauksissa, ja havainnollista piirroksin (käsin hahmotellen).

Piirrä faasipiirroksia kolmeen luonteeltaan erilaiseen tapaukseen
liittyen (Maplella). Mainitse KRP:n $O$ luonne ja stabiilisuus.

Ota mallia {\tt L/harj12maple.mws-työarkista nyt erityisesti.}


\item %4
Laske kynää ja paperia ja tarvittaessa esim.
Maplea matriisilaskimena käyttäen:
$e^{A t}$ matriiseille \\ 
a) $A=\left [\begin {array}{cc} 1&0\\
\noalign{\medskip}0&2\\\noalign{\medskip} \end{array}\right]$
b) $A=\left [\begin {array}{cc} 0&1\\\noalign{\medskip}0&0\\
\noalign{\medskip} \end{array}\right]$
c) $A=\left [\begin {array}{cc} 0&1\\\noalign{\medskip}-1&0\\
\noalign{\medskip} \end{array}\right]$

Kommentteja: a) no comments, b) on ``nilpotentti'', korota
potenssiin, niin näet, mitä tarkoittaa, c):ssä näkyy tietty jaksollisuus.

Voit tarkistaa {\tt linalg[exponential]}-komennolla.

Muotoile tehtävät, joiden ratkaisut voit näiden avulla kirjoittaa.


%Vast:
%a) $e^{A t}=\left [\begin {array}{cc} e^t&0\\
%\noalign{\medskip}0&e^{2t}\\\noalign{\medskip} \end{array}\right]$
%b)$e^{A t}= \left [\begin {array}{cc} 1&t\\\noalign{\medskip}0&1\\
%\noalign{\medskip} \end{array}\right]$
%c) $e^{A t}=\left [\begin {array}{cc} \cos t&\sin t\\
%\noalign{\medskip}-\sin t&\cos t\\
%\noalign{\medskip} \end{array}\right]$

\item
Epähomogeenisen systeemin
$$ \V{x}' = A \V{x} + \V{g}(t), \ \ \V{x}(0)=\V{x}^{(0)}$$
ratkaisu on
$$\V{x}(t)=e^{At}\V{x}^{(0)} + \int_0^t e^{A(t-s)}\V{g}(s) ds .$$

Vektorifunktion integrointi tarkoittaa yksinkertaisesti integrointia
komponenteittain.

Ratkaise t\"am\"an avulla systeemi

$\V{x}'=\left [\begin {array}{cc} 1&4\\\noalign{\medskip}-1&-3\end {array}
\right ] \V{x} + \left [\begin {array}{c} t\\\noalign{\medskip}1\end {array}
\right ],\ \  \V{x}(0)=(1,-1) $

Voit katsoa mallia tästä: \\
\verb_ www.math.hut.fi/teaching/k3/luentomateriaali/L20maple.html _

Mallista poiketen saat tässä oikaista niin, että käytät valmista
{\tt linalg[exponential]}-funktiota (tätäpä ei näy uudesta 
{\tt LinearAlgebra})-kirjastosta löytyvän.)

\item %6
Määritä yleistettyjä ominaisvektoreita käyttäen $e^{A t}$, kun \\ 
$A=
\left[ \begin {array}{cccc} 1&0&1&2\\\noalign{\medskip}1&1&2&1\\
\noalign{\medskip}0&0&2&0\\\noalign{\medskip}0&0&1&1\end {array}
 \right].
$
\\Ratkaise AA-tehtävä 
$\V{y}'=A \V{y}$
alkuehdolla $\V{y(0)}=[1,0,1,0]^T$.



\end{enumerate}


%\end{document}



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

Hyvä tiivis esitys aiheesta on Pekka Alestalon kirjoittamassa tekstissä.
Siihen on viite www-sivuillamme. 

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. TE-prujusta. 

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 (viimeistään) 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: 


\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 \matl:ssa {\tt expm}. 

(Toki myös yleiset numeeriset ratkaisumenetelmät toimivat yleensä hyvin linearisiin 
sovellettuna.)


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


\end{document}

\begin{enumerate}

\item
AV-tehtävä 6 Maple-avusteisesti havainnollisine kuvineen.

\item
Johda diffyhtälösysteemi oheisen kuvan 
%~\ref{jouset} 
massa-jousisysteemille.
Voit olettaa: $x_1 > x_2 > x_3.$

\begin{figure}
  \begin{center}
    \leavevmode
     \epsfig{file=kuvat/harj6k2.eps, height=4cm, width=7.5cm}
     \caption{}
  \end{center}
\end{figure}


Ratkaise tehtävä Laplace-muunnostekniikalla tietenkin.

Massat ja jousivakio(t): $m_1=1, m_2=2, m_3=2, k=1$ \\
Alkuehdot:$x_1(0)=1, x_2(0)=2, x_3(0)=3, x_1'(0)=x_2'(0)=x_3'(0)=0 .$ \\
Ulkoinen voima (heräte): $f(t)=\cos \omega t$, olkoon $\omega=1 .$

Vihje: Tätä ei tarvitse kirjoittaa 1. kl. $6\times6$-systeemiksi vaan
voidaan hyvin käsitellä 2. kl. $3\times 3$-systeeminä.
Kirjoita kukin yhtälö erikseen tyyliin
\begin{verbatim}
> dyht1:=diff(y1(t),t,t)=...);dyht2:=...;dyht3:=...;
\end{verbatim}
L-muunna ja ratkaise.

\item
Määritä diffyhtälörymän
$$\begin{cases}
y_1'=y_1+y_2 \\
y_2'=5y_1-3y_2
\end{cases}$$
yleinen ratkaisu.

Käytä matriisitekniikkaa (ominaisarvoja). Määritä (AA)-ratkaisu $y_1(0)=0,
y_2(0)=1$.


\item
Muunna seuraavat diffyhtälö(ryhmät)t 1. kertaluvun systeemeiksi:
(Ei tarvitse ratkaista.)

(a) $y''-3y'+y=0$, \quad (b) $y'''+e^t y' + y=0$, \quad 
(c) $\Theta''+k\sin \Theta=0$, \quad 
(d) $\begin{cases}
     y_1''-y_1'-2y_1=t^2 \\
     y_2''-y_2-3y_1=0 
    \end{cases}$

Huomaa, että (c) on epälineaarinen (heiluriyhtälö), sille on turha etsiä
matriisia.

\item
Ratkaise systeemi $\V{x}'=A \V{x}$, missä

$A=\left [\begin {array}{cc}           
                            2&10\\      
                            -1&4    
          \end {array}\right ] $
\kohta{a}{alkuehdolla $\V{x(0)}=[1,0]'$}
\kohta{b}{alkuehdolla $\V{x(0)}=[0,1]'$}
\kohta{c}{Piirrä pplane5-piirroksia, katsele niin faasitaso
%\footnote{Faasitasokäyrät ovat muotoa $(y_1(t),y_2(t)),$
%kun $t$ on käyräparametrina.} 
kuin aikakuvaa
(time series)}


\item
Kaksi s\"aili\"ot\"a, $A$ ja $B$ sis\"alt\"a\"a $50 l$ nestett\"a kumpikin. 
Niit\"a yhdist\"a\"a
kaksi putkea siten, ett\"a $A \to B$ virtaa nestett\"a 
nopeudella $4$ l/min ja
$B \to A$ $1$ l/min ja oletetetaan, ett\"a neste sekoittuu heti. 
Puhdasta vett\"a virtaa s\"aili\"o\"on $A$ nopeudella $3$ l/min
ja s\"aili\"ost\"a $B$ poistuu nestett\"a niinik\"a\"an nopeudella $3$ l/min.
Oletetetaan, ett\"a alkuhe
tkell\"a s\"aili\"o $A$ 
sis\"alt\"a\"a $25$ kg suolaa ja
s\"aili\"o $B$ pelkk\"a\"a vett\"a. M\"a\"arit\"a 
suolam\"a\"ar\"at kummassakin s\"aili\"oss\"a
ajan funktiona. 

%\vfill\eject

\begin{verbatim}

3 l/min 
__________          __________
-->       |         |         |
__        |_________|         |
 |         4 l/min ->         |
 |        -----------         |
 |        |         |         |
 |    A   |         |    B    |
 |        |---------|         |
 |         <-- 1 l/min        |
 |        |---------|         ----
 |________|         | ____________-> 3 l/min

\end{verbatim}

Miten pitkän ajan kuluttua suolamäärät säiliöissä ovat yhtäsuuret?
Piirrä $x_1(t)$ ja $x_2(t)$ aikariippuvuuskuvana (``time series'')
ja faasitasokäyränä eli $(x_1,x_2)$-tasossa, pitämällä aikaa $t$ 
käyräparametrina . (Ilman pplane5:n apua, muuten saat käyttää Matlabia
tai Maplea.)

\item
Muodosta yleinen ratkaisu yhtälölle $x' = A x $, kun
$
\left [\begin {array}{cc} 8&-14\\\noalign{\medskip}3&-4\end {array}
\right ]
$
Piirrä trajektoreita ja suuntakenttä pplane5:n avulla.

\item
Osoita, että 
matriisieksponenttifunktion määritelmäsarja

$$e^A=I +A+\frac{A^2}{2!}+ \ldots + \frac{A^n}{n!} + \ldots$$

suppenee mielivaltaisella neliömatriisilla $A$.

Todistus perustuu kahteen asiaan:

1) $e^x=1+x+\frac{x^2}{2!}+ \ldots + \frac{x^n}{n!} + \ldots$
suppenee kaikilla $x \in \R$.

2) $||A^k||$ ei kasva liian lujaa.

Jälkimmäiseen on kätevää käyttää matriisinormiepäyhtälöä: 
$||A B|| \leq || A || || B || . $

\end{enumerate}



\subsubsection*{Loppuviikko (LV)}


\begin{enumerate}


\item
Laske kynää ja paperia ja tarvittaessa esim.
Maplea matriisilaskimena käyttäen:
$e^{A t}$ matriiseille \\ 
a) $A=\left [\begin {array}{cc} 1&0\\
\noalign{\medskip}0&2\\\noalign{\medskip} \end{array}\right]$
b) $A=\left [\begin {array}{cc} 0&1\\\noalign{\medskip}0&0\\
\noalign{\medskip} \end{array}\right]$
c) $A=\left [\begin {array}{cc} 0&1\\\noalign{\medskip}-1&0\\
\noalign{\medskip} \end{array}\right]$

Kirjoita sitten vastaavan (HY):n yleinen ratkaisu.

Kommentteja: a) on ääliömäisen helppo, b) on ``nilpotentti'', korota
potenssiin, niin näet, mitä tarkoittaa, c):ssä näkyy tietty jaksollisuus.

%Vast:
%a) $e^{A t}=\left [\begin {array}{cc} e^t&0\\
%\noalign{\medskip}0&e^{2t}\\\noalign{\medskip} \end{array}\right]$
%b)$e^{A t}= \left [\begin {array}{cc} 1&t\\\noalign{\medskip}0&1\\
%\noalign{\medskip} \end{array}\right]$
%c) $e^{A t}=\left [\begin {array}{cc} \cos t&\sin t\\
%\noalign{\medskip}-\sin t&\cos t\\
%\noalign{\medskip} \end{array}\right]$


%Tee ohjetyöarkin mukaan. Katso myös ``.../luentomateriaali/L21.html'' 
%(teoria) ja erityisesti
%``L21maple.html'' sekä `L21maple.mws''. 

Tarkista Maplella ja kokeile lisäksi tyyliä

\verb_ seq(evalm((t*A)&^k)/k!,k=0..10); map(series,eat,t=0,10); _

olevia temppuja joidenkin matriisien kohdalla.
(Huomaat, että matriisieksponenttifunktio ei ole ollenkaan kummallinen olio,
ehkä siihen kehittyy riippuvuus!)



%URL:ssä
%\verb_ http://www.math.hut.fi/~apiola/maple/opas/ohjelmat.mpl _
%on kerättynä kurssilla tarvittuja ohjelmia. Listaus on HAM-kirjan lopussa,
%mutta tiedostoa on päivitetty (mm. 3.11.) 
%Kyseessä on tekstitiedosto,
%jonka voi ottaa ensin omaan hakemistoon ja lukea \maple-istuntoon komennolla

%\verb_ read("ohjelmat.mpl");_

%Siitä voi vaihtoehtoisesti leikata/liimata haluamiaan koodeja.

%HAM  ss. 161--176 \footnote{[HAM] {\it Heikki Apiola}, 
%Symbolista ja numeerista matematiikkaa \maple-ohjelmalla, ilmestyi 20.10. 98}

%Ohjetyöarkki joko www:stä tai
%{\tt cp /p/edu/mat-1.443/H98/harj8pohja.mws .}

\item
Ratkaise $\V{x}'=A \V{x}$,  $\V{x}(0)=[0, 1, 1]$ kun 
\[
A=\left [\begin {array}{ccc}  -14&-160&-40\\            
          \noalign{\medskip} 181&5&2\\                  
           \noalign{\medskip} 96&84&18\end {array}\right ]
\]
%Ominaisarvot ovat $\lambda_{1,2}=6\pm 180 i$ ja $\lambda_3=-3$ sekä
%tätä vastaava ominaisvektori $w_3=(0,-1,4)$.
Käytä matlab-hakemiston {\tt linsys}-funktiota (jonka toiminta perustuu
Matlabin sisäänrakennettuun $e^{At}$-funktioon nimeltään {\tt expm}, kts. myös
{\tt ...L/L\_ldys.html}).

Piirrä ratkaisukäyräparvi (aikariippuvuuskuva) sekä faasiavaruuskuva.
Piirrä samaan kuvaan reaalinen ominaisvektori ja kompleksisen 
ominaisvektorin re- ja im- osien muodostamat vektorit. Matlabissa
{\em real, imag }. Näiden ``ominaisakselien'' piirto voisi sujua
seuraavaan tyyliin:
\begin{verbatim}
sopivakerroin=3 ; % tms. 
aks=sopivakerroin*v1;plot3([0,aks(1)],[0,aks(2)],[0,aks(3)],'r')
...
\end{verbatim} 

Kokeile myös sellaisilla alkuarvoilla, jotka sijaitsevat kahden ominaisvektorin
(Re/Im-osien) määräämillä tasoilla.

%Voit vielä huviksesi 
%piirtää faasikuvan projektiot koordinaattitasoille. Kts. linsysys-
%koodin helppitekstejä.

%Määritä $\lambda_1$:tä vastaava ominaisvektori $\bar w$ ja lausu yleinen
%reaalinen ratkaisu muodossa $\bar y=a  z_1(t)\bar u + b  z_2(t)\bar v
%+ c z_3(t)\bar w_3$.
%Mieti trajektoreita 3d-avaruudessa ja varustautu piirtoon LV-harj. mennessä.

\item
Johdantoesimerkissä on massa-jousisysteemi, kalvokopiot jaettiin, kuva on
myös KRE s. 158, {\em Fig. 77}. Olkoot $m_1=m_2=1, k_1=3, k_2=2. $

Muunna systeemi 1. kertaluvun ($4\times 4$)-systeemiksi ja ratkaise.
Saat käyttää ratkaisussa Maplea tai Matlabia, etenkin piirtämiseen 
voisi olla suositeltavaa käyttää
jälleen Matlab-funktiomme {\tt linsys}. Piirrä trajektoreita 
$(y_1(t),y_2(t))$. Miten näet trajektorista, milloin massat kulkevat samaan
ja milloin eri suuntaan.

Toisaalta analyyttisen ratkaisulausekkeen saat Maplella mukavasti. Sovitaan,
että saat käyttää {\tt linalg[exponential]}- funktiota.


\item
Luennolla osoitettiin, ett\"a ep\"ahomogeenisen systeemin
$$ \V{x}' = A \V{x} + \V{g}(t), \ \ \V{x}(0)=\V{x}^{(0)}$$
ratkaisu on
$$\V{x}(t)=e^{At}\V{x}^{(0)} + \int_0^t e^{A(t-s)}\V{g}(s) ds .$$

Vektorifunktion integrointi tarkoittaa yksinkertaisesti integrointia
komponenteittain.

Ratkaise t\"am\"an avulla systeemi

$\V{x}'=\left [\begin {array}{cc} 1&4\\\noalign{\medskip}-1&-3\end {array}
\right ] \V{x} + \left [\begin {array}{c} t\\\noalign{\medskip}1\end {array}
\right ],\ \  \V{x}(0)=(1,-1) $

Voit katsoa mallia tästä: \\
\verb_ www.math.hut.fi/teaching/k3/luentomateriaali/L20maple.html _

Mallista poiketen saat tässä oikaista niin, että käytät valmista
{\tt linalg[exponential]}-funktiota (tätäpä ei näy uudesta 
{\tt LinearAlgebra})-kirjastosta löytyvän.)

{\bf Eiköhän tässä ole jo tarpeeksi.}


\end{enumerate}



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


\textbf{Lineaaristen faasitasonimityksiä}

- Lähdenoodi (epästabiili) ``source node'', nielunoodi (vahvasti stabiili)
\qquad \qquad \quad ``sink node''  \\
- Satulapiste ``saddle point'' (epästabiili)\\
- Keskus (stabiili) ``center''\\
- Lähdespiraalipiste (epästabiili spiraalip.), nieluspiraali (vahvasti stabiili)
 

{\bf Linearisointi ja (KRP):t = (TP):t}. Linearisointi tapahtuu 
laskemalla ensin kriittiset pisteet (KRP), eli oikean puolen 0-kohdat.
(Käytetään myös nimitystä tasapainopiste (TP).)
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}))_{i,j}$
pisteessä $ {\bf p_0}$.


%%%%%%%%%%%%%
\documentclass[finnish,11pt,twocolumn]{article}
\usepackage{/home/apiola/tex/KP3harj}
\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\Ln#1{\textit{Ln #1}}
\def\Log#1{\textit{Log #1}}

\begin{document}

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







