\documentclass[finnish,11pt,twocolumn]{article}
% PC su 20.3.05
\usepackage{/home/apiola/tex/numsym05}
%\usepackage{numsym05}
\usepackage{ae,aecompl} 
\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}{12--13}{31.3.}

%{\bf Versio 1} ma 21.3.05

\begin{verbatim}
 http://www.math.hut.fi/teaching/numsym/05/maple/ns05.mpl
 http://www.math.hut.fi/teaching/numsym/05/L/
 http://www.math.hut.fi/teaching/numsym/05/matlab/
 http://www.math.hut.fi/teaching/numsym/05/maple/
 Cooper: Ch 4 ss. 111--153 (jaossa ti. 22.3.)

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

%Muista: uusi kehote: CTR-J tai CTR-K, miss� J, niinkuin "j�lkeen".


Vapautuaksemme Maple-Matlab-yhdistelm�n ``steady-state''-ongelmien 
mahd. 
aiheuttamista turhautumista, otamme nyt jo vakiintuneet
v�lineet hy�tyk�ytt��n. Dokumenttiosat ovat alla.
Tarkoitus on, ett� kukin voisi peruspalikoista (kuten writedata)
l�htien rakennella omia v�lineit��n joustavasti. 

{\bf N�m� k�yd��n l�pi ti 5.4}

\begin{enumerate}

\item 
Arvioitaeesa l�mp�yht�l�n analyyttisen ratkaisun tarkkuutta summausrajan $N$ funktiona, joudutaan pohtimaan alkuehtofunktion Fourier-sarjan
suppenemisnopeutta. Se riippuu AA-funktion $f$ sileydest� ja siit�,
toteuttaako $f$ my�s reunaehdot.

(a) Osoita (kahdella osittaisintegroinnilla), ett� jos 
$f:\lbrack 0,L\rbrack \to \R$ on $C^2$ (kaksi jatkuvaa
derivaattaa), niin sinisarjan kertoimet ovat $O(\frac{1}{n^2}).$
(Voit ottaa $L=\pi.$)

(b) Osoita, ett� sama k�yt�s saattaa synty� lievemminkin oletuksin
tarkastelemalla esimerkki�.
$ f(x)=\begin{cases}
x, \ \ x \leq L/2\\
L-x, \ \ x > L/2
\end{cases}
$

(c) Olkoon $f(x)=\sin (\frac{\pi x}{5})\, e^{x/2}$, v�lill� 
$\lbrack 0,10\rbrack.$
Muodosta numeerisesti sinisarjan kertoimien jono (ainakin arvoon $50$
saakka, ehk� jopa $100$, kunhan jaksat v�h�n odottaa) 
funktion {\tt siniSkertoimet} avulla. (T�t� teht�v�� varten ei tarvitsis kirjoittaa tiedostoon, mutta kirjoitapa valmiiksi teht�v�n 3 tarpeisiin.)\\
Tutki kokeellisesti (a)-kohdan v�itett�. N�ytt��k� k�yt�s pikemminkin
jopa $O(\frac{1}{n^3})-$tyyppiselt�?

\item

Olkoon sauvan pituus $L=10$, l�mp�yht�l�ss� $0$- reunaehdot ja alkuehtona
$f(x)=-\sin (\frac{\pi x}{5})\, e^{x/2}$ %e^{\frac{x}{2}}.$

(a) Muodosta \maple:lla analyyttinen ratkaisu. 
Huomaa, ett� {\tt lamposin}-funktio on onnen omiaan, koska integroinnit eiv�t onnistu
muuten kuin numeerisesti. 

(b) 
Luennolla (ti 22.3.) johdettiin virhekaava:

%\end{document}

$$| u(x,t)-u_{N-1}(x,t) | 
\leq \max_{n\ge N}
|b_n|
\frac{\exp(-(\frac{N \pi}{L})^2\, t)}
{1-\exp(-N(\frac{\pi}{L})^2\, t)}$$

Selvit� t�m�n arvion avulla, (a)-kohdassa laskettuja kertoimia
hy�dynt�en, miten monta termi� summaan kannattaa ottaa, jotta
maksimivirhe olisi vaikkapa $\leq 10^{-6}$, kun \\
(1) $t\ge t_0 $ ja (2) $0 \leq t \leq t_0$, miss� $t_0$ on sopiva
pieni aika-arvo (vaikkapa $t_0 = 0.1$).
(Voit s��dell� virhetoleranssia ja kynnysaikaa $t_0$ sopivammiksi
tilanteen mukaan.)





%Tutki kokeellisesti, miten monta termi� osasummaan on hyv� ottaa, jotta ratkaisu saataisiin "riitt�v�ll�" 
%tarkkuudella ainakin sopivan ajan kuluttua. My�s kerroinjonon pienenemist� on
%hyv� tarkkailla. \\

