%\chapter{Ohjelmointi}\label{ohjelmointi}
\section*{Matlab- ja Maple- ohjelmointi}\label{ohjelmointi}

\textbf{Perusasioita}

\today

%\index{ohjelmointi}

%Koska osa teekkarikansastakin n�yt�� vieraantuneen ohjelmoinnista,
%pyrin selvitt�m��n juurta jaksaen, mist� (ja vain mist�) yleens�
%ohjelmoinnissa (korkean tason kielell�) on kyse.

Yleist� filosofiaa ja erityisesti Numsym05-kurssin tarpeita
palvellee parhaiten, jos esit�n asian rinnakkain \maple:n ja \mlb:n 
kannalta. 

Ohjelmien kirjoittaminen merkitsee systeemin laajentamista
omilla funktioilla, jotka periaatteessa toimivat aivan kuten
systeemin (\maple\/,\mlb\/, ym.) omat funktiot.

Sek� \maple\/ ett� \mlb\/ sis�lt�v�t suuren m��r�n kirjastofunktioita, joiden koodi on avointa. Hyv� tapa saada ohjelmointioppia 
on tutkia ohjelman valmistajan omia koodeja. Sit� ennen kannattaa
harjoitella lukemaan ja kirjoitamaan helppoja ja riisuttuja
muutaman rivin funktioita.

T�ss� samaistamme termit ``ohjelma'', ``proseduuri'' ja ``funktio''. Sek� \maple\/ ett� \mlb\/ ovat {\em funktionaalisia} ohjelmointikieli�, eli yleens� ohjelmat ``ovat funktioita''.

No, asiaan:

\textbf{Maple}


Olemme jo harrastaneet ohjelmointia kirjoittaessamme 
funktiom��rityksi�.

Tyyppi� \verb_ f:=x->x^2 _  olevat m��ritykset ovat jo tuttuja.
\begin{verbatim}
> type(f,procedure);
                                     true
\end{verbatim}
Toden totta, kirjoitimme jo ohjelman.

Vastaavasti voidaan m��ritell� useamman muuttujan funktio 
vaikkapa $f(x,y,z)=x^2+y^2+z^2$ n�in
\begin{verbatim}
   > F:=(x,y,z)->x^2+y^2+z^2;
\end{verbatim}

Jos halutaan lokaaleja muuttujia ja yleens� laajempaa koodim��r��, kirjoitetaan
m��rittely k�ytt�en {\tt proc}-syntaksia.

\index{proseduuri} \index{proc}
Edell� oleva $F$ voitaisiin m��ritell� t�t� pitemp�� syntaksia
k�ytt�en n�in:

\begin{verbatim}
  > F:=proc(x,y,z)
    x^2+y^2+z^2
    end:
\end{verbatim}

T�llaisessa {\tt proc}-m��rittelyss� palautetaan tuloksena viimeiseksi 
suoritetun komennon tulos (ellei {\tt RETURN}-lauseella toisin 
m��r�t�).

{\it {\bf K�yt�nn�n neuvo}:} Kun kirjoitat koodirivej� \maple-
ty�arkille, k�yt� \texttt{SHIFT-ENTER-} yhdistelm�� rivinvaihdoissa!}

Proseduuri on \maple-objekti siin�, miss� yksinkertainen muuttuja, lista,
matriisi jne. Se siis sijoitetaan arvoksi jollekin symbolille, jos sit�
halutaan nimell� kutsua.

Proseduurin yleinen muoto on t�llainen

\begin{verbatim}
   proc(argjono)
      local muuttujajono1;
      global muuttujajono2;
      options muuttujajono3;
      komento_1;
      komento_2;
      ...
      komento_n
   end:
\end{verbatim}

\index{local} \index{proc!local}
Sanoilla \texttt{local, global} ja \texttt{options} alkavat rivit 
eiv�t ole pakollisia.

{\bf Viimeiseksi laskettu} tulos on se, jonka {\bf proseduuri palauttaa}, 
ellei {\tt RETURN}-
lauseella toisin m��r�t�.
Kaikki lauseet p��tet��n puolipisteeseen, paitsi viimeinen, {\tt end}-sanaa
edelt�v� (siihenkin saa laittaa puolipisteen). Loppumerkkin� on normaaliin
tapaan puolipiste tai kaksoispiste.


\textbf{Matlab}


Sama \mlb:ll�:

1) ns. inline-m��rittely tehd��n \mlb-istunnossa:

\begin{verbatim}
>> F=inline('x^2+y^2+z^2','x','y','z')
 
F =
 
     Inline function:
     F(x,y,z) = x^2+y^2+z^2
>> F(1,2,3)
 
ans =
 
    14
\end{verbatim}

\mlb:ssa \maple:n proc:ia vastaava funktiom��ritys kirjoitetaan
suoraan teksititedostoon, ei siis \mlb-istuntoon.
Tiedosto on tyyppi� {\tt .m} ja sen tulee sijaita joko
oletushakemistossa (jonne voidaan siirty� \mlb:n sis�ll�
\unix-komennolla {\tt cd}). Yleisp�tev�mpi tapa on k�ytt��
{\tt addpath}-komentoa, jolla \mlb:n hakupolkuun lis�t��n
kyseinen hakemisto.

Tehd��np� tuo yll� sanottu nyt t�ll� tyylill�: Tiedostoon
{\tt F.m} kirjoitetaan jollain tekstieditorilla (luontevimmin
\mlb:n omalla editorilla) rivit:


\begin{verbatim} 
function w=F(x,y,z)
w=x^2+y^2+z^2;
\end{verbatim}

T�ss� on vastaava \mlb-istunto:
\begin{verbatim} 
>> pwd
 
ans =
 
/mnt/kiekko/home/apiola
 
>> addpath opetus/materiaali03/maple
>> F(1,2,3)
 
ans =
 
    14 
\end{verbatim}

