\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{7}{15--16}{6 -- 15.4. 2004}

\textbf{Lisäyksiä, korjauksia ja tarkennuksia 6.4.}

Niitä on tehty edelliseen versioon verrattuna, mainittakoon
tässä nämä:

Teht. 2 jätetään Greenin funktio, (c)-kohta pois. \\
Tehtävässä 3 on tarkoitus ratkaista 11.5 (a),
sitä täsmennettiin (helpotettiin). \\
Lisättiin tehtävä, joka numeroitiin 4:ksi, siirrettiin vanha
4 tehtäväksi 5.

Tehtäviä 3 ja 4  voidaan pitää päätehtävinä.



%PDE-aiheista irtoaisi helposti projektityö vaikka Molerin 
%tehtävien ja materiaalin pohjalta.

\hrule

Jatketaan yhdessä tekemistä ti 6.4,
ja käydään ratkaisut läpi pääsiäisen jälkeen (to 15.4.)

Sovitaan projekteista ja aikataulusta mahdollisuuksien mukaan 
ti 6.4.
Osa projekteista voi hyvin olla osdiffyht. aiheeseen liittyvää,
vaikkapa pääasiassa tässä oleviin lähteisiin pohjautuvaa 
``projektiainesta''. Siitä riittää useammallekin tekijälle.

Lähteitä: \\
Moler: Ch 11 \\
Tweito--Winther: Intro to PDE's, A computational 
approach, Springer 1998 \\ %(Jaossa osia Ch2:sta ja Ch3:sta)
Kreyszig

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

%{\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
Muodosta 1-ulotteisen Laplace-operaattorin diskretointioperaattorin $\Delta_h$
matriisiesitys, olkoon se $\frac{1}{h^2}D.$ Tässä $h=\frac{1}{m+1},$ jos välin 
pituus $=1$ ja jakopisteitä on reunat mukaan lukien $m+2$ kpl 
(aivan kuten Molerissa).

Muodosta matriisi $D$ {\tt spdiags}-komennolla, kokeile
jollain pienellä n:llä.
(Muistahan, että {\tt full} näyttää matriisin täytenä.)

Muodosta neliöalueen 2-ulotteisen Laplace-operaattorin $\Delta_h$
matriisiesitys $\frac{1}{h^2}A.$ Matriisin $A$ muodostamiseen on elegantteja
tapoja, eräs tyylikkäimmistä pohjautuu ilmaisuun

{\tt kron(D,I)+kron(I,D)}, missä {\tt I} on yksikkömatriisi.

Selvitä huolellisesti, miksi.

Muodosta matriisi A myös suoraan {\tt spdiags}:lla. Tarvitset yhtenä 
sivulävistäjänä ykkösvektorin, jossa joka n:s (sopiva n) alkio nollataan.
Tämä on kätevää tehdä vaikkapa tyylin \verb_ ~mod(v,n)==0_, missä
{\tt v} on indeksivektori {\tt v=1:pituus}. (Varo matoa, jos käytät 
desktoppia Unixissa!)

\item %2
Moler 11.2. Käytä (b)-kohtaan tietysti mieluummin Maplea suoraan.

Jätetään tämä väliin, jos joku on ehtinyt tehdä, annetaan plussaa:\\
(c) Muodosta tehtävälle Greenin funktio, määrittele se Maple-funktioksi
ja kirjoita ratkaisu sen avulla.


\item %3
Tutki ja täydennä seuraavaa Matlab-koodia.

\begin{verbatim}
n=4;
h = ...
[x,y] = meshgrid(-3:h:3);
f = peaks(x,y);
contourf(x,y,f)

e = ones(n,1);
D = spdiags(...);
I = eye(n,n);
A = (1/h^2)*...
j = 2:n+1;
F = f(j,j);
U = A\F(:);
u = zeros(n+2,n+2);
u(j,j) = reshape(U,n,n);
\end{verbatim}


Rakenna näistä aineksista yleinen Poissonin yhtälön ratkaisja
neliöalueessa. Riittää ottaa $0$-reunaehdot, mikä tekee tehtävän
varsin helpoksi. (Yleiset (Dirichlet'n) reunaehdot tulevat mukaan
tehtävässä 5.) Ratkaise sitten tehtävä 11.5 (a).



\item %4
Selvittele luentokalvojen {\tt ../L/differenssikalvoja.pdf}
avulla Poissonin yhtälön oikeaoppista Matlab-tyylistä
ratkaisutapaa.
Voit tehdä koneen ääressä kokeiluja {..L/PDE.m}-tiedostoa
mukaillen.

Pohjana on Molerin koodi funktiolle {\tt delsq}.

Selvittele koodin toiminta ja muokkaa siitä yleinen
Poissonin yhtälön ratkaisija alueessa, jonka geometrian
määrää indeksimatriisi {\tt G}, joka saadaan muutamalle
tyyppigeometrialle {\tt numgrid}-funktiolle, mutta yleisemmin
käyttämällä {\tt inregion}-funktiota
\footnote{{\tt inregion} löytyy: {\tt ../matlab/moler/:sta.}}
 [Moler-odes]-prujun
s. 6 mukaisesti.

Huomaa, että koodissa olevan for-silmukan sisällä tarvitsee
vain valita reunasolmujen indeksit ja päivittää niiden avulla
b-vektoria (joka alustetaan loopin ulkopuolella ``lähdefunktion'' f arvoilla.)

Sovella muutamaan erityyppiseen geometriaan ja lähdetrmiin
sekä reunaehtoihin.

\item %5

11.6. ``Method of lines''. Opettavaisempaa olisi aloittaa (c)-
kohdasta, mutta tee missä järjestyksessä haluat.


\subsubsection*{Ohjeita}


{\bf Teht. 3.}

Jos teet peaks-funktiolle, voisi alku olla tällainen:
\begin{verbatim}
function u=peaks_poisson(n);
% 2-D Poissonin yhtälö, lähde :  peaks(x,y), reuna-arvot = 0.

if nargin < 1, n = 31; end
h = 6/(n+1);
\end{verbatim}

Lopussa voisi piirrellä vaikka näin:
\begin{verbatim}
subplot(1,2,1)
contourf(x,y,f,-8:1:8)
caxis([-8 8])
axis square
colorbar
set(gca,'ytick',get(gca,'xtick'))

subplot(1,2,2)
a = min(min(u));
b = max(max(u));
c = a + (b-a)*(0:.1:10)';
contourf(x,y,u,c)
caxis([a b])
axis square
colorbar
set(gca,'ytick',get(gca,'xtick'))

figure(2)
subplot(1,2,1)
mesh(x,y,f)

subplot(1,2,2)
mesh(x,y,u)

\end{verbatim}

\textbf{Teht. 4}

Tässä on yleisen yhdesti yhtenäisen, polygonin rajoittaman alueen
käsittelyä selventävä
skripti (Molerin prujun ja help inregion:n yhdistelmä.)

\begin{verbatim}
clear ; close all
xv=[0 0 1 1 -1 -1 0]; yv=[0 -1 -1 1 1 0 0]
fill(xv,yv,'y')
n=4; h=1/(n+1);
[x,y]=meshgrid(-1:h:1);
%help inregion
[in,on] = inregion(x,y,xv,yv);
hold on
p = find(in-on);   % Sisäsolmujen indeksit (jonouttettuina)
q = find(on);      % Reunasolmujen indeksit
plot(xv,yv,'-',x(p),y(p),'ko',x(q),y(q),'ro')
\end{verbatim}
 
 
\end{enumerate}

\end{document}