%Huomaa, ett� lamposin-funktio on onnen omiaan, koska integroinnit eiv�t onnistu
%muuten kuin numeerisesti. 
%(Tuskin kannattaa ottaa paljon enemp�� kuin $10$ termi�.)

%(b) Kuvassa kannattaa kiinnitt�� huomiota ep�symmetrisen alkuehdon 
%symmetrisoitumiseen, kun aika kasvaa. Olisko t�h�n jokin (fysikaalinen) 
%selitys tarjolla?

\item
Jatka edellist� teht�v�� \matl:ssa.\\
K�yt� edell� \maple:lla kirjoittamaasi kerroinvektoria hyv�ksesi 
{\tt lamposin.m}-funktiossa. 
Laske analyyttinen ratkaisu ottamalla osasummaan termej� virhearvioiden mukaan. (Tarvitset luonnollisesti enemm�n termej�, kun haluat tarkastella tilannetta
pienill� $t:$n arvoilla.)

Laske lis�ksi numeerinen ratkaisu {\tt CrNich.m}-funktiollamme ja 
vertaa analyyttiseen.

Laske maksimivirhe, ja katso, miten sen k�y ajan kasvaessa. 
Visualisoi.



\item
Olkoon t�ll� kertaa alkuehtona $f(x)=0$ ja reunaehtoina:\\
$u(0,t)=\frac{t}{1+t} + t e^{-\frac{t}{5}}$, $u(10,t)=0.$

(a) Muodosta ratkaisu {\tt CrNich}-funktiolla ainakin ajanhetkeen $t=10$ saakka,
ehk� pitemp��nkin. N�ytt��k� silt�, ett� ratkaisut l�henev�t ajasta riippumatonta, ``steady-state''-ratkaisua.

(b) Tutki reuna-arvofunktion raja-arvoa, kun $t\to\infty$. P��ttele sen
perusteella, mik� on tuo ``steady state''-ratkaisu $U(x)$.

(c) Totea aikasiivukuvista maksimiperiaatteen voimassaolo.

%(d) Miten suuri $t$:n pit�� olla, jotta $\max_x |u(x,t)-U(x)| < 0.1$ ?

\item
{\em Neumannin RE:t} Tarkastellaan l�mp�yht�l�� sauvassa, jossa on
eristetyt p��t, ts. $u_x=0$ sauvan p�iss�.
Palauta mieleen ratkaisu
$$u(x,t)=\sum_{n=0}^\infty 
a_n\cos\frac{n\pi x}{L}e^{-\lambda_n^2 t}, \ \ 
\lambda_n=\frac{c\,n\,\pi}{L}\ \ (k=c^2).$$

(a)
Muodosta analyyttinen ratkaisu, kun alkuarvofunktio\\
$f(x)=\exp(-(x-5)^2)$. $L=10, c=1.$
Voit mielell��n ottaa mallia {\tt siniSkertoimet} ja/tai
{\tt lamposin}-funktioista. 