Peruserot:
\begin{itemize}
\item
\maple:n proc:ssa ei tarvitse sijoittaa funktion palauttamaa
arvoa millek��n muuttujalle. (Haluttaessa sen voi tehd�, jolloin
kyseinen muuttuja kirjoitetaan {\tt local}-sanan per��n paikallisten muuttujien jonoon.
\item
\mlb:n funktiom��ritys alkaa aina varatulla sanalla 
{\tt function}, jonka j�lkeen kirjoitetaan otsikkoriville paluumuuttujan nimi 
(esimerkissa {\tt w}. Se
on vapaasti valittavissa (ei tietenk��n saa olla sama kuin
jokin argumenttilistan j�sen). Viimeiseksi t�lle muuttujalle sijoitettu arvo on se, jonka funktio palauttaa. 
\end{itemize}

Vastaavasti kuin \mlb:ss�, \maple-proseduurit kannattaa kirjoittaa
tekstitiedostoon, kun niit� alkaa kerty�. Niiden hakeminen
eri ty�arkeilta k�y pian tuskalliseksi.
Toisin kuin \mlb:ss�, \maple-funktiot kannattaa koota samaan
tekstitiedostoon, joka luetaan \maple:een tarvittaessa.
Esim: tiedosto {\tt omatohjelmat.mpl}

\begin{verbatim} 
# Tiedosto omatohjelmat.m  1.2.2005
# T�nne alan ker�t� omia hengentuotteitani.
# Aloitetaan vaikka t�ll�:

minmax:=proc(lista)
#Haetaan listan pienin ja suurin alkio, palautetaan
# tulos listana muodossa [m,M]
local jono;
jono:=op(lista);
[min(jono),max(jono)];
end:
\end{verbatim} 

\textbf{Maple-istunto}

Luetaan nyt �skeinen tiedosto \maple:en ja kokeillaan siin�
olevan (toistaiseksi ainoan) funktion  
toimintaa.

\begin{verbatim} 
> read("/home/apiola/opetus/materiaali03/maple/omatohjelmat.mpl");
> print(minmax);

proc(lista)
local jono;
    jono := op(lista); [min(jono), max(jono)]
end proc

> minmax([seq(i,i=-10..10)]);

                           [-10, 10]
> minmax([seq(sin(i),i=-10..10)]);

                           [-sin(8), sin(8)]
\end{verbatim}



 
%\end{document}


\section{Ohjauslauseita}\label{ohjaus} 

Tulkkaavassa j�rjestelm�ss�, kuten molemmat \maple\/ ja \mlb\/ voidaan k�ytt��
ohjauslauseita suoraan istunnossa, tarvitsematta aina kirjoittaa
funktioita. 

{\tt for  do  end do}-
rakenteen
yleinen muoto on
\begin{verbatim}
 |for <name>| |from <expr>| |by <expr>| |to <expr>| |while <expr>|
                         do <statement sequence> end do;
 tai

   |for <name>| |in <expr>| |while <expr>|
                         do <statement sequence> end do;
 
\end{verbatim}
\index{do $\ldots$ end do}
N�iss� kunkin pystyviivaparin rajaama osa on valinnainen. Esimerkiksi

\verb_ for i to 10 do a[i]:=x^i end do; _

tarkoittaa, ett� k�ytet��n oletusarvoja {\tt from 1 by 1} ja j�tet��n
{\tt while}-osa pois.

{\bf Huomaa} {\tt do <statement sequence> end do;}, miss� voidaan ajatella,
ett� {\tt do} on alkusulku ja {\tt end do} loppusulku.
Vastaava pari on {\tt if} --- {\tt end if} .


{\bf Esim.}  {\em Tsebyshevin polynomit} voidaan m��ritell� rekursiokaavalla

$$T_0(x)=1, T_1(x)=x, T_{n}=2xT_{n-1}(x)-T_{n-2}(x)$$

Teht�v�n� on m��ritell� proseduuri n�iden laskemiseksi.
Huomautamme aluksi, ett� laskenta voidaan tehd� suoraan 
interaktiivisesti t�h�n tapaan:
\begin{verbatim}
  > t0:=1;t1:=x;n:=5:
  > for i from 2 to n do
    tn:=expand(2*x*t1-t0);t0:=t1;t1:=tn;end do;
  > i;  # katsotaan i:n arvo for-silmukan j�lkeen.
\end{verbatim}

Interaktiivisessa k�yt�ss� {\tt for -- do -- end do --} rakenteella on
indeksimuuttujaan liittyv� sivuvaikutus, se j�� arvoon $\textit{yl�raja}+1$


Proseduuriksi kirjoitettuna toteutus voisi olla t�llainen.

\begin{verbatim}
  T:=proc(n::integer,x)
  local t0,t1,tn,i;
  t0:=1;t1:=x;
  for i from 2 to n do
    tn:=expand(2*x*t1-t0);
    t0:=t1;t1:=tn;
  end do;
  tn;  # T�m� palautetaan tuloksena.
  end:
\end{verbatim}

Nyt voidaan kutsua vaikkapa \verb_ T(5,x); _
Lis�simme parametrille n tyyppitarkistuksen.
Siten esimerkiksi kutsu \verb_ T(x,5);_ tai \verb_ T(a,10);_
antaa asianmukaisen virheilmoituksen.

Koska silmukkaindeksi {\tt i} on m��ritelty lokaaliksi,
proseduurin kutsu ei muuta sen arvoa. T�t� voi kokeilla
vaikka n�in \verb_ > i:='i':T(7,z);i; _

\textbf{Sama \mlb:lla}

\mlb:ll� ei voida tehd� symbolista. Tehd��n nyt oikean tyylinen
toteutus. Alkukommenteissa n�kyy funktion toiminta. \mlb\/ on 
siit�(kin) hyv�, ett� alkukommentit tulevat n�kyvin t�ss� tapauksessa komennolla {\tt help tseb}.

Tyypillisesti \matl-funktio rakennetaan toimimaan 
vektoriargumentilla. Kun lasketaan useiden polynomian arvoja
annetun x-vektorin pisteiss�, on luonnollista palauttaa
tulos matriisina, jonka sarakkeen pituus on $x$ ja kunkin
sarakkeen arvot ovat ao. asteisen polynomin arvoja $x$-vektorin
pisteiss�. T�ss� muodossa tulos voidaan samalla antaa heti 
{\tt plot}:lle.

\begin{verbatim}
function Y=tseb(x,n)
% tseb  Tsebychevin poynomit
%   Y=tseb(X,n) laskee n ensimm�ist� Tseb-polynomia
%   vektorin x pisteiss�.
%   Tulosmatriisin Y k:s sarake antaa astetta k-1 olevan Tseb-
%   polynomin arvot.

Y=ones(length(x),n); % Ykk�smatriisi, n saraketta.
x=x(:);            % Varmistetaan, ett� x on sarakevektori.
  if n==1, return, end  % Pelkk� 1. Tseb. pol =1.

Y(:,2)=x;
for k=3:n
  Y(:,k)=2*x.*Y(:,k-1)-Y(:,k-2);
end
\end{verbatim}
T�h�n loppui funktion koodi. Huomaa, ett� {\tt end} ei ole
funktion loppusana, vaan se kuuluu {\tt for}-silmukkaan.
(Funktiolla ei ole mit��n ``lopputervehdyst�''.) 

Nyt sitten innolla kustumaan:
\begin{verbatim}
>> help tseb
 
  tseb  Tsebychevin poynomit
    Y=tseb(X,n) laskee n ensimm�ist� Tseb-polynomia
    vektorin x pisteiss�.
    Tulosmatriisin Y k:s sarake antaa astetta k-1 olevan Teseb-
    polynomin arvot.
 
>> x=linspace(-1,1);
>> Y=tseb(x,4);
>> plot(x,Y)
\end{verbatim}
Matriisin $Y$ ($100 \times 4$) sarakkeina ovat nyt asteita $0,1,2,3$ olevien
Tseb-polynomien arvot $x$-vektorin pisteiss�.


Sek� \maple\/- ett� \mlb\/- ratkaisut ovat omalla tavallaan tyylikk�it�. \maple:lla voidaan suorittaa symbolisesti polynomien k�sittely�. \mlb:lla taas saadaan numeeriset arvot tehokkaasti k�ytt��n,
ja numeerinen jatkok�sittely on ��rimm�isen helppoa ja tehokasta.

\end{document}


\section{Evaluaatios��nt�j�}

Ohjelman sis�isiin muuttujiin ei sovelleta t�ysevaluaatiota, t�m�
johtuu ennen kaikkea tehokkusn�k�kohdista.
Luettelemme joitakin s��nt�j�, jotka voivat v�hent�� yll�tystilanteita.
\begin{itemize}
\item
Lokaaleille muuttujille suoritetaan vain yhden tason evaluaaatio.
T�t� voi simuloida komennolla {\tt eval(nimi,1)}.
\item
Parametrilistan symboleja ei voi k�ytt�� lokaaleina muuttujina
\item
Globaaleille muuttujille tehd��n t�ysevaluaatio.
\item
Edelliseen tulokseen viittaava \%-merkki k�ytt�ytyy evaluoinnin suhteen
kuten globaali muuttuja, siis sille suoritetaan t�ysevaluaatio, siit�
huolimatta, ett� se on selv�sti paikallinen.
\end{itemize}

Tarkat s��nn�t k�yv�t ilmi kirjasta ~\cite[ss. 42--49]{maple:progr}.


N�ill� tiedoilla voidaan jo tehd� paljon hauskaa ohjelmistoa.
Esittelemme joitakin esimerkkikoodeja, jotka samalla n�ytt�v�t
\maple:n mahdollisuuksia eri matematiikan aloilla. 

\section{Diskretointi}


Edell� k�sittelimme funktion diskretointia ja taulukointia. Vaikka
sellainen onkin helppo rakentaa {\tt seq}-funktion avulla, voi tuntua
helpommalta k�ytt�� funktiota, joka laskee kaiken suoraan. \matl-k�ytt�jille
on hyvin tuttua k�ytt�� {\tt linspace}-funktiota. \label{linspace}
Katsotaan, miten 
ohjelmoisimme tarkalleen samanlaisen \maple:lla. (\maple-versio toimikoon 
my�s symbolisella askelpituudella.)

Ensimm�inen versio olisi suora funktiom��rittely:
\begin{verbatim}
 linspace:=(a,b,n)->[seq(a+i*(b-a)/n,i=0..n) ];
\end{verbatim}

\index{linspace}
Ainoa varsinainen kauneusvirhe on indeksimuuttuja {\tt i}, jota t�ss�
tyyliss� ei voitane lokalisoida.

Kun lis�ksi haluamme version, joka k�ytt��
oletusarvona $100$:aan osaan jakoa, niin joudumme tutustumaan muutamaan
uuteen asiaan. T�m� koodi lienee hyvin itsens� selitt�v�. Huomaamme, ett�
proseduurissa voidaan viitata tiettyihin varattuihin nimiin, kuten 
\index{proc} \index{proc!args}  \index{proc!nargs}
{\tt args} \index{args} 
ja {\tt nargs}, \index{nargs} 
jotka mahdollistavat vaihtuvanpituisen 
parametrilistan.

\begin{verbatim}
linspace:=proc() local i,n,a,b; 
          a:=args[1];b:=args[2];
          if nargs=2 then n:=100 else n:=args[3] end if;
          [seq(a+i*(b-a)/(n-1),i=0..n-1) ] end;
\end{verbatim}

\section{Numeerinen integrointi}\label{numint}

Aloitamme puolisuunnikass��nn�st�, joka lasketaan kaavalla
\index{numeerinen!integrointi evalf(Int())} 
\index{evalf(Int(...))}
$$\int_a^b f(x) dx = h(1/2f(a)+f(x_1)+f(x_2)+\ldots+
 f(x_{n-1}) + 1/2f(b))+R,$$

miss� $h={\frac{b-a}{n}}$  ja 

$R=-{\frac{(b-a)h^2}{12}}f''(\xi), \xi \in (a,b)$.


%** Evaluointis��nt�, yhden tason evaluaatio, ei t�ysevaluaatiota

%** Operaattorit, ota differenssi ja PeriodicExtender

%** Ohjelmankehitys: Tarkista ensin vuorovaikutteisesti oman idean toimivuus.

Kirjoitetaan proseduuriksi, joka paluttaa listan $\left[S,R\right]$, jossa
$S$ on approksimoiva summa ja $R$ virhetermi.
\label{trapetsi}
\begin{verbatim}
trapetsi:=proc(f,a,b,n)
  local h,S,R;
  h:=(b-a)/n;
  S:=evalf((h/2)*(f(a)+2*sum(f(a+j*h),j=1..n-1)+f(b)));
  R:=-(h^2*(b-a)/12)*(D@D)(f)(xi);
  [S,R]  # Palutetaan lista: [Intappr, virhe]
end: 
\end{verbatim}

\huom Integraalia approksimoivan painotetun summan ohjelmoiminen k�y
kutakuinkin yht� helposti numeerisella matriisikielell�, kuten \matl.
Sensijaan virhetermin ohjelmoiminen numeerisella kielell� ei onnistu.

\index{trapetsi}
\index{puolisuunnikaskaava}

Kokeilemme {\tt trapetsia} integraalin $\int_a^b e^x \cos x dx$ 
laskemiseen v�lill�
$\left[0,\pi \right]$         

\label{trapetsi}
\begin{figure}
  \begin{center}
    \leavevmode
     \epsfig{file=trapetsi.ps, height=16cm, width=14cm}
     \caption{\tt www.../maple/opas/trapetsi.mws }
    \label{fig:trapetsi}
  \end{center}
\end{figure}

Virheen yl�raja on kuvasta katsoen luokkaa $0.6$

T�t� voitaisiin puolestaan k�ytt�� l�ht�kohtana {\it Romberg}-integroitiin.


\begin{verbatim}
> IRlista:=seq(trapetsi(f,0,1,2^k),k=1..9); 
> Ilista:=map(z->z[1],IRlista);
\end{verbatim} 


J�t�mme aiheen edelleenkehittelyn lukijalle.

Huomautamme, ett� kirjastopakkaus {\tt student} sis�lt�� mm. funktiot\rv
{\tt middlesum, trapezoid, simpson}. Niiden kutsu poikkeaa edellisest�
siin�, ett� integroitava annetaan lausekkeena eik� funktiona, eik�
virhetermi� ole mukana.
\index{packages!student}

\index{trapezoid}
\index{simpson}


\maple:n omaa numeerista integraattoria kutsutaan tyyliin\rv
\verb_ evalf(Int(lauseke,x=a..b));_

\Ind{evalf(Int(...))}

%evalf(Int(...));
%infolevel[int]:=5:

%T�ss� voisi olla hyv� l�ht�kohta ns. Rombergin menetelm�lle.
%Viittaus \cite[s. 107]{ISR}

\section{Interpolaatio} 

Annettuun dataan $\left\{x_0,\ldots, x_n\right\}$ liittyv�t {\it Lagrangen} 
kertojapolynomit ovat astetta n olevia polynomeja $L_k$, jotka toteuttavat
ehdot $L_k(x_j)=\delta_{k,j}$, miss� {\it Kroneckerin } $\delta$ on
$1$, kun indeksit ovat samat ja $0$ muuten. \Ind{Lagrangen interpolointi}
\Ind{interpolointi}

Voimme kirjoittaa funktion {\tt lag}, joka palauttaa annettuun dataan
liittyvien {\it Lagrangen polynomien} jonon.

\begin{verbatim}
lag:=(xd,x)->seq(product((x-xd[k])/(xd[j]-xd[k]),k=1..j-1)*
product((x-xd[k])/(xd[j]-xd[k]),k=j+1..nops(xd)),j=1..nops(xd));
\end{verbatim}

Testataan t�t� kolmen datapisteen muodostamaan joukkoon.

\begin{verbatim}
> lagjono:=lag([x[0],x[1],x[2]],t);
\end{verbatim}
$$\textit{lagjono:= }
{\frac {\left (t-x_{{1}}\right )\left (t-x_{{2}}\right )}{\left (x_{{0
}}-x_{{1}}\right )\left (x_{{0}}-x_{{2}}\right )}},{\frac {\left (t-x_
{{0}}\right )\left (t-x_{{2}}\right )}{\left (x_{{1}}-x_{{0}}\right )
\left (x_{{1}}-x_{{2}}\right )}},{\frac {\left (t-x_{{0}}\right )
\left (t-x_{{1}}\right )}{\left (x_{{2}}-x_{{0}}\right )\left (x_{{2}}
-x_{{1}}\right )}}
$$
Interpolaatioteht�v�n $p(x_k)=y_k,\  \ k=0,\ldots,2$ ratkaiseva polynomi
voidaan nyt kirjoittaa:
\begin{verbatim}
p:=sum(y[i-1]*lagjono[i],i=1..3):
\end{verbatim}
Tarkistetaan:
\begin{verbatim}
> subs(t=x[0],p),subs(t=x[1],p),subs(t=x[2],p); 
\end{verbatim}
$$y_{{0}},y_{{1}},y_{{2}}$$

\textbf{Vaihtoehtoinen Lagrange-koodi}

Ohjelmointiteht�ville on yleens� monia erilaisia ratkaisuja. Katsotaan
t�ss� kenties hieman elegantimpaa, joskaan ei koodim��r�lt��m lyhyemp��
tapaa. T�ss�
kirjoitetaan ensin lauseke $(x-x_0)(x-x_1)\ldots (x-x_n)$ ja jaetaan se
   $(x-x_j)$:ll�. 
N�in saadaan osoittaja. Nimitt�j� saadaan t�st� korvaamalla x arvolla
   $x_j$. Koodi voitaisiin kirjoittaa t�h�n tapaan:
\begin{verbatim}
L:=proc(j,xd,x)
  local oso,nimi,i,j1;
  j1:=j+1;
  oso:=product((x-xd[i]),i=1..nops(xd))/(x-xd[j1]);
  nimi:=subs(x=xd[j1],oso);
  oso/nimi;
end:
\end{verbatim}

T�ss� siis kirjoitimme yhden Lagrangen polynomin kerrallaan ja
aloitimme indeksoinnin 0:sta.
Lagrangen kertojapolynomijono liittyen x-arvojen listaan {\tt xd} 
muodostetaan nyt n�in:

\verb_  > seq(L(j,xd,x),j=0..nops(xd)-1); _

\textbf{Valmis funktio {\tt interp}}

Huomautamme, ett� \maple:ssa on my�s valmiina funktio {\tt interp}, jota
t�ss� tapauksessa kutsuttaisiin n�in:

\verb_ > interp([x[0],x[1],x[2]],[y0,y1,y2],t);_

\index{interp}
T�m� funktio kertoo polynomin auki, mik� k�yt�nn�ss� ei useinkaan ole 
tarkoituksenmukaista. Se saattaa huonontaa numeerista laskentatarkkuutta,
jos polynomi on korkea-asteinen. Laskentatehokkuuden kannalta kannattaa
ainakin saattaa polynomi ennen numeerista arvojen laskentaa 
ns. {\it Hornerin muotoon} komennolla {\tt horner} \index{horner}.

\section{Newtonin menetelm�}\label{newton}

Tarkastelemme klassista {\it Newtonin menetelm��} yht�l�iden
ja yht�l�ryhmien ratkaisemiseen.

Newtonin iteraatioaskel lasketaan kaavasta

$$x_{n+1}=x_n - {\frac{f(x_n)}{f'(x_n)}}$$

Kyseess� on siten funktion $x\to x-{\frac{f(x)}{f'(x)}}$ iterointi.

\maple-toteutus on lyhyt. Kirjoitetaan iterointifunktio {\tt N}


\begin{verbatim}
N:=x->evalf(x-f(x)/D(f)(x)):
\end{verbatim}

Jos haluaisimme laskea vaikkapa funktion
$$x{e^{-x}}-\left |\sin(10\,x)\right |$$
l�hell� pistett� $x_0=2.5$ olevan juuren, 
m��rittelisimme ensin funktion
\begin{verbatim}
> f:=x->x*exp(-x)-abs(sin(10*x));
\end{verbatim}
Sitten voisimme iteroida muutaman kierroksen:
\begin{verbatim}
> x1:=(N@@5)(0.25);
                               x1 := .2921679186
> x1:=(N@@2)(x1);  
                               x1 := .2921679186
\end{verbatim}

Suppeneminen oli hyvin nopeaa. Huomaamme, ett� \maple:n derivoimisoperaattori
osasi k�sitell� my�s itseisarvofunktiota.

T�m� oli helppo, mutta ohjelmointia ajatellen hieman v��r�oppinen tapa, 
k�ytimme globaalia symbolia $f$, joka tosin n�in yksinkertaisessa 
vuorovaikutteisessa ohjelmointitilanteessa on riskit�nt�.

Elegantti, oikeaoppinen tapa esitet��n kohdassa {\em operaattoreita}

\section{Differentiaaliyht�l�iden numeriikkaa}

Tarkastelemme alkuarvoteht�v��
$y' = f(t,y), \ \ y(t_0)=y_0$

\subsection{Eulerin menetelm� ja Runge--Kutta}\label{euler}

Differentiaaliyht�l�iden alkuarvoteht�vien ensimm�inen numeerinen menetelm�
on {\it Eulerin menetelm�}. Siin� edet��n l�ht�pisteest� tangentin
(tunnettuun) suuntaan askeleen $h$ verran. Jatketaan sitten samalla tavoin
ottaen l�ht�pisteeksi edellisen askeleen maali. Kaavan muodossa se
etenee n�in:
$$y_{n+1}=y_n+h f(t_n,y_n)$$

Kirjoitetaan t�m� proseduuriksi {\tt Euler}, joka palauttaa tasav�liseen
diskretointiin liittyv�n ``Eulerin murtoviivan''.

\begin{verbatim}
Euler:=proc(f,a,b,ya,m) 
local n,h,t,y;
h:=evalf((b-a)/m);
t[0]:=a;y[0]:=ya;
for n from 0 to m do
  y[n+1]:=y[n]+h*f(t[n],y[n]);
  t[n+1]:=t[n]+h;
od;
[seq([t[n],y[n]],n=0..m)];
end:
\end{verbatim}
Jos soveltaisimme {\it Eulerin menetelm��} vaikkapa teht�v��n

$$y'=t-y^2, \\ y(0)=1,$$

niin voisimme edet� n�in:

\begin{verbatim}
f:=(t,y)->t-y^2;
e3:=Euler(f,0,5,1,3);
plot(e3);
\end{verbatim}

Kenties haluaisimme rakentaa animaation {\it Eulerin menetelm�n} 
k�ytt�ytymisest� askelpituuden puolittuessa.
Se k�y tosi helposti.
\begin{verbatim}
> with(plots):
> display({seq(plot(Euler(f,0,5,1,2^k)),k=1..5)},insequence=true);
\end{verbatim}
Piirr�mme kaksi otosta animaatiosta samaan kuvaan.
%% Viittaus kuvaan 1jerry
\begin{verbatim}
> sini:=plot(Euler(f,0,5,1,2^2),color=blue):
> puna:=plot(Euler(f,0,5,1,2^5),color=red):
> display({sini,puna});
\end{verbatim}
\label{euler}
\begin{figure}
  \begin{center}
    \leavevmode
     \epsfig{file=euler.ps, height=5cm, width=5cm}
     \caption{}
    \label{fig:euler}
  \end{center}
\end{figure}


{\bf Teht�vi�}

\begin{itemize}
\item
Implementoi Runge-Kutta menetelm� ja testaa, tee vertailuja Euleriin.

\item
Yleist� edell� olleet yht�l�systeemille.

\end{itemize}

Huomautamme, ett� niin Eulerin menetelm� kuin Runge--Kutta voidaan
yht� helposti toteuttaa esimerkiksi \matl-ohjelmalla. Katso {\tt www}-
sivultamme kohtaa \matl-toteutuksia. Yll� esitetyiss�
ei ole mit��n erityist� symbolilaskentahy�ty�.
Toki listojen lista tietorakenteena tekee k�sittelyn joustavaksi.

\Ind{Runge--Kutta}

%%%* Vertailuksi \matl-Euler *

%###################################################################
%# Runge-Kutta 4
%###################################################################
%RK4:=proc(f,a,b,ya,m)
%local n,h,t,y,k1,k2,k3,k4;
%h:=evalf((b-a)/m);
%t[0]:=a;y[0]:=ya;
%for n from 0 to m do
%  k1:=h*f(t[n],y[n]);
%  k2:=h*f(t[n]+h/2,y[n]+k1/2);
%  k3:=h*f(t[n]+h/2,y[n]+k2/2);
%  k4:=h*f(t[n]+h,y[n]+k3);
%  y[n+1]:=y[n]+(1/6)*(k1+2*k2+2*k3+k4);
%  t[n+1]:=t[n]+h;
%od:
%seq([t[n],y[n]],n=0..m);
%end:
%\end{verbatim}

\subsection{Implisiittinen Eulerin menetelm�}

\index{Eulerin menetelm�} \index{Eulerin menetelm�!implisiittinen} 
Yksinkertaisin {\em implisiittinen menetelm�} on ``takaperoinen'' Eulerin
menetelm� (``Backward Euler'').

Idea on lyhyk�isyydess��n se, ett� ei edet�k��n alkupisteen 
vaan loppupisteen tangentin suuntaan. Ongelmana on, ett� t�t�
suuntaa ei tunneta, vaan se on haettava iteroimalla.
Analyyttisesti on kyse algoritmista

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

miss� uusi arvo $y_{n+1}$ esiintyy yht�l�n kummallakin puolella.

Implisiittiset menetelm�t ovat stabiilisuusominaisuuksiltaan yleens�
parempia kuin eksplisiittiset, mutta vastaavasti mutkikkaampia
ohjelmoida.

%T�ss� on siis tuntematon $y_{n+1}$ ja muut ovat tunnettuja.
Ratkaistavana on siis kullakin $n$:n arvolla yht�l�

$F(y)=0$, miss� $F(y)=y-y_n-h f(t_{n+1},y).$

J�timme $y$:n indeksin $n+1$ pois helpottaaksemme yht�l�n vakioiden
ja tuntemattoman ($y$) erottumista.

Vaikka t�m�kin on puhtaasti numeerinen menetelm�, sen implementointi
sujuu helpommin ja tyydytyst� tuottavalla eleganssilla, kun voimme
m��ritell� sopivat paikalliset iterointifunktiot. J�t�mme my�s 
parametrilistaan luvun {\tt niter}, jolla voimme antaa kullakin
askeleella suoritettavien {\it Newton}-iteraatioiden lukum��r�n.

\begin{verbatim}
 BE:=proc(f,a,b,Y0,N,niter)
     local y0,t0,y1,polyg,G,F,h,i,t1;
     h:=evalf((b-a)/N);
     y0:=Y0;t0:=a;
     polyg:=[a,y0];
     for i from 0 to N do
       t1:=t0+h;
       F:=unapply(y-y0-h*eval(f(t1,y)),y);
       G:=unapply(simplify(y-F(y)/D(F)(y)),y);
       y1:=y0+h*eval(f(t0,y0));    # Iteraation alkuarvo
                                   # tavallisella Eulerilla
       y1:=(G@@niter)(y1);         # Newtonin iteraatio
       polyg:=polyg,[t1,y1];
       t0:=t1;y0:=y1;
     od:  
 [polyg];          # Palautamme pisteiden (listojen) listan.
 end:
\end{verbatim}

Lukija voi kokeilla seuraavan \maple-komentojonen suorittamista, kunhan
ensin {\tt BE}-funktio on lis�tty mukaan.
Istunto n�ytt�� dramaattisen eron tavallisen ja implisiittisen Eulerin
menetelm�n v�lill� t�ss� teht�v�ss�. Kyse on tavallisen Eulerin menetelm�n
stabiilisuusalueen rajallisuudesta.

\begin{verbatim} 
 > f:=(t,y)->-10*(t-1)*y;
 > ratk:=t->exp(-5*(t-1)^2); # analyyttinen ratkaisu (k�sin tai 
                             # dsolve-komennolla)
 > map(ratk,[seq(1+i*.25,i=0..4)]);
 > [BE(f,1,2,1,4,5),Euler(f,1,2,1,4)];map(matrix,%);
 > with(plots):Eulerkuva:=plot(Euler(f,1,4,1,4),color=GREEN):
 > BEkuva:=plot(BE(f,1,4,1,4,5),color=BLACK):
 > ratkkuva:=plot(ratk,1..5):
 > display({BEkuva,ratkkuva});
 > display(Eulerkuva);
\end{verbatim}

Huomautamme, ett� \maple:n omat numeeriset differentiaaliyht�l�ratkaisijat
sis�ltyv�t komennon {\tt dsolve(...,numeric);} piiriin.
Komennolla on lukuisia lis�valitsimia mm. erilaisten menetelmien 
k�ytt�miseksi.

Se on toteutettu edell� kuvatulla operaattorityylill�, mik� n�kyy
kutsussa syntyv�st� ilmoituksesta, joka n�ytt�� uuden proseduurin
syntyneen. K�sittelemme aihetta tarkemmin kohdassa ~\ref{odenum}.


%%%%-----------------------------------------------------------------


\section{Kompleksitason k�yr�integraali}
%\esim
Olkoon teht�v�n� kirjoittaa proseduuri, joka laskee annetun
kompleksifunktion k�yr�integraalin annettua kompleksitason
polkua pitkin. Suoritetaan integrointi numeerisesti.
Integraali palautuu yhden reaalimuuttujan integrointiin kaavalla
$$\int_a^b f(p(t))p'(t) dt,$$
miss� $p$ on polun parametriesitys ja $t$ on k�yr�parametri.
Ajatus voidaan toteuttaa t�h�n tapaan:

\begin{verbatim}
   cintf:=proc(f,p,a,b)
   # f on integroitava kompleksifunktio, m��ritelt�v� funktioksi
   # p on polun m��rittelev� kompleksifunktio, kuten yll�
   # a ja b ovat parametriv�lin alku- ja loppupisteet.
   #
    local intfun,t,dp;
    dp:=diff(p(t),t);
    intfun:=evalc(f(p(t))*dp;
    evalf(Int(intfun,t=a..b)));
    end:

\end{verbatim}

\begin{itemize}
\item
Integroitava funktio $f(p(t))p'(t)$ rakennetaan lausekkeena. Erityisen
vaivattomaksi t�m�n tekee se, ett� $f$ annetaan funktiona.
\item
Valitsimme numeerisen integroinnin kirjoittamalla {\tt evalf(Int(...));}
\end{itemize}

Jos haluaisimme laskea vaikkapa integraalin $$\oint_c z^2 dz$$ pitkin 
yksikk�ympyr�n keh��, voisimme menetell� n�in:

\begin{verbatim}
  > f:=z->z^2 ; p:=t->exp(I*t);
  > cintf(f,p,0,2*Pi);
\end{verbatim}
Tulokseksi saatava $0$ on sopusoinnussa {\it Cauchyn 
integraalilauseen} kanssa.

Kutsu voitaisiin suorittaa my�s kirjoittamalla nimett�m�t funktiom��ritykset
suoraan parametrilistaan:
\verb_ > cintf(z->z^2,t->exp(I*t),0,2*Pi); _ )

Avaruuden $\R^3$ k�yr�- ja pintegraaleja on k�sitelty kiintoisalla ja
opettavaisella tavalla {\it R. Israel}in kirjassa  
\cite[Lab 28, ss. 222--226]{maple:israel}

\section{Picard--Lindel�fin iteraatio}

\Ind{Picard--Lindel�f}
T�m� on tyypillinen symbolilaskenta-aihe joka osoittaa t�llaisen kielen
suuren matemaattisen ilmaisuvoiman. Implementoimme iteraation samaan
tyyliin kuin {\it Newtonin} menetelm�n yll�. Elegantimpaa olisi tehd�
operaattorityylill� {\tt MakePL}, kuten seuraavassa kohdassa 
\ref{operaattorit}.
T�m�kin on kyll� eleganttia. Toteutamme sen vaihteeksi lausekelinjalla,
mik� tekee koodista hieman vaikealukuisempaa, mutta k�yt�st� taas
hieman helpompaa. Kaikenlis�ksi selvi�mme mainiosti pelkill� 
funktiom��rityksill� ilman {\tt proc}-lausetta.

J�t�mme t�m�n v�lineen kokeilun ja edelleenkehittelyn lukijalle.
Annamme koodit ja ajoesimerkin.

\begin{verbatim}
picklind1:=(fxt,t,x,xfun,t0,x0)->
          x0+int(subs(x=xfun,t=tau,fxt),tau=t0..t):

PL:=xfun->picklind1(x^2+t^2,t,x,xfun,0,1);

(PL@@3)(1); 
seq((PL@@k)(1),k=1..4);
plot([%],t=0..1);

\end{verbatim}

Iteraation l�ht�funktiona on vakiofunktio $1$, eli alkuarvo.


\section{Operaattoreita}\label{operaattorit}

\maple:ssa on mahdollista m��ritell� proseduureja, jotka palauttavat
tuloksenaan proseduureja. 

%Otamme kaksi esimerkki�. Lis�� ***

\subsection{Newton-operaattori}

Edell� k�sitelty Newtonin iteraatioaskel voitaisiin toteuttaa ilman
globaalina symbolina roikkuvaa funktiota t�h�n tapaan.

\begin{verbatim}
> MakeNewtonstep:=proc(flauseke::algebraic,x::name)
>    local askel;
>    askel:=simplify(x-flauseke/diff(flauseke,x));
>    unapply(askel,x);
> end:
\end{verbatim}
Jos hakisimme funktion
$f(x)=x-2\,\sqrt {x}$ nollakohtaa, suorittaisimme ensin Newtonaskelfunktion
luomisen t�h�n nimenomaiseen funktioon liittyen, eli m��rittelisimme
funktion lausekkeen ja antaisimme 
ohjelmointiapulaisellemme siihen liittyv�n teht�v�n.
\begin{verbatim}
> f:=x-2*sqrt(x);
> N:=MakeNewtonstep(f,x);                          
\end{verbatim}
$$
\textit{N:= }x\mapsto {\frac {x}{\sqrt {x}-1}}
$$
Iteroimme ensin $5$ kertaa l�htien alkupisteest� $2.0$ (huomaa desimaalipiste).
Sitten lis��mme laskentatarkkuutta ja iteroimme toiset $5$ kertaa.
\begin{verbatim}
> x0:=2.0:x0:=(N@@5)(x0);
\end{verbatim}
$$                   \textit{x0 :=} 3.999999999  $$
\begin{verbatim}
> Digits:=20:
> x0:=(N@@5)(x0);
\end{verbatim}
                     $$\textit{x0 :=} 4.0000000000000000000 $$

\subsection{Jaksollinen jatko}

Seuraavaksi esittelemme jaksollisen laajentajan, joka on erityisen
hy�dyllinen Fourier-sarjojen yhteydess�. Kts. lukua ~\ref{fourier}.
Se on toteutettu yleisluonteisella, vaikkakin ensialkuun hieman
vaikeaselkoisen n�k�isell� tavalla. Siihen voi suhtautua ``mustana
laatikkona''. Toisaalta, jos hetkiseksi pys�htyy koodia katselemaan,
niin huomaa, ett� siin� on sovellettu {\tt subs} komentoa {\tt proc}-
objektiin. Komennossa siirret��n yl�tason ``ohjelmointiapulaiselle''
annetut parametrit merkkijonotyylisesti palautettavan proseduurin koodiin.
T�ll� tavoin v�ltet��n kaikki evaluointij�rjestykseen liittyv�t 
ongelmat. (Kyse ei ole varsinaisesti merkkijonojen k�sittelyst�, k�ytetyt
lainausmerkit est�v�t evaluoinnin.)

Koodissa on suorastaan ylti�p�inen yksityiskohta, liek� tahallinen.
Siin� k�ytet��n varattua symbolia {\tt D}, mutta siihen ei kertaakaan
sijoiteta mit��n.

\index{PeriodicExtender}
\index{jaksollinen laajennus}
\begin{verbatim}
##########################
#
# Jaksollinen laajennus: MapleTech Vol 3 No. 3 1996
#
PeriodicExtender:=proc(f,d::range)
  subs({'F'=f,'L'=lhs(d),'D'=rhs(d)-lhs(d)},
  proc(x::algebraic) local y;
    y:=floor((x-L)/D);
    F(x-y*D);
  end)
end:
\end{verbatim}

T�m� on todella k�tev� ty�kalu. Ohjelmointiapulaisellemme siis annetaan
funktio ja haluttu jaksov�li. Apulaisemme ohjelmoi meille koko reaaliakselilla
m��ritellyn funktion t�lle nimenomaiselle l�ht�funktiolle ja jaksov�lille.

Esimerkiksi
\begin{verbatim}
kanttiaalto:=PeriodicExtender(signum,-1..1);
plot(kanttiaalto,-4..4);
\end{verbatim}
Saadaan kuva ~\ref{fig:kantti}.

\label{kantti}
\begin{figure}
  \begin{center}
    \leavevmode
     \epsfig{file=kantti.ps, height=5cm, width=5cm}
     \caption{}
    \label{fig:kantti}
  \end{center}
\end{figure}



%Esim. 3  (ehk� my�hemmin) T�ss� voisi Newton olla sopivampi
%Rekursiivinen funktio, option remember
%\begin{verbatim}
%
%  Fibonacci:=proc(n::nonnegint)
%    option remember;
%    if n<2 then RETURN(n) fi;
%    Fibonacci(n-2)+Fibonacci(n-1);
%  end:

%\end{verbatim}
%Ohjelmointia k�sittelemme hieman laajemmin \ref{}

%Ohjelmointia k�sitell��n perusteellisesti
%mm. viitteiss�  \cite{[progrgui]} ja \cite{maple:heck}.
%Joitakin ohjelmakoodeja on ~\ref{..}.


\section{Ohjelmankehitys} \label{ohjelmankehitys}

Jos koodi kirjoitetaan ty�arkille,
kannattaa ennen uudelle riville siirtymist� n�pp�ill�
\fbox{\tt Shift-Enter} (tai jatkaa kirjoittamista rivin lopun yli). 
Lopuksi sitten pelkk� \fbox{\tt Enter} suorittaa m��rittelyn (tai kertoo
syntaksivirheest�).

\index{ohjelmankehitys}
\index{ohjelmankehitys!debug}
\index{ohjelmankehitys!stopat}
\index{ohjelmankehitys!unstopat}
Koodin virheiden etsint��n ja tutkimiseen on komennot {\tt stopat}, 
\index{debug} 
{\tt unstopat}, {\tt DEBUG} ym. \Indmp{stopat}

%*** ota jokin debug-esimerkki ***


Toinen ohjelmankehitystekniikka on pit�� koodi tekstieditorissa 
{\tt koodi.txt} ja
kirjoittaa \maple-ty�arkille komento % {\tt read "koodi.txt;"}, 
\verb_ read "koodi.txt";_ 
joka voidaan
aina k�yd� suorittamassa uudelleen.

\index{read}
\index{luku ja kirjoitus!read}

\paragraph*{Ohjelmakoodien talletus} kannattaa yleens� tehd� tavalliseen
\asc-tiedostoon. T�ll�in ohjelmat voi vaivattomasti ja nopeasti lukea 
tarpeen mukaan. 
Toki koodi tallettuu my�s ty�arkin talletuksen yhteydess� ja
saadaan k�ytt��n avattaessa ty�arkki, kun muistetaan suorittaa
ty�arkilla olevat m��rittelyt. T�llainen k�ytt�tapa soveltuu 
tyypilliseti muutaman
pikaistuksissa kirjoitetun proseduurin talletukseen.

\asc-tiedostoon talletus on aina turvallisin tapa taata koodin s�ilyvyys
ja luoksep��sy tulevaisuutta ajatellen.





