\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{6}{13--14}{23 -- 30.3 2004}

P�ivitetty 26.3.04 (BeulerFixP korjattu)

{\bf Huom:} Teht�v��n 5 lis�iltiin hiukan.

N�m� teht�v�t k�yd��n l�pi {\bf ti 30.3.}

Annetaan viel� teht�v�t 7 ennen p��si�ist�, sitten onkin loppuprojektien vuoro.

Sopivia harj7-teht�vi� ovat mm. 
Moler: Ch 11: PDE: teht. 11.1, 11.2, 11.6
%{\tt http://www.math.hut.fi/teaching/numsym/04/matlab/moler/odes.pdf} \\

%{\tt http://www.math.hut.fi/teaching/numsym/04/matlab/CV/CHAPTER.9} \\

%{\tt http://www.math.hut.fi/teaching/numsym/04/L/ODE.html} ODE-"portaali" \\
{\tt http://www.math.hut.fi/teaching/numsym/04/maple/ns04.mpl} Maple-kokoelma  \\
{\tt http://www.math.hut.fi/teaching/numsym/04/matlab/dfield5.m} LAODE-funktio
(muistathan viel� {\tt pplane5}:n) \\



%Sek� [CvL]:ss� ett� Moler:ssa on laaja teht�v�valikoima.
%Erityisesti Molerista on mahdollisuus valita my�s laajempia oppilaskohtaisia teht�vi�, joissa
%p��see mukavasti k�ytt�m��n eri alueiden laskentav�lineit� yhteisty�ss�.

%Olkoot lyhyesti C ja M nuo viitteet alla olevien teht�v�numeroiden kohdalla.



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



%\subsubsection*{Aikataulu}


%Ti 16.3: Neuvontaharjoitus harj. 4 ja harj5 teht�vien
%edist�miseen. 
 
%To 18.3: K�yd��n l�pi teht�v�t 4 ja 5. 
%T�ydennet��n teht�vi� 5: ``stability, stiffness'', annetaan my�s teht�vi�
%osdy-aiheista. Alla on jo osittain t�ydennetty, mutta torstaiksi riitt��
%tehd� aiemman, puhtaasti Matlab-tekniikkaa opettavan tavan mukaan.



\begin{enumerate}

\item %1

Aiheena kankeat diffyht�l�t. Jatketaan harj.5 teht. 3:a.
Yht�l�h�n on
$$y'=-1000(y-\sin t) + \cos t .$$

Jotta voitaisiin paremmin havainnollistaa kuvilla, lievennet��n kankeutta
ottamalla vakion $1000$ tilalle vaikkapa $50 .$

Muodosta Maplella yleinen ratkaisu antamalla alkuehto muodossa $y(a)=ya$.

Piirr� ratkaisufunktioita vaikka ohjeen mukaan.

M��rittele sopiva grafiikkafunktio ohjeen tyylisesti ja piirtele sill� joitakin
 kuvia 
 ja rakenna sen avulla
 ratkaisuk�yr�parvi. (Piirr� vaikka alapuoli sinisell� ja yl�puoli punaisella.)
 
Tee suunakentt�- ja ratkaisuparvipiirros Matlab-(LAODE)-funktiolla {\tt dfield5}.

\item %2

(a)
Ratkaise edellinen yht�l� Euerin menetelm�ll� v�lill� $[0,\pi]$.
Tiedostossa {\tt .../maple/ns04.mpl} on {\tt Eulerh}, listaus {\tt harj5}-paperissa.
Kokeile erityisesti askelpituuksia $h_{\textit{max}}$, $h_{\textit{max}}/2$ ja 
esim. $1.2 h{\textit{max}}$.


(b) 
Tutustu {\tt ../maple/ns04.mpl}-tiedoston funktioon {\tt BeulerFixP}, selvit� itsellesi
sen toiminta ja sovella sit� samaan teht�v��n.
Suppenemista voit tarkkailla vaikka tyyliin:\\
\verb_hmax:=...; BEF:=BeulerFixP(F,[0,Pi],1,0.9*hmax,5): %[-1];_
Kun saat j�rkev��, voit piirt��: \verb_plot(BEF);_

\textbf{Varoitus!} Edell� olevan tyylisiss� kutsuissa on syyt� varmistaa, ett�
kaikilla parametreilla on j�rkev�t arvot. Jos vaikka F:n m��rittely unohtuu, niin
Maple ryhtyy ikuiseen jauhantaan, 
jolloin ollaan t�ysin {\bf stop}-nappulan armoilla.

(c) 
Selit� teht�vien 1) ja 2) kuvien perusteella, mist� t�ss� on kyse!

\item %3
Kirjoita Maple-funktio {\tt BeulerNewt}, johon on hyv� pohja ja
ohjeet olemassa.

Tee ensin versio, jossa on {\tt niter} parametrin� aivan kuten "Fix":ss�.

Voit ter�st�� funktiosta ``aidon implisiittimenetelm�n'' iteroimalla
toleranssiin saakka. Lopetusehtona vosi olla esim.

$$|y_{n+1}^{(k+1)}-y_{n+1}^{(k)}|\leq 10^{-6}|y_{n+1}^{(k+1)}|. $$

(jos suht. toleranssina on $10^{-6}$.) Olkoon ter�stys "vapaaehtoista", "niter-versiokin"
toimii varsin hyvin.

Testaa t�t� menetelm�� eri askelpituuksilla.

Kirjoita vaikka grafiikkafunktio:

\verb_BenKuva:=h->plot(BeulerNewt(F,[0,Pi],1,h,5)) _
(jos k�yt�t ``niter-versiota'')




Selvit� taas. mist� on kyse.

\item %4
Tarkastele lineaarista yht�l�systeemi� $y' = A y$, miss�
$A= \left[ \begin {array}{cc} 998&1998\\\noalign{\medskip}-999&-1999
\end {array} \right] 
$

Alkuarvoteht�v�n ratkaisu voidaan kirjoittaa muotoon $y(t)=e^{A t}y(0),$ joka
Maplella voidaan tehd� alla kerrattuun tyyliin.

T�ss� on havainnollisempaa laskea ratkaisu ominaisarvojen avulla ominaisvektorikannassa:

$y(t)=C_1 e^{\lambda_1 t}v_1 + C_2 e^{\lambda_2 t}v_2.$ 

Selvit� t�m�n muodon perusteella, mist� systeemin kankeus johtuu. 

Laske ja piirr� (faasikuvana) jokin Euler-askel AA-teht�v�lle, $y(0)=[1,0]^T$.
T�m�n voit tehd� joko Maplella tai Matlabilla. Maple-ratkaisuun voit etsi� valmiin
{\tt EulerV}:n {\tt ..maple/ns04.mpl:}st�. (Ei ole pahitteeksi, vaikka kirjoittaisit
koodin yhdenmukaisesti edell� esiintyneiden Euler-versioiden kanssa.)
Mielell��n voit piirt�� my�s faasikuvan {\tt DEplot}:lla \\

Matlab-ratkaisussa voit piirt�� faasitasokuvan {\tt pplane5}:ll� ja samaan
kuvaan joitakin Eulerin polkuja. ( Kokeile {\tt figure(1); figure(2) ...} sitten
{\tt hold on} ja piirtoa rohkeasti.\\
Kuvissa voit taas skaalata siten, ett� arvon $1000$ sijasta k�yt�t vaikka arvoa $50$.

Jos kokeilet arvoa $1000$ pplane5-piirroksessa, niin kankeuden pit�isi n�ky� hitautena
ratkaisuk�yrien piirrossa, koska pplane5 k�ytt�� ei-kankeaa ratkaisijaa ode45.



\item
Tutki kolmen aineen kemiallisia reaktioita kuvaavaa yht�l�systeemi�:

$$\begin{cases}
y_1'=-0.04 y_1 + 10^4 y_2 y_3 \\
y_2'=0.04y_1 -10^4 y_2y_3-3\times 10^7y_2\\
y_3'=3\times 10^7y_2^2
\end{cases}
$$

Muodosta Maplella systeemin Jakobiaani y-muuttujien suhteen ja m��rit� sen ominaisarvot.
P��ttele niiden perusteella systeemi kankeaksi (ainakin alla olevassa alkuarvopisteess�). 

Ratkaise AA-teht�v�, jossa alkuarvona on vaikkapa $[1,0,0]^T$, v�lill� $[0,3]$.
Vertaa ratkaisijoita {\tt ode45} ja {\tt ode15s}.
Tutki askelten lukum��ri�, havainnollista piirroksin, voit k�ytt�� my�s optiota:
{\tt odeset('Stats','on')}.

(Jos tutkiskelet jaettuja materiaaleja, voit tehd� hyvi� l�yt�j�.)

\textbf{Huom!} ke 24.3. tein pieni� lis�yksi� t�h�n teht�v��n:
Ota Maplella generoitu jakobiaani Matlab-editoriin ja tee siit� Matlab-funktio vaikkapa
"jakobiaani". Suorita {\tt optiot=odeset('Jacobian',@jakobiaani)}.

Jakobiaania ei tarvita ode45:ss�, mutta ei se haittakaan. Sensijaan saattaa k�yd� niin,
ett� ode45 laskee toivottoman kauan. Jotta lasekennasta saataisiin jotain tietoa
(ennen {\tt CTR-C}:t�),
kannattaa laittaa globaali muuttuja, vaikkapa {\tt kutsut} diffyht�l�systeemiin
laskuriksi.
\begin{verbatim}
function yp=kemikaalit(t,y)
global kutsut
kutsut=kutsut+1
...
\end{verbatim}
Istunnossa (p��ohjelmassa) ennen ode45-kutsua kirjoitetaan {\tt global kutsut}
ja {\tt kutsut=0}.



\end{enumerate}

\subsection*{Ohjeita}

Yleisesti ottaen kannattaa palautella mieliin harj2:n ja harj3:n teht�vi� ja ohjeita.

Teht. 1: Diffyht�l�n analyyttinen ratkaisuhan saadaan tyyliin \\
\verb_ dsolve({dyht,y(0)=ya},y(t));_ \\
Ratkaisulausekkeen poimiminen on kaikkein yksinkertaisinta n�in: \\
\verb_ Y:=rhs(%); _ Funktioksi muuttaminen: \verb_Yf:=unapply(Y,t);_

"Mielenkiintoinen" lis�piirre tulee t�ss�, kun Maple haluaa lausua tuloksen $\cosh:$n
ja $\sinh:$n avulla. Nimitt�j�n $\cosh(kt) - \sinh(kt)$ aiheuttaa nollalla jaon 
pienill�kin $t$:n arvoilla ainakin, jos $k=1000$.

Kannattaa leikata/liimata nimitt�j� ja sanoa sille: {\tt nim:=convert(nim,exp);}.
Outputista saa helposti leikkaus/liimaus-tekniikalla input-muodon, jota sitten
sopivasti editoidaan. (Varmasti helpompaa kuin sopivan {\tt subs}-komennon miettiminen.)
Pieni nimitt�j� poistuu (itse asiassa nimitt�j��n ei j�� mit��n).
 
Kankeusasian ymm�rt�miseksi kannattaa ensin piirt�� ratkaisuk�yr�parvi analyyttist� ratkaisua
k�ytt�en. T�h�n tyyliin voisit edet�:

\begin{verbatim}
 plot(subs(a=0.10,ya=1,Y),t=0..0.5,y=0..2);
 plot(subs(a=0.10,ya=1,Y),t=0..Pi,y=-1.5..1.5);
\end{verbatim}

Muistellaan tehokasta Maple-tekniikkaa sopivia grafiikka-arvoisia 
(``kertak�ytt�'')funktioita m��ritellen. (Funktion m��rittely�
editoidaan aina kulloisenkin modifikaatiotarpeen mukaan.)

\begin{verbatim}
 kuvaf:=(A,Tvali,Ya,vari)->plot(subs(a=A,ya=Ya,Y),t=Tvali[1]..Tvali[2],
y=-1.5..1.5,color=vari);

display(seq(kuvaf(A,[0,Pi],1,red),A=[seq(k*0.1,k=0..30)]));
... blue ...
display(%,%%);
\end{verbatim}

Kuvajonosta saadaan animaatio lis��m�ll� display-funktion valitsimeksi {\tt insequence=true}.

\subsection*{Implisiittinen Euler}

$y_{n+1}=y_n+h f(t_{n+1},y_{n+1})$.

Ratkaistava $y=y_{n+1}$ yht�l�st� $F(y)=y$, miss�
$F(y)=y_n+h f(t_{n+1},y).$

Teht�v� on suoraan kiintopisteiteraatiomuotoa (kts. [HAM] ss. 66--68.)

Menetelm�n Maple-implementointi on nautinnollista, m��ritell��n
lokaali iterointifunktio $F$. Kiintopisteiteraation suppenemisehto
on $|F'(x)| < 1$.

T�ss� ote tiedostosta {\tt .../maple/ns04.mpl} (p�ivitetty 22.3.04)
Kommenteissa on selitys menetelm�n rajoituksesta, joka johtuu yll�
mainitusta suppenemisehdosta.


\hrule
{\it Backward Euler, eli implisiittinen Euler. T�ss� iteroidaan kiintopisteiteraatiolla.
 Suppeneminen vaatii askelpituusrajoituksen $|h f_y| < 1$
   K�yt�nn�ss� aika kelvoton, eik� omaa implisiittisen menetelm�n p��etua:
   askelpituusrajoituksesta vapautumista.
    Helppo impelemntoida.}
\hrule
\begin{verbatim}
 BeulerFixP:=proc(f,vali,y0,h,niter)
local a,b,T,Y,F,s,n,N;
F:=unapply(Y[n]+h*f(T[n+1],s),s);
a:=evalf(vali[1]); b:=evalf(vali[2]);
N:=ceil((b-a)/h);
Y[0]:=evalf(y0);
T[0]:=a;
for n from 0 to N do
  T[n+1]:=T[n]+h;
  Y[n+1]:= (F@@niter)(Y[n]+h*f(T[n],Y[n]));
end do;
[seq([T[n],Y[n]],n=0..N)];
end:
\end{verbatim}

Implisiittisest� Eulerista on juttua [HAM] ss. 123 -- 124.

\textbf{Iterointi Newtonin menetelm�ll�}

Newtonin menetelm� ei k�rsi yll� olevasta rajoituksesta, ainoa
pieni ep�m��r�isyys on, ett� alkuarvauksen on oltava 
``riitt�v�n hyv�'', jotta menetelm� suppenee, 
jolloin se yleens� suppenee nopeasti.
``Oikea'' implisiittimenetelm� saadaan iteroimalla siihen saakka,
kunnes toleranssiraja saavutetaan.

4)
Newtonia k�ytt�v� implisiittinen Euler saadaan kiintopisteEulerista
helposti.
[HAM]-versiossa ({\tt BEN}) (joka on my�s {\tt ns04.mpl}:ss�), on turhaan viety 
paikallisten funktioiden m��rittelyt
for-silmukan sis��n. Tee yll� olevan {\tt BeulerFixP}:n 
henkinen versio.

Systeemin $y'=A y$ ratkaisu: $y=e^{A t}y_0$:
\begin{verbatim}
> with(LinearAlgebra):with(linalg):
> A:=<<a11 | a12>,<a21 | a22>>;
> eAt:=Matrix(exponential(A,t));
\end{verbatim}
Huomaa, ett� {\tt exponential} on {\tt linalg-} funktio, jos/kun haluamme k�sitell� tulosta {\tt LinearAlgebra}-
rakenteena, t�ytyy soveltaa {\tt Matrix}:ia. (��neen lukien tuo kuullostaa juuri silt�, kuin pit��kin.)

T�ss� on erityisen hy�dyllist� k�ytt�� {\tt Eigenvectors}-
komentoa, jolloin saadaan johtop��t�ksiin parhaiten sopiva
muoto (esitys ominaisvektorikannassa, kun sellainen on 
tarjolla).

5)
Jakobiaanin muodostamiseen tarvitaan vanhaa linalg:ia, siis {\tt with(linalg):}
{\tt ?jacobian}.
{\tt lprint} antaa Matlabiin siirrett�v�n muodon. K�tev�� on ensin muuntaa matriisi
listojen listaksi: {\tt convert(J,listlist);}\\
Jos leikkau/liimaus ei onnistu, voit vaikka komentaa:
\begin{verbatim}
> writeto("jac.txt"):
> lprint(...);
> writeto(terminal):
\end{verbatim}

\begin{verbatim}


\end{verbatim}



%Tee suunakentt�- ja ratkaisuparvipiirros Matlab-(LAODE)-funktiolla {\tt dfield5}.





\end{document}

Seuraavana asiakokonaisuutena k�sitell��n osittaisdiffyht�l�iden numeriikkaa.
Sopivana siirtym�vaiheena otetaan ``method of lines'', jossa ``laskentakoneena'' on ODE-ratkaisija (suuri tavallinen differentiaaliyht�l�ryhm�).

[Moler]-prujuissa on koko joukko kiintoisia sovellusteht�vi�, kuten t�st�
ODE-osasta huomataan. 
N�ihin liittyen on my�s mahdollista kehitell� loppuprojektiksi jokin sopiva
kokoelma teht�vi� mahd. t�ydennettyn� ja/modifioituina soveliaasti.

\subsubsection*{Stabiilisuus, kankeus (stiffness)}

Eulerin menetelm�n stabiilisuusehto stabiilille yht�l�lle 
($f_y <0$) on $h < 2/f_y.$



\begin{verbatim}
Eulerh:=proc(f,vali,ya,h)
local a,b,n,m,t,y;
a:=vali[1]; b:=vali[2];
t[0]:=evalf(a);y[0]:=evalf(ya);
m:=ceil((b-a)/h);
for n from 0 to m do
  y[n+1]:=y[n]+h*f(t[n],y[n]);
  t[n+1]:=t[n]+h;
end do;
[seq([t[n],y[n]],n=0..m)];
end:
\end{verbatim}


\end{document}