%T�ss� voit toisaalta k�ytt�� hyv�ksesi {\tt erf}-funktiota
%(vaikka niinkin, ett� muodostat \maple:lla kaavan ja halutessasi
%jatkat ty�t� samalla {\tt erf}-kaavalla \matl:ssa.

(b) Ratkaise teht�v� {\tt pdepe}-funktiolla. \\
Malli:
{\tt http://www.math.hut.fi/opetus/numsym/05/matlab/pdelampoNeumann.m}\\
P��ttele analyyttisen ratkaisun perusteella, mit� arvoa $u(x,t)$
l�henee, kun $t\to \infty$ ja selvit�, ett� numeerinen ratkaisu
n�ytt�� k�ytt�ytyv�n samoin.



\item
Tarkastellaan l�mp�yht�l�� reunaehdoin \\ $u_x(0,t)=0, \ \ 
u_x(L,t)+ h\, u(L,t)=0,$ miss� $h > 0.$

Muuttujien erottelu johtaa $x-$ yht�l�n osalta ominaisarvoteht�v��n

$$-X''(x)=\lambda X(x),\ \ X'(0)=0,\ \ X'(L)+h X(L)=0.$$

T�ss� on merkitty $-\lambda:$lla muuttujanerotteluprosessin yhteist�
vakiota.

(a)
Osoita, ett� reunaehdot johtavat yht�l��n
$\tan s = \frac{L h}{s},$ miss� $s=L\sqrt{\lambda}$. Ominaisarvot ovat
siis muotoa $\lambda_n =(\frac{s_n}{L})^2,$ miss� $s_n$ on tuon yht�l�n
ratkaisu.  Ominaisfunktiot ovat muotoa $X_n(x)=\cos \sqrt{\lambda_n}x.$

(b)
Olkoon $L=10$ ja $h=1.$
M��rit� $10$ ensimm�ist� ominaisarvoa ja piirr� vastaavat ominaisfunktiot.
Yht�l�n numeeriseen ratkaisuun voit k�ytt�� \maple-funktota {\tt fsolve}
tai \matl-funktiota {\tt fzero}. Alkuarvot yht�l�n numeeriselle ratkaisijalle
valitaan kuvaa hy�dynt�en.

(c) Kirjoita Maplella laskemasi ominaisarvot tiedostoon Matlab-
jatkok�sittelymahdollisuutta varten.


\end{enumerate}


\subsection*{Ohjeita, koodeja}

L�mp�yht�l�n analyyttiseen ratkaisuun tarvittavat Maple/Matlab-v�lineet ovat
{\tt ns05.mpl}-tiedoston {\tt lamposin} sek� samanniminen \matl-funktio.
N�m� ovat nyt {\em Fri Mar 18 16:13:42 EET 2005} sek� 
{\tt www.../matlab}:ssa ett� {\tt /p/edu/mat-1.192/}:ssa. 
(Kyll�p� onkin vaikeaa pit�� yll� monessa paikassa. 
Jos tuo polku sattuu lipsahtamaan v��r�ksi, ole yst�v�llinen ja vaihda se itse.)

\begin{verbatim}
function Z=lamposin(c,x,t,Lb)
% HA 18.3.05
%
% Lasketaan l�mp�yht�l�n u_t = c^2 u_xx sinisarjaratkaisu korkeusmatriisina
% Kutsu: U=lamposin(c,x,t,Lb)
% Input: 
%        c:   c (L�mp�yht�l�n c)  (vakio = c^2)
%        x:   x-pisteet
%        t:   aikapisteet
%        Lb:   Joko tiedoston nimi tai muuttuja Lb, jossa sauvan pituus ja 
%              kertoimet muodossa [L,b1,b2,...]
%             ( Tyypillinen muoto 'Lb' tai Lb)
% Output: U -  l�mp�tilamatriisi, rivit ovat x:n, sarakkeet t:n pituisia.
% Huom! Jos kertoimet luetaan tiedostosta, on tarkistettava ja tarpeen
%       mukaan muutettava tiedostopolku ja nimi yhteensopivaksi sen kanssa,
%       jonne Maple-funktiolla 'siniSkertoimet' kirjoitettiin.
% Esim:
% clear; close all
% 
% t=linspace(0,0.5,20); x=linspace(0,pi,50);
% U=lamposin(1,x,t,'HALb'); % b-kertoimet luetaan tiedostosta 'HALb'.
% surfc(x,t,U); colorbar % Pintapiirros
% figure
% plot(x,U)      % L�mp�tilaprofiilit t-vektorin 
%                % ilmaisemilla ajanhetkill�.
\end{verbatim}

\begin{verbatim}
lamposin:=proc(f,L,c,N,tiedosto)
description "L�mp�yht�l� u_t = c^2 u_xx, 0-RE't, u(x,0)=f(x)";
local b,lambda,x,t,n;
# Palautetaan osasummafunktio, jossa termit 1..N.
# Jos annetaan "tiedosto", kirjoitetaan /tmp/tiedosto:on [L,b1,b2,...,bN]
# Kirjoittamisen suorittaa siniSkertoimet. Jos haluat muuttaa polkua, editoi
#                           ^^t�t�^^
# Esim:
# S1:=lamposin(x->x,1,1,10):
# S2:=lamposin(x->x,1,1,20,"b"):
# plot3d(S1(x,t),x=0..L,t=0..0.1)  # sama S2:lle.
\end{verbatim}


\end{document}


\subsubsection*{Teht. 1}





Maple-malli visualisointiin ym. : {\tt L/L6.mws} (ei sis�ll� max/min-laskentaa).
Matlab:ssa voit vaikka m��ritell� inline-funktioiksi, joiden oltava x:n ja t:n suhteen
"array smart" (vrt. aiemmin jaettu Cooperin kirjan Matlab-appendix.). Sitten 
{\tt meshgrid}, ... \\
Kts. my�s  {\tt http://www.math.hut.fi/teaching/numsym/} "Luennot"-kohdan ohjeita
Luento 6:n kohdalta.


\subsubsection*{Teht. 2}

T�ss� on funktion {\tt sinkert} kutsu ja alkukommentit. Koodia t�ytyy varautua
muuttamaan ainakin tiedostopolun suhteen (Unix:ssa tuo koodissa oleva lienee aina hyv�.)
Voidaan k�ytt�� pelkk�n� Maple-funktiona tai Maple-Matlab-yhdistelm�n�.
\begin{verbatim}
siniSkertoimet:=proc(f,L,N,tiedosto)
# b:=siniSkertoimet(f,L,N,tiedosto); Funktion f sinisarjan kertoimet.
# b:=siniSkertoimet(f,L,N); palauttaa b-kertoimet, mutta ei kirjoita 
#                           tiedostoon.
# f -- funktio (esim. x->x^2) Oltava oikea Maple-funktio, ei lauseke.
# L -- jakson puolikas
# N -- Mihin saakka kertoimia lasketaan
# tiedoston nimi merkkijononona (ilman polkua), esim. "bkert.txt"
# Kertoimet palautetaan liukulukutaulukkona (array) ja kirjoitetaan tiedostoon,
# josta ne voi suoraan lukea Matlabiin load-komennolla: 
#       Lb=load('/tmp/bkert.txt')
# Kirjoitetaan vektori [L,b1,b2,...,bN].
# Ty�t� voi jatkaa joko Maplessa tai Matlabissa.

local B,n,t,omega,polkujatiedosto;
B:=array(1..N);
omega:=Pi/L;
for n from 1 to N do
   B[n]:=evalf(2/L*int(f(t)*sin(n*omega*t),t=0..L));
end do;
if nargs=3 then RETURN(B) end if;

#polkujatiedosto:=cat("C:/MATLAB6p1/work/",tiedosto);  # HA:n kannettava   
 polkujatiedosto:=cat("/tmp/",tiedosto);                 # Unix

writedata(polkujatiedosto,[L,seq(B[j],j=1..N)]); # Kirjoitetan my�s L.
print("Min�, siniSkertoimet kirjoitin [L,b]:n tiedostoon ",polkujatiedosto);
#print("Palautan myos b-kertoimien taulukon (array) t�h�n sinulle");
RETURN(B)
end proc:
\end{verbatim}

Vastaava Matlab-funktio on \\
{\tt www.math.hut.fi/numsym/05/matlab/Msinsarja}. Voidaan siis k�ytt��
pelkk�n� Matlab-funktiona tai Maple-Matlab-yhdistelm�n�. 
\begin{verbatim}
function [summa,osasummat] = Msinsarja(x,LtaiFile,b)
%  matlab/Msinsarja.m   
% Kutsut: 
% 1) [summa,osasummat] = Msinsarja(x,tiedosto)
%    Tiedostosta luetaan L ja b
% 2) [summa,osasummat] = Msinsarja(x,L,b)
%    Kerroinvektori annetaan Matlabissa.
% Koko suoritus:
% 1) Maplessa: ns05.mpl-funktio siniSkertoimet(f,L,N,tiedosto);
% Esim: B:=siniSkertoimet(x->x^2,1,10,"bkert");
% 2) Matlab:ssa: 
% x=linspace(-2,2); 
% [summa,osasummat]=Msinsarja(x,'bkert');
% plot(x,summa); figure; plot(x,osasummat)
polku='/tmp/';        % UNIX  % Muuta tarpeen mukaan.
%polku='C:\MATLAB6p1\work\';  % WIN  % Muuta tarpeen mukaan.
if nargin==2
    tiedosto=[polku LtaiFile];
    Lb=load(tiedosto);    
    L=Lb(1);     % ensimm�isen� on jakso/2.
    b=Lb(2:end); % Loput muodostavat b-vektorin.
else
    L=LtaiFile;
end;
N=length(b);
kertoimet=b*ones(size(x)); % Er�s monista tavoista tehd� matriisi, jonka
                            % 1. rivi on vakio a0
                            % 2. rivi on vakio a1,
                            % ...
termit=kertoimet.*sin(pi/L*(1:N)'*x);          
summa=sum(termit);
osasummat=cumsum(termit);
\end{verbatim}

\subsection*{Teht. 3}

\begin{verbatim}
lamposin:=proc(f,L,c,N,tiedosto)
description "L�mp�yht�l� u_t = c^2 u_xx, 0-RE't, u(x,0)=f(x)";
local b,lambda,x,t,n;
# Palautetaan osasummafunktio, jossa termit 1..N.
# Esim:
# S1:=lamposin(x->x,1,1,10):
# S2:=lamposin(x->x,1,1,20,"b"):
# plot3d(S1(x,t),x=0..L,t=0..0.1)  # sama S2:lle.
if nargin = 4 then
  b:=siniSkertoimet(f,L,N);  
else
  b:=siniSkertoimet(f,L,N,tiedosto);  
end if;
 lambda:=n->c*n*Pi/L;
 RETURN(evalf@unapply(sum(b[n]*sin(n*Pi*x/L)*exp(-lambda(n)^2*t),n=1..N),x,t));
 end proc:
\end{verbatim}

\subsection*{Teht. 4}

\begin{verbatim}
lampodiffE:=proc(f,nx,nt,dt)
 description"L�mp�yht, v�li: [0,1], l�mp�vakio =1 (vrt. [T-W])";
# f -- alkuarvofkt.
# nx,nt,dt ilmeisi�
# Tulos: Listojen lista, m:s lista edustaa ao. aikatasoa.
# Tulostetaan lis�ksi r=dt/dx^2. sen pit�� olla alle 1/2, jotta stabiili.
\end{verbatim}

\subsection*{Teht. 6}

Koodissa {\tt heat3.m} on ehk� selke�mp�� ja tehokkaampaa korvata A:n rakentaminen 
ennest��n k�ytetyll� tyylill�:
\begin{verbatim}
e=ones(J-1,1);
A=spdiags([-r*e (1+2*r)*e -r*e],-1:1,J-1,J-1);
\end{verbatim}
Sin�ns� on opettavaista tutustua sparse komennon moninaisiin muotoihin.

\end{document}


\subsubsection*{Teht. 1}
\matl:ssa voit jakaa grafiikkaruudun osiin {\tt subplot}:lla.
Toisaalta voit my�s piirt�� eri ikkunoihin komentamalla {\tt figure} 
(seuraava)
tai {\tt figure(n)}.
Piirrett�v�n datan voit ker�t� matriisiin {\tt for}-silmukan avulla
sarake tai rivi kerrallaan.
(Silmukkaa ei tehd� pitkin vektorin komponentteja --- sellainen saa
opettajat hyppim��n seinille!)

Toisaalta voidaan my�s tehd� ilman yht��n for-silmukkaa:

\textbf{Ulkotulo, sarake {\tt *} rivi, meshgrid}
Monissa yhteyksiss� tarvitaan vektorien $x$ ja $y$ ulkotuloa, joka on
matriisi $M,$ miss� $m_{i\,j}=x_i y_j.$
Jos $x$ ja  $y$ ovat vaakavektoreita, saadaan $M$ matriisitulona
{\tt M=x'*y}. Toinen tapa (esim)
\begin{verbatim}
 >> x=linspace(0,1,8);
 >> y=linspace(-1,1,5);
 >> [X,Y]=meshgrid(x,y);
 >> M=X.*Y;   % Katso:  X, Y, M
 >> mesh(x,y,M); surf(x,y,M)  % Kokeile!
\end{verbatim}
Huomaa: edellisess� matriisikertolasku {\tt(*)}, j�lkimm�isess� {\tt (.*)}

T�ll� voidaan korvata for-silmukka niin jatkuvassa kuin diskreetiss�kin.
(Jatkuva on Matlabin kannalta oikeasti my�s diskreetti, pisteit� on vain
enemm�n.)
Lue meshgridin k�yt�st� vaikka "pikku-oppaasta" tms. ja \\
{\tt help meshgrid, help mesh, help surf, help contour}.


\subsubsection*{Teht. 2}
\matl:ssa komento {\tt eig} laskee ominaisarvot (ja pyydett�ess� my�s vektorit). \verb_ >> [V,D]=eig(A)_ muodostaa "diagonalisoinnin", kts. {\tt help eig}.

Jos haluat verrata kaavasta laskettuja ja eig:n antamia, niin ominaisarvojen 
muodostaman vektorin saat edell� olevasta D:st�: {\tt diag(D)}.\\
Ominaisvektorit ovat matriisin $V$ sarakkeet. T�ss� kaikki ominaisavaruudet
ovat $1-$ulotteisia (miksi?). Jos k�yt�ss�si on kaksi matriisia ($V, W$), 
joiden
sarakkeet koostuvat ominaisvektorikandidaateista samassa j�rjestyksess�, 
voit selvitt�� vertailun suorittamalla {\tt V./W} (miksi). Ja viel� yksi
"miksi": Niin siis miksi on t�rke�� t�m�n toimivuudelle, ett� ominaisarvot 
ovat yksinkertaisia?

\subsubsection*{Teht. 3}

\matl-ratkaisussa voitaisiin vaikka \maple:n avulla muodostaa:

\verb_ b=inline('..','n')  _ n on skalaari tai indeksivektori, k�yt� piste-operaatoita
tai funktiota {\tt vectorize}.

Osasummia laskettaessa kannattaa muodostaa sarjan termien matriisi siten,
ett� termi $n$ muodostaa rivin n (eri x:n arvoilla).
Funktio {\tt sum} laskee matriisin sarakesummien vektorin, joka antaa nyt 
osasumman $s_N(x_1,...x_{100})$. ({\tt sum} on esimerkki ns. vektorifunktiosta,
muita on {\tt min, max, ...}. (Kts. laajempaa opasta (kunhan korjaan linkit).)
Esim: 
\begin{verbatim}
x=linspace(-2,2);
N=10;
termit=zeros(N,length(x));
for n=1:N
   termit(n,:)=b(n)*sin(...);
end;
summa=sum(termit);
\end{verbatim}

Huom! T�ss�kin for-silmukka voidaan korvata ulkotulolla, kuten yll�.
Tehokkuuden kannalta voi olla merkityst�, t�ss�h�n voidaan haluta laskea
isoillakin N:n arvoilla (ainakin Gibbsiss�).

Jos haluat osasummaparven, saat senkin samasta termimatriisista funktion
{\tt cumsum} avulla. Voit poimia j�lkimm�isest� sopivasti osia.

\subsubsection*{Teht. 5}
Funktion v�litt�minen argumenttina \matl-funktiolle: \\
Esim ({\tt [Hig-Hig]} Ch 10): 
Funktio {\tt eros} palauttakoon parametrina annetun funktion $f$ 
erotusosam��r�n: $\frac{f(x+h)-f(x)}{h}$ Voitaisiin kirjoittaa 
(tiedostoon {\tt eros.m}) n�in:
\begin{verbatim}
function y=eros(f,x,h)
y=(feval(f,x+h)-feval(f,x))/h;
\end{verbatim}
Siis se, mik� kirjoitettaisiin istunnossa {\tt f(x)}, kirjoitetaan funktion
rungossa {\tt feval(f,x)}. {\em {\bf Huom:} T�m� koskee vain tapausta, jossa {\tt f} 
esiintyy argumenttilistassa.}\\
Kutsuesim:
\verb_ eros(@sqrt,1,0.1)_  Kutsussa voi esiinty� $sqrt:$n paikalla 
yht� hyvin inline-m��rittely
tai m-tiedosto. (My�s \verb_ 'sqrt'_ toimii, mutta se alkaa n�ytt�� "vanhanaikaiselta", muoto \verb_@fun_ on ns. "function handle", olisko vaikka "kahva".)


\subsubsection*{Fourier-sarjat} 

Kts. yll� olevia referenssej� \maple- ja \matl-tekniikoihin.

$2 L$-jaksoisen funktion $f$ Fourier-sarja:

$$a_0+\sum_{n=1}^\infty (a_n \cos {\frac{n \pi x}{L}}+
 b_n \sin {\frac{n \pi x}{L}}),$$

$a_0={\frac{1}{2L}}\int_{-L}^L f(x)dx,\ \
a_n={\frac{1}{L}}\int_{-L}^L f(x)\cos {\frac{n \pi x}{L}} dx,\ \ 
b_n={\frac{1}{L}}\int_{-L}^L f(x)\sin {\frac{n \pi x}{L}} dx .$


\subsubsection*{L�mp�yht�l�n Fourier-sarjaratkaisu \matl:lla}

(Valtaosa seuraavan funktion riveist� on kommentteja.)
\begin{verbatim}
function [Z,x,T]=lamposin(L,c,b,T,M)
% HA 25.11.03 
% Lasketaan l�mp�yht�l�n u_t = c^2 u_xx sinisarjaratkaisu korkeusmatriisina
%
% Input: L:   jakso/2
%        c:   c (L�mp�yht�l�n c)  
%        b:   sinisarjan kerroinvektori (nollat mukana)
%        T:   aikapisteet
%        M:   x-diskrertointi (oletus 50)
%        
% Output: [Z,x,t] -  l�mp�tilamatriisi, x-vektori, t-vektori
%
% Esim:
% clear; close all
% L=10;c=sqrt(1.14);nb=10;
% b=1./(1:nb);b=rem((1:nb),2).*b;
% kerroin=400/pi; b=kerroin*b;
% T=linspace(0,10,10);
% [U,x,t]=lamposin(L,c,b,T);
% surfc(x,t,U);  % Pintapiirros
% figure
% plot(x,U)      % L�mp�tilaprofiilit t-vektorin 
%                %  ilmaisemilla ajanhetkill�.
%
N=length(T);
if nargin==4, M=50; end 
nb=length(b);
x=linspace(0,L,M);
Z=zeros(M,N);
for n=1:nb;
   Z=Z+(b(n)*sin(n*pi*x./L))'*exp(-((n*pi*c/L)^2)*T);
   % Vanha tuttavamme ulkotulo on t�ss� korvaamaton.
end;
Z=Z';
\end{verbatim}

\end{document}




\begin{enumerate}

\item
Suorita jokunen Fourier-sarjaksi kehitt�minen. Alla pari esimerkkifunktiota,
voit valita yht� hyvin jotain muuta.

Piirr� $f$:n kuvaaja muutaman jakson alueella. Muodosta $f$:n Fourier-sarja.
Selvit�, miss� pisteiss� sarja suppenee ja miss� se ei suppene kohti kyseist�
funktion arvoa. Mit� arvoa kohti sarja n�iss� ``ei-pisteiss�'' suppenee?
Piirr� jaksollinen laajennus ja osasummafunktioiden kuvaajia.

Tutki kokeellisesti {\em Gibbsin ilmi�t�} 
(n. 9\%:n "yli/ali-ampuminen" ep�jatkuvuuskohdan ymp�rist�ss�).

(Luultavasti ei kannata piirt�� hirve�� m��r�� osasummia, vaan valittuja
muutamia harvakseen, isoilla n:n arvoilla.)

Esim:

(a)
$f(t)=t, \ \ t\in (-\pi,\pi],\ \ f(t+2\pi)=f(t) .$

(b)
$f(x)=\begin{cases}
         100, \ \ -4 < x < -2 \\
         0, \ \ -2 \le x \le 4 \end{cases}$,
$f(x+8)=f(x)$.

\item
Tarkastellaan vaimennettua v�r�htely� kuvaavaa yht�l��

$my''+cy'+ky=r(t)$.

Olkoon $m=1, c=0.02, k=25$ ja $$r(t)=
\begin{cases}t+\pi/2, \ \ -\pi<t<0,\\
             -t+\pi/2, \ \ 0<t <-\pi  \end{cases},\ \ 
r(t+2\pi)=r(t)$$ .

M��rit� ``pitk�n aikav�lin'' (``steady state'') ratkaisu $y(t)$.
(T\"ast\"a on homogeeniosan (ja siten alkuehtojen) vaikutus vaimentunut
pois.)

Periaate: Kehitet\"a\"an $r(t)$ Fourier-sarjaksi. Etsit\"a\"an kutakin sarjan
termi\"a vastaten erityisratkaisua $y_n$ muodossa $y_n=A_n\cos nt+
B_n\sin nt$. Lineaarisuuden perusteella (jos sarjan termeitt\"ain 
derivoiminen
on sallittua) p\"a\"atell\"a\"an, ett\"a $\sum y_n$ on koko yht\"al\"on ratkaisu.

Laske ainakin 5 termi\"a ratkaisusarjasta ja piirr\"a niiden summa sek\"a 
termit erikseen. Jokin niist\"a on dominoiva, identifioi se. Piirr\"a
samaan kuvaan her\"ate $r(t)$ ja "steady state" ratkaisu v\"alill\"a 
$[0,4\pi]$.

Tutki, mit� tapahtuu, jos jousivakiolle annetaan arvo $k=9$ sek� toisaalta
j�yk�n jousen, vaikkapa $k=49$ tapauksessa. Ent�, jos vaimennusta lis�t��n?


\item 
Ratkaise AA-teht\"av\"a $x''+2x=r,\ \ x(0)=x'(0)=0$,miss\"a $r(t)$ on $2\pi$-
jaksoinen neli\"oaalto: $$r(t)=\begin{cases}-1,  0 < t < \pi \\
                                   1, \pi < t < 2\pi \end{cases}.$$
Suorita ratkaisu seuraavilla kahdella eri tavalla:
\kohta{a}{Muodosta her�tefunktio vaikkapa 3:n jakson alueella ja
ratkaise yht�l� Laplace-muunnoksella. Saat k�ytt�� \maple-funktiota
{\tt dsolve(...,method=laplace);}.

      }

\kohta{b}{
Kehit� $r(t)$ Fourier-sinisarjaksi ja k�yt� osasummaa $S_7$ 
(eli mukana ovat termit $\sin 7t$-termiin saakka.) oikean her�tteen
approksimointiin.
   }

\kohta{c}{
Vertaa molemmilla menetelmill� saamiasi tuloksia, esit� graafisesti
ja katso, miten tarkkuus paranee, kun otat korkeamman asteisia Fourier-
osasummia.
     }
\huom\/ Hy�dynn� tietysti ennen mainittua {\tt PeriodicExtender} (alias
 JJ) -funktiota.



\end{enumerate}


Referenssej�:

[Kof] Ch 25 pp. 351 -- 363  \\
[HAM] s. 177 ... (kts. HAM-sivun fourier-ty�arkki) Ota k�ytt��n
jaksollinen laajennus.\\
[KRE] `` Forced oscillations ``


\end{document}


\item %5
Tarkastellaan $2:$n pituusyksik�n mittaista sauvaa, jonka materiaali on
sellainen, ett� l�mm�njohtumisyht�l�n vakio $c=1$. 
 
Katkaistaan sauva keskelt� kahtia ja sijoitetaan toinen puolisko
kiehuvaan vesis�ili��n ($100^\circ C$). 
Toinen  puolisko laitetaan j��veteen ($0^\circ C$). Kun sauvan 
puolikkaat ovat saavuttaneet kauttaaltaan ymp�rist�ns� l�mp�tilan,
ne
sijoitetaan 
l�mp�eristeeseen, liitet��n takaisin 
yhteen, ja sijoitetaan vapaaksi j��v�t p��t j��vesis�ili��n ($0^\circ C$).

 
Suorita vaiheittain muuttujanerotteluprosessi sauvan l�mp�tilafunktion
m��ritt�miseksi.

 
Laske sauvan l�mp�tila hetkell� $t=0.1$ kummankin sauvan puolikkaan
keskipisteess� k�ytt�en sarjan kahta ensimm�ist� termi�.




\subsection*{Fourier-kaavoja}

\subsubsection*{T-jaksoisen funktion Fourier-sarja}

\paragraph*{Reaalimuoto}

Olkoon $f$ T-jaksoinen reaalifunktio, miss� $T=\frac{2\pi}{\omega}$. $f$:n
Fourier-sarjaesitys: \\
$f(t)=\frac{1}{2}a_0 + \sum_{n=1}^\infty a_n \cos n\omega t+
\sum_{n=1}^\infty b_n \sin n\omega t$,
miss�

$$a_n=\frac{2}{T}\int_d^{d+T}f(t)\cos n\omega t dt, \ \ n= 0, 1, 2, \ldots $$
$$b_n=\frac{2}{T}\int_d^{d+T}f(t)\sin n\omega t dt , \ \ n= 1, 2, \ldots$$

\paragraph*{Kompleksimuoto}
$f(t)=\sum_{n=-\infty}^{\infty} c_n e^{i n \omega t}$, miss�

$$c_n=\frac{1}{T}\int_d^{d+T}f(t)e^{-i n \omega t}dt,\  \ 
n=0,\pm 1,\pm 2 \ldots$$


\end{document}
\subsection*{Fourier-kaavoja}

\subsubsection*{Yleinen tapaus, $f$ m��ritelty v�lill� $(-L,L)$}
$f(x)=a_0+\sum_{n=1}^\infty a_n \cos {\frac{n \pi x}{L}}+
\sum_{n=1}^\infty b_n \sin {\frac{n \pi x}{L}}$ \rv
$a_0={\frac{1}{2L}}\int_{-L}^L f(x)dx, \ \
a_n={\frac{1}{L}}\int_{-L}^L f(x)\cos {\frac{n \pi x}{L}} dx ,\ \
b_n={\frac{1}{L}}\int_{-L}^L f(x)\sin {\frac{n \pi x}{L}} dx$

\subsubsection*{Puolen jakson laajennukset,$f$ m��ritelty v�lill� $(0,L)$}

\paragraph*{Sinisarja}

$$f(x)=\sum_{n=1}^\infty b_n \sin {\frac{n \pi x}{L}}, \ \
b_n={\frac{2}{L}}\int_0^L f(x)\sin {\frac{n \pi x}{L}} dx 
\qquad (HS_1, HS_2) $$

\paragraph*{Kosinisarja}

$$f(x)=a_0+\sum_{n=1}^\infty a_n \cos {\frac{n \pi x}{L}}, \ \
a_0={\frac{1}{L}}\int_0^L f(x)dx, \ \
a_n={\frac{2}{L}}\int_0^L f(x)\cos {\frac{n \pi x}{L}} dx 
\qquad (HC_1 - HC_3) $$


%\end{enumerate}





\end{document}





\item
Ratkaise edellinen teht�v� Laplace-muunnoksella katkaisemalla her�te
vaikkapa 3:n jakson alueelle (tai j�t� jaksojen lukum��r� muuteltavaksi
parametriksi).

Vertaa ratkaisua edelliseen.





\item
M��rit� tasapainotilaratkaisu ("steady state") v�r�htelysysteemille
 $m x'' + c x' + kx = r(t)$, miss� $m=1, c=0.02, k=25$,
 ja $r(t)=\begin{cases}t+\pi/2, \ \ -\pi<t<0,\\
                                      -t+\pi/2, \ \ 0<t <-\pi 
\end{cases},\ \ r(t+2\pi)=r(t)$ .

Tutki lis�ksi, miten vasteen amplitudit k�ytt�ytyv�t, jos jousivakio pienenee
aina arvoon $k=9$ ja toisaalta, jos jousta j�ykistet��n arvoon $k=49$ saakka.
(Vrt. Kreyszig-ote).
