\input{kapittel}

\kapittel{10}{Metoder for partielle differensiallikninger}

\section*{\Large Numerisk derivasjon}

Taylorrekker kan brukes til å utlede en stor klasse av numeriske metoder for beregne tilnærminger til den deriverte.



I M1 definerte vi den deriverte til en funksjon $f$ som 
\[
f'(x)=\lim_{h \to 0} \frac{f(x+h)-f(x)}{h}.
\]
Uttrykket 
\[
\frac{f(x+h)-f(x)}{h}
\]
er stigningstallet til sekanten til $f$ mellom punktene $x$ og $x+h$. 
For små $h$ er denne sekanten en grei tilnærming til stigningstallet $f(x)$:
\[
f'(x)\approx \frac{f(x+h)-f(x)}{h}.
\]
På papiret er det slik at jo mindre $h$, desto bedre tilnærming. 
\begin{ex}
La $f(x)=e^x$, og la $h=0.1$. Vi beregner
\[
f'(1.5)\approx \frac{e^{1.6}-e^{1.5}}{0.1}=4.7134. 
\]
Merk at 
\[
f'(1.5)=e^{1.5}=4.4817,
\]
så denne tilnærmingen bommer med rundt $2 \cdot 10^{-1}$. 
Vi kan også prøve $h=0.01$. 
Da får vi
\[
f'(1.5)\approx \frac{e^{1.51}-e^{1.5}}{0.01}=4.5042. 
\]
som er noe bedre, nå er feilen på rundt $2 \cdot 10^{-2}$. Vi knekker til med enda en:
\[
f'(1.5)\approx \frac{e^{1.501}-e^{1.5}}{0.001}=4.4839,
\]
og får en feil på rundt $2 \cdot 10^{-3}$. 
\end{ex}

\begin{merkx}
Feilen i forrige eksempel er tydelig proporsjonal med $h$ - 
deler du $h$ på 10, deler du feilen på 10. 
En god illustrasjon av lineær feil. 
\end{merkx}


Det går an å utlede formelen
\[
f'(x)\approx \frac{f(x+h)-f(x)}{h}.
\]
ved å bruke taylorutviklingen til $f$ i punktet $x$:
\[
f(x+h)=f(x)+f'(x)h+\frac{f''(x)}{2}h^2+\cdots.
\]
Stigningen til sekanten kan skrives
\begin{equation*}
\frac{f(x+h)-f(x)}{h}=f'(x)+\frac{f''(x)}{2}h+\cdots
\end{equation*}
og vi ser at dette  stigninstallet består av den eksakte verdien
for $f'(x)$ pluss resten av taylorrekken til $f$
\begin{equation*}
\frac{f''(x)}{2}h+\frac{f'''(x)}{6}h^2+\cdots
\end{equation*}
Denne halen forteller oss noe om feilen.
Dersom  $h$ er liten nok, vil $h$ være mye større enn $h^2$, og vi skriver
\begin{equation*}
\frac{f''(x)}{2}h+\frac{f'''(x)}{6}h^2+\cdots=O(h)
\end{equation*}
for å signalisere at feilen er proporsjonal med $h$.



Vi kan relativt lett forbedre den lineære tilnærmingen fra forrige avsnitt ved å skrive
\begin{equation*}
f(x-h)=f(x)-f'(x)h+\frac{f''(x)}{2}h^2+\cdots
\end{equation*}
 og sette opp tilnærmingen
\begin{equation*}
\frac{f(x+h)-f(x-h)}{2h}=f'(x)+\frac{f'''(x)}{3!}h^2+\frac{f^5(x)}{5!}h^4+\cdots
\end{equation*} 
Feilen i denne tilnærmingen er 
\begin{equation*}
\frac{f'''(x)}{3!}h^2+\frac{f^5(x)}{5!}h^4+\cdots
\end{equation*}
Hvis $h$ er liten, er det rimelig å anta denne feilen er mye mindre enn for den første tilnærmingen, siden $h^2<<h$. 
\begin{ex}
La $f(x)=e^x$, og la $h=0.1$. Vi beregner
\[
f'(1.5)\approx \frac{e^{1.6}-e^{1.4}}{0.2}=4.489162287752202,
\]
som gir en feil på rundt $7.5\cdot 10^{-3}$. 
Mye bedre enn i sted. Vi kan prøver også $h=0.01$:
\[
f'(1.5)\approx \frac{e^{1.51}-e^{1.49}}{0.02}=4.481763765529401
\]
som er noe bedre, nå er feilen på rundt $7.5 \cdot 10^{-5}$. 
Knekker vi til med $h=0.001$, får vi 
\[
f'(1.5)\approx \frac{e^{1.501}-e^{1.499}}{0.002}= 4.481689817286139
\]
som gir en feil på $7.5 \cdot 10^{-7}$.
\end{ex}
\begin{merkx}
Legg merke til hvordan feilen i forrige eksempel er nærmest perfekt kvadratisk - 
vi får to nye desimaler hver gang vi deler $h$ på 10. Feilen deles altså på 100 når $h$ deles på 10.
\end{merkx}
Hvis du virkelig vil slå på stortrommen, kan du bruke formelen
\begin{align*}
f'(x)&=\\ &\frac{f(x-2h)-8f(x-h)+8f(x+h)-f(x+2h)}{12h}
\end{align*}
\begin{ex}
Vi bruker den store formelen med $h=0.1$,  $h=0.01$ og   $h=0.001$. 
Då får vi feil på $10^{-5}$, $10^{-9}$, og $10^{-13}$. 
Fire nye desimaler hver gang $h$ deles på 10. Prøv selv. 
\end{ex}

%\section*{Konsistens}


Dersom du trenger en tilnærming for $f''(x)$, 
kan du bruke den andre ordens sentraldifferansen
\begin{equation*}
\frac{f(x+h)-2f(x)+f(x-h)}{h^2}=f''(x)+O(h^2)%\frac{f''''(x)}{12}h^2+\cdots
\end{equation*}
Nå bør det etterhvert være klart hvordan en slik derivasjonsformel konstrueres - 
man søker en lineærkombinasjon av funksjonsverdier $f(x)$, ~$f(x-h)$ og  ~$f(x+h)$ og lignende ledd,
for å oppnå to ting:
\begin{itemize}
\item Korrekt tilnærming av den $n$-te deriverte.
\item Så høy orden som mulig.
\end{itemize}



Det går an å kombinere tilnærminger av forskjellig orden til å oppnå høyere ordens tilnærminger.
Vi definerer 
\begin{equation}\notag
\phi(h)= \frac{f(x+h)-f(x-h)}{2h}.
\end{equation}
Da har vi at 
\begin{equation}\notag
\phi(h)= f'(x)+h^{2}\frac{f'''(x)}{6}+h^{4}\frac{f^{5}(x)}{120}+\cdots
\end{equation}
 og  
\begin{equation}\notag
\phi\left(\frac{h}{2}\right)= f'(x)+\left(\frac{h}{2}\right)^{2}\frac{f'''(x)}{6}+\left(\frac{h}{2}\right)^{4}\frac{f^{5}(x)}{120}+\cdots
\end{equation}
 Her er trikset:
\begin{equation}\notag
\frac{4\phi\left(\frac{h}{2}\right)-\phi(h)}{3}=f'(x)-h^{4}\frac{f^{5}(x)}{480}+\cdots
\end{equation}
Med andre ord: den rette line{\ae}rkombinasjonen av to estimater med forskjellige gitterfinheter kan f{\aa} et (eller flere) ledd i feilutviklingen til {\aa} forsvinne, og da f{\aa}r vi en h{\o}yere ordens tiln{\ae}rming.
\begin{ex}
Hvis vi setter $h=0.1$, og tar to tidligere approksimasjoner, 
\[
\phi(0.1) = \frac{e^{1.6}-e^{1.4}}{0.2}=4.489162287752202,
\]
med en feil på rundt $7.5\cdot 10^{-3}$, og 
\[
\phi(0.01) = \frac{e^{1.51}-e^{1.49}}{0.02}=4.481763765529401
\]
med en feil på rundt $7.5\cdot 10^{-5}$, og 
og kombinerer dem, får vi
\begin{equation}\notag
\frac{100\phi\left(\frac{h}{10}\right)-\phi(h)}{99}=4.481689032981695
\end{equation}
som gir en feil på $-3.73 \cdot 10^{-8} \approx \left(\frac{h}{10}\right)^4$.
\end{ex}

Det går an å formulere presise teoremer som forteller hvordan man skal lineærkombinere 
tilnærminger til høyere ordens tilnærminger, 
men vi nøyer oss med denne lille smakebiten. 








Vi skal lage numerisk metoder for varmeligningen
\begin{equation*}
u_t=u_{xx}
\end{equation*}
med randkrav
\begin{equation*}
u(0,t)=u(1,t)=0
\end{equation*}
og 
\begin{equation*}
u(x,0)=f(x).
\end{equation*}
Metodene vi skal lage er basert på ting vi har gjort tidligere, 
og varmelikningen er kun en illustrasjon. 
Det er relativt lett å konstruere liknende metoder for andre likninger
med de teknikkene vi skal gå gjennom for varmelikningen.


\section*{Gitteret}
Når vi skal løse en partiell differensiallikning, 
må vi holde styr på to gitre - et i $x$-retningen, og et i $t$-retningen.
Vi skal finne approksimasjoner til løsningen $u$ i alle gitterpunktene.
Vi gitrer opp intervallet $[0,1]$ på $x$-aksen med gitteravstanden $h$, 
og nummererer punktene slik at $x_0=0$ og $x_n=1$.
Den positive $t$-aksen gitrer vi opp med gitteravstanden $k$, 
og nummerer slik at $t_0=0$.
\begin{center}
\hspace*{-3em}
    \begin{tikzpicture}[thick]
%      \draw[->] (-1.5,0) -- (1.5,0) node[right]{$x$};
          \draw[->] (-3,-1) -- (4,-1) node[right]{$x$};
          \draw[->] (-5*.6,-1) -- (-5*.6,3) node[right]{$t$};
           \draw[-] (5*.6,-1) -- (5*.6,3) ;

%          \draw[->] (-2.5,-1) -- (-2.5,1.3) node[above]{$y$};
%          \draw[line width = 1mm] (-2.35,.05) -- (-2.35,-.05) node[below]{$a$};
          \draw[line width = 1mm] (5*.6,-.95) -- (5*.6,-1.05) node[below]{$x=1$};
           \draw[line width = 1mm] (4*.6,-.95) -- (4*.6,-1.05) ;
           \draw[line width = 1mm] (3*.6,-.95) -- (3*.6,-1.05) ;
           \draw[line width = 1mm] (2*.6,-.95) -- (2*.6,-1.05) ;
           \draw[line width = 1.4mm] (1*.6,-.93) -- (1*.6,-1.07) node[below]{$x_i=ih$};
           \draw[line width = 1.4mm] (-5*.586,1) -- (-5*.614,1) node[left]{$t_j=jk$};
           \draw[line width = 1.4mm] (1*.586,1) -- (1*.614,1) node[above]{$(x_i,t_j)=(ih,jk)$};
           \draw[line width = 1mm] (0,-.95) -- (0,-1.05) ;
           \draw[line width = 1mm] (-1*.6,-.95) -- (-1*.6,-1.05) ;
           \draw[line width = 1mm] (-2*.6,-.95) -- (-2*.6,-1.05) ;
           \draw[line width = 1mm] (-3*.6,-.95) -- (-3*.6,-1.05) ;
           \draw[line width = 1mm] (-4*.6,-.95) -- (-4*.6,-1.05) ;        
           \draw[line width = 1mm] (-5*.6,-.95) -- (-5*.6,-1.05) ;
          \draw[line width = 1mm] (5*.6,.05) -- (5*.6,-.05) ;
           \draw[line width = 1mm] (4*.6,.05) -- (4*.6,-.05) ;
           \draw[line width = 1mm] (3*.6,.05) -- (3*.6,-.05) ;
           \draw[line width = 1mm] (2*.6,.05) -- (2*.6,-.05) ;
           \draw[line width = 1mm] (1*.6,.05) -- (1*.6,-.05) ;
           \draw[line width = 1mm] (0,.05) -- (0,-.05) ;
           \draw[line width = 1mm] (-1*.6,.05) -- (-1*.6,-.05) ;
           \draw[line width = 1mm] (-2*.6,.05) -- (-2*.6,-.05) ;
           \draw[line width = 1mm] (-3*.6,.05) -- (-3*.6,-.05) ;
           \draw[line width = 1mm] (-4*.6,.05) -- (-4*.6,-.05) ;        
           \draw[line width = 1mm] (-5*.6,.05) -- (-5*.6,-.05) ;

          \draw[line width = 1mm] (5*.6,.95) -- (5*.6,1.05) ;
           \draw[line width = 1mm] (4*.6,.95) -- (4*.6,1.05) ;
           \draw[line width = 1mm] (3*.6,.95) -- (3*.6,1.05) ;
           \draw[line width = 1mm] (2*.6,.95) -- (2*.6,1.05) ;
           \draw[line width = 1.4mm] (1*.6,.93) -- (1*.6,1.07) ;
           \draw[line width = 1mm] (0,.95) -- (0,1.05) ;
           \draw[line width = 1mm] (-1*.6,.95) -- (-1*.6,1.05) ;
           \draw[line width = 1mm] (-2*.6,.95) -- (-2*.6,1.05) ;
           \draw[line width = 1mm] (-3*.6,.95) -- (-3*.6,1.05) ;
           \draw[line width = 1mm] (-4*.6,.95) -- (-4*.6,1.05) ;        
           \draw[line width = 1mm] (-5*.6,.95) -- (-5*.6,1.05) ;

          \draw[line width = 1mm] (5*.6,1.95) -- (5*.6,2.05) ;
           \draw[line width = 1mm] (4*.6,1.95) -- (4*.6,2.05) ;
           \draw[line width = 1mm] (3*.6,1.95) -- (3*.6,2.05) ;
           \draw[line width = 1mm] (2*.6,1.95) -- (2*.6,2.05) ;
           \draw[line width = 1mm] (1*.6,1.95) -- (1*.6,2.05) ;
           \draw[line width = 1mm] (0,1.95) -- (0,2.05) ;
           \draw[line width = 1mm] (-1*.6,1.95) -- (-1*.6,2.05) ;
           \draw[line width = 1mm] (-2*.6,1.95) -- (-2*.6,2.05) ;
           \draw[line width = 1mm] (-3*.6,1.95) -- (-3*.6,2.05) ;
           \draw[line width = 1mm] (-4*.6,1.95) -- (-4*.6,2.05) ;        
           \draw[line width = 1mm] (-5*.6,1.95) -- (-5*.6,2.05) ;


            \draw[line width = 1mm] (-5*.59,0) -- (-5*.61,0) ;
            \draw[line width = 1mm] (-5*.59,2) -- (-5*.61,2) ;
            \end{tikzpicture}
  \end{center}

\section*{Diskretisering i $x$}
Vi setter først opp en tilnærming for $u_{xx}(x,t)$, 
basert på den andre ordens endelige differanseformelen for $x$:
\begin{equation*}
u_{xx}(x,t) \approx \frac{u(x+h,t)-2u(x,t)+u(x-h,t)}{h^2}
\end{equation*}
Man kan fint bruke høyere ordens differanseformler, 
men det skal ikke vi gjøre. Vi kan nå tenke at vi erstatter 
$u(x,t)$ med $n+1$ envariable funksjoner $u_i(t)$, 
som beskriver temperaturendringen i hvert sitt punkt $x_i$ på stangen:
\begin{equation*}
u_{xx}(x_i,t) \approx \frac{u_{i+1}(t)-2u_{ij}(t)+u_{i-1}(t)}{h^2},
\end{equation*}
Nå setter vi inn dette uttrykket i varmelikningen: 
\begin{equation*}
u_t(x_i,t)=u_{xx}(x_i,t) \approx \frac{u_{i+1}(t)-2u_{ij}(t)+u_{i-1}(t)}{h^2},
\end{equation*}
Siden $u_i$ kun er en funksjon av $t$, passer det å skrive 
\begin{equation*}
u_i'(t) = \frac{u_{i+1}(t)-2u_{ij}(t)+u_{i-1}(t)}{h^2},
\end{equation*}
som vi kjenner igjen som et system av ordinære differensiallikninger.
Dette systemet har $n-1$ likninger, 
for i $u_0=0$ og $u_n=0$ kjenner vi; 
disse er gitt av randbetingelsene.

\section*{Diskretisering i $t$}
Nå kan vi løse systemet med en ønsket metode fra forrige kapittel, 
og vi skal ta i betraktning tre valg. Etter diskretiseringen i $t$, 
skriver vi approksimasjonen i punktet $(x_i,t_j)$ som $u_{ij}$:
\[
u(x_i,t_j)\approx u_{ij}
\]
De tre metodene vi skal bruke, er eksplisitt Euler, implisitt Euler, og trapesmetoden. 
De korresponderende skjemaene for varmelikningen blir: \\

\noindent  Eksplisitt:
\[
\frac{u_{i,j+1}-u_{ij}}{k} =\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^2}  \qquad \quad 
\]

\noindent  Implisitt:
\[
\frac{u_{i,j+1}-u_{ij}}{k} =\frac{u_{i+1,j+1}-2u_{i,j+1}+u_{i-1,j+1}}{h^2} 
\]

\noindent Crank-Nicholson:
\begin{align*}
\frac{u_{i,j+1}-u_{ij}}{k} &= \frac{u_{i+1,j}-2u_{i,j}  + u_{i-1,j}}{2h^2} \\[3pt] &+\frac{u_{i+1,j+1}-2u_{i,j+1}+u_{i-1,j+1}}{2h^2}
\end{align*}

%\begin{align*}
% \text{Eksplisitt:}  &\\
%\frac{u_{i,j+1}-u_{ij}}{k} &=\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^2}  \\[5pt]
%\noindent  \text{Implisitt:}\\
%\frac{u_{i,j+1}-u_{ij}}{k} &=\frac{u_{i+1,j+1}-2u_{i,j+1}+u_{i-1,j+1}}{h^2} \\[5pt]
%\noindent \text{Crank-Nicholson:} \\
%\frac{u_{i,j+1}-u_{ij}}{k} &= \frac{u_{i+1,j}-2u_{i,j}  + u_{i-1,j}}{2h^2} \\[3pt] 
%&+\frac{u_{i+1,j+1}-2u_{i,j+1}+u_{i-1,j+1}}{2h^2}
%\end{align*}

\begin{thm}
Det explisitte og implisitte skjema er av orden 1 i $t$, og av orden 2 i $x$. 
Crank-Nicolson er av orden 2 i begge variable.
\end{thm}

\noindent Beviset er en ganske hårete variant av argumentet for ordenen til eksplisitt Euler i forrige kapittel.
Hver av disse skjemaene har sine fordeler, 
og vi skal behandle dem i tur og orden.

\subsubsection*{Eksplisitt skjema}
Det eksplisitte skjemaet 
\[
\frac{u_{i,j+1}-u_{ij}}{k} =\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^2}  
\]
har egentlig bare en fordel: 
det er så lett å programmere opp.
Vi løser for $u_{i,j+1}$, og får
\[
u_{i,j+1} = u_{ij} + \frac{k}{h^2}\left( u_{i+1,j}-2u_{i,j}+u_{i-1,j}\right)
\]
Nå er det vanlig å tegne opp noe som kalles \defterm{stensilen}.
  Dette er en figur som illustrerer hvilke gitterpunkter som er involvert i ligningen man bruker for å beregne nye approksimasjoner. 
\begin{center}
    \begin{tikzpicture}[thick]
%      \draw[->] (-1.5,0) -- (1.5,0) node[right]{$x$};
          \draw[->] (-3,-1) -- (4,-1) node[right]{$x$};
          \draw[->] (-5*.6,-1) -- (-5*.6,3) node[right]{$t$};
           \draw[->] (5*.6,-1) -- (5*.6,3) ;
%          \draw[->] (-2.5,-1) -- (-2.5,1.3) node[above]{$y$};
%          \draw[line width = 1mm] (-2.35,.05) -- (-2.35,-.05) node[below]{$a$};
          \draw[line width = 1mm] (5*.6,-.95) -- (5*.6,-1.05) node[below]{$x=1$};
           \draw[line width = 1mm] (4*.6,-.95) -- (4*.6,-1.05) ;
           \draw[line width = 1mm] (3*.6,-.95) -- (3*.6,-1.05) ;
           \draw[line width = 1mm] (2*.6,-.95) -- (2*.6,-1.05) ;
           \draw[line width = 1mm] (1*.6,-.95) -- (1*.6,-1.05) ;
           \draw[line width = 1mm] (0,-.95) -- (0,-1.05) ;
           \draw[line width = 1mm] (-1*.6,-.95) -- (-1*.6,-1.05) node[below]{$x_i$} ;
           \draw[line width = 1mm] (-2*.6,-.95) -- (-2*.6,-1.05) ;
           \draw[line width = 1mm] (-3*.6,-.95) -- (-3*.6,-1.05) ;
           \draw[line width = 1mm] (-4*.6,-.95) -- (-4*.6,-1.05) ;        
           \draw[line width = 1mm] (-5*.6,-.95) -- (-5*.6,-1.05) ;

           \draw[line width = 1mm] (0,.95) -- (0,1.05) ;
           \draw[line width = 1mm] (-1*.6,.95) -- (-1*.6,1.05) ;
           \draw[line width = 1mm] (-2*.6,.95) -- (-2*.6,1.05) ;

%           \draw[line width = 1mm] (0,1.95) -- (0,2.05) ;
           \draw[line width = 1mm] (-1*.6,1.95) -- (-1*.6,2.05) ;
 %          \draw[line width = 1mm] (-2*.6,1.95) -- (-2*.6,2.05) ;


            \draw[line width = 1mm] (-5*.59,0) -- (-5*.61,0) ;
            \draw[line width = 1mm] (-5*.59,1) -- (-5*.61,1) node[left]{$t_j$};
            \draw[line width = 1mm] (-5*.59,2) -- (-5*.61,2) ;
            \end{tikzpicture}
  \end{center}
Anta at du har beregnet $u_{ij}$ for alle $i$ og opp til og med en bestemt $j$. 
Vi kan da enkelt beregne $u_{i,j+1}$ for alle $i$ ved å bruke formelen. 

For å analysere mer, må vi sette opp et matrise-vektorprodukt som beskriver iterasjonen. 
Vi definerer 
\[
\V u_j=
\begin{bmatrix}
u_{1j} \\
u_{2j} \\
\vdots \\
u_{n-2,j} \\
u_{n-1,j}
\end{bmatrix}
\]
Det eksplisitte skjemaet kan skrives kompakt som 
\[
\V u_{j+1}= (I-\frac{k}{h^2} A) \V u_j
\]
der
\[
A=
\begin{bmatrix}
2      & -1      &  &  &   &  \\
-1      & 2      & -1 &   &  & \\[-2pt]
    & -1      & \ddots    & \ddots  & &  \\[-2pt]
    &  &\ddots & \ddots &-1  &    \\
      &   &    & -1   & 2  & -1  \\
      &    &   &    &  -1 & 2 
\end{bmatrix} 
\]
Det kan vises at denne matrisen har $n-1$ forskjellige egenverdier gitt ved
\[
\lambda_k= 4\sin^2 \frac{k\pi}{2n} \quad \text{for} \quad  1\leq k \leq n-1.
\]
Egenverdiene til 
\[
I-\frac{k}{h^2} A
\]
 er dermed gitt ved 
 \[
 1-4\frac{k}{h^2}\sin^2 \frac{k\pi}{2n}.
 \]
Dersom 
\[
\frac{k}{h^2}<\frac{1}{2}.
\]
går dette veldig bra, for alle egenverdiene er mellom 0 og 1. Da vil 
\[
\V u_{j+1}= (I-\frac{k}{h^2} A) \V u_j
\]
være en kontraksjon, og konvergere mot et fikspunkt. 
Dersom 
\[
\frac{k}{h^2}\geq \frac{1}{2}
\]
vil egenverdiene sende $\V u_i$ mot uendelig, 
og dette er ikke oppførsel vi ønsker fra et numerisk skjema for varmelikningen. 
Vi vet jo at temperaturen skal synke mot 0 i denne situasjonen!

Dette er ganske restriktivt. Kanskje noen andre skjema er enklere å ha med å gjøre.


\subsubsection*{Implisitt skjema}
Problemet med det eksplisitte skjemaet er at det går til h... med mindre $k/h^2<1/2$, og dette gjør at du må ha veldig tett mellom punktene $t$-aksen. 
Akkurat som for ordinære differensialliknigner, kan vi bøte på dette med å bruke implisitt skjema. 
Det implisitte skjemaet skriver vi
  \begin{equation*}
\left(1+2\frac{k}{h^2}\right)u_{i,j+1}-\frac{k}{h^2}u_{i+1,j+1}-\frac{k}{h^2}u_{i-1,j+1}=u_{ij}.
\end{equation*}
Stensilen er
\begin{center}
    \begin{tikzpicture}[thick]
%      \draw[->] (-1.5,0) -- (1.5,0) node[right]{$x$};
          \draw[->] (-3,-1) -- (4,-1) node[right]{$x$};
          \draw[->] (-5*.6,-1) -- (-5*.6,3) node[right]{$t$};
           \draw[->] (5*.6,-1) -- (5*.6,3) ;
%          \draw[->] (-2.5,-1) -- (-2.5,1.3) node[above]{$y$};
%          \draw[line width = 1mm] (-2.35,.05) -- (-2.35,-.05) node[below]{$a$};
          \draw[line width = 1mm] (5*.6,-.95) -- (5*.6,-1.05) node[below]{$x=1$};
           \draw[line width = 1mm] (4*.6,-.95) -- (4*.6,-1.05) ;
           \draw[line width = 1mm] (3*.6,-.95) -- (3*.6,-1.05) ;
           \draw[line width = 1mm] (2*.6,-.95) -- (2*.6,-1.05) ;
           \draw[line width = 1mm] (1*.6,-.95) -- (1*.6,-1.05) ;
           \draw[line width = 1mm] (0,-.95) -- (0,-1.05) ;
           \draw[line width = 1mm] (-1*.6,-.95) -- (-1*.6,-1.05) node[below]{$x_i$} ;
           \draw[line width = 1mm] (-2*.6,-.95) -- (-2*.6,-1.05) ;
           \draw[line width = 1mm] (-3*.6,-.95) -- (-3*.6,-1.05) ;
           \draw[line width = 1mm] (-4*.6,-.95) -- (-4*.6,-1.05) ;        
           \draw[line width = 1mm] (-5*.6,-.95) -- (-5*.6,-1.05) ;

  %         \draw[line width = 1mm] (0,.95) -- (0,1.05) ;
           \draw[line width = 1mm] (-1*.6,.95) -- (-1*.6,1.05) ;
%           \draw[line width = 1mm] (-2*.6,.95) -- (-2*.6,1.05) ;

           \draw[line width = 1mm] (0,1.95) -- (0,2.05) ;
           \draw[line width = 1mm] (-1*.6,1.95) -- (-1*.6,2.05) ;
          \draw[line width = 1mm] (-2*.6,1.95) -- (-2*.6,2.05) ;


            \draw[line width = 1mm] (-5*.59,0) -- (-5*.61,0) ;
            \draw[line width = 1mm] (-5*.59,1) -- (-5*.61,1) node[left]{$t_j$};
            \draw[line width = 1mm] (-5*.59,2) -- (-5*.61,2) ;
            \end{tikzpicture}
  \end{center}
I hver tidssteg må vi løse et lineært likningssystem
\begin{equation*}
 (I+\frac{k}{h^2} A) \mathbf{u}_{j+1}=\mathbf{u}_{j}
\end{equation*}
der $A$ er den samme matrisen som i sted. 

Egenverdiene til $ (I+\frac{k}{h^2} A)$ er 
 \[
 1+4\frac{k}{h^2}\sin^2 \frac{k\pi}{2n},
 \]
som åpenbart er større enn 1 uansett hva $\frac{k}{h^2}$ er.
Dermed vil egenverdiene til $ (I+\frac{k}{h^2} A)^{-1}$ alltid være mindre enn 1, 
og følgelig vil skjemaet konvergere mot 0 uansett. 

%der $\mathbf{u_{j}}$ er en søylevektor med approksimasjonene på tidssteg $j$, og matrisen er gitt ved 
%\begin{equation*}
%A = \begin{bmatrix} 
%    1+2\frac{k}{h^2} & -\frac{k}{h^2} & 0 & 0 \\
%    -\frac{k}{h^2} & 1+2\frac{k}{h^2} & -\frac{k}{h^2} &0  \\
%   0 &    -\frac{k}{h^2}    &  1+2\frac{k}{h^2} & -\frac{k}{h^2} \\
%   0 &  0 &    -\frac{k}{h^2}    &  1+2\frac{k}{h^2}
%    \end{bmatrix}
%\end{equation*}
%eller 
%\[
%A=I+\frac{k}{h^2}C
%\]



\subsubsection*{Crank-Nicolson}
Crank-Nicolson skriver vi 
  \begin{align*}
2u_{i,j+1}+\frac{k}{h^2}&\left(-u_{i+1,j+1}+2u_{i,j+1}-u_{i-1,j+1}\right)=\\2u_{ij} &- \frac{k}{h^2}\left(-u_{i+1,j}+2u_{i,j}-u_{i-1,j}\right).
\end{align*}
Stensilen er 
\begin{center}
    \begin{tikzpicture}[thick]
%      \draw[->] (-1.5,0) -- (1.5,0) node[right]{$x$};
          \draw[->] (-3,-1) -- (4,-1) node[right]{$x$};
          \draw[->] (-5*.6,-1) -- (-5*.6,3) node[right]{$t$};
           \draw[->] (5*.6,-1) -- (5*.6,3) ;
%          \draw[->] (-2.5,-1) -- (-2.5,1.3) node[above]{$y$};
%          \draw[line width = 1mm] (-2.35,.05) -- (-2.35,-.05) node[below]{$a$};
          \draw[line width = 1mm] (5*.6,-.95) -- (5*.6,-1.05) node[below]{$x=1$};
           \draw[line width = 1mm] (4*.6,-.95) -- (4*.6,-1.05) ;
           \draw[line width = 1mm] (3*.6,-.95) -- (3*.6,-1.05) ;
           \draw[line width = 1mm] (2*.6,-.95) -- (2*.6,-1.05) ;
           \draw[line width = 1mm] (1*.6,-.95) -- (1*.6,-1.05) ;
           \draw[line width = 1mm] (0,-.95) -- (0,-1.05) ;
           \draw[line width = 1mm] (-1*.6,-.95) -- (-1*.6,-1.05) node[below]{$x_i$};
           \draw[line width = 1mm] (-2*.6,-.95) -- (-2*.6,-1.05) ;
           \draw[line width = 1mm] (-3*.6,-.95) -- (-3*.6,-1.05) ;
           \draw[line width = 1mm] (-4*.6,-.95) -- (-4*.6,-1.05) ;        
           \draw[line width = 1mm] (-5*.6,-.95) -- (-5*.6,-1.05) ;

          \draw[line width = 1mm] (0,.95) -- (0,1.05) ;
           \draw[line width = 1mm] (-1*.6,.95) -- (-1*.6,1.05) ;
           \draw[line width = 1mm] (-2*.6,.95) -- (-2*.6,1.05) ;

           \draw[line width = 1mm] (0,1.95) -- (0,2.05) ;
           \draw[line width = 1mm] (-1*.6,1.95) -- (-1*.6,2.05) ;
          \draw[line width = 1mm] (-2*.6,1.95) -- (-2*.6,2.05) ;


            \draw[line width = 1mm] (-5*.59,0) -- (-5*.61,0) ;
            \draw[line width = 1mm] (-5*.59,1) -- (-5*.61,1) node[left]{$t_j$};
            \draw[line width = 1mm] (-5*.59,2) -- (-5*.61,2) ;
            \end{tikzpicture}
  \end{center}
Også her blir det et lineært ligningssystem 
\begin{equation*}
 (2I+\frac{k}{h^2} A)\mathbf{u}_{j+1}= (2I-\frac{k}{h^2} A)\mathbf{u}_{j}
\end{equation*}
å løse på hvert tidssteg. %Egenverdiene til 
%\[
% (2I+\frac{k}{h^2} A)^{-1} (2I-\frac{k}{h^2} A)
% \]
% blir 
% \[
%\frac{ 2-4\frac{k}{h^2}\sin^2 \frac{k\pi}{2n}}{ 2+4\frac{k}{h^2}\sin^2 \frac{k\pi}{2n}}.
% \]
% Siden 
% \[
% -2\leq 2-4\frac{k}{h^2}\sin^2 \frac{k\pi}{2n}\leq 2
% \]
% og 
% \[
% 2+4\frac{k}{h^2}\sin^2 \frac{k\pi}{2n}\geq 2
% \]
% ser vi at 
% \[
%\frac{ 2-4\frac{k}{h^2}\sin^2 \frac{k\pi}{2n}}{ 2+4\frac{k}{h^2}\sin^2 \frac{k\pi}{2n}}\leq 1,
% \]
 Det kan vises at Crank-Nicolson er en stabil metode for varmelikningen, men det dropper vi.
 
 
 \section*{Laplaces likning}

Vi kan også sette opp et numerisk skjema for Laplaces likning
\begin{equation*}
u_{xx}+u_{yy}=0,
\end{equation*}
på kvadratet $[0,1] \times [0,1]$ med randkrav 
\begin{equation*}
u(x,0)=u(0,y)=u(1,y)=0 \qquad \text{og} \qquad u(x,1)=f(x).
\end{equation*}
Vi gitrer opp på følgende vis:
\begin{center}
\hspace*{-3em}
    \begin{tikzpicture}[thick]
%      \draw[->] (-1.5,0) -- (1.5,0) node[right]{$x$};
          \draw[->] (-3,-1) -- (3,-1) node[right]{$x$};
          \draw[->] (-5*.6,-1) -- (-5*.6,4) node[right]{$y$};
           \draw[-] (3*.6,-1) -- (3*.6,3) ;
           \draw[-] (-3,3) -- (3,3) node[right]{$y=1$};
%          \draw[->] (-2.5,-1) -- (-2.5,1.3) node[above]{$y$};
%          \draw[line width = 1mm] (-2.35,.05) -- (-2.35,-.05) node[below]{$a$};
%          \draw[line width = 1mm] (5*.6,-.95) -- (5*.6,-1.05) node[below]{$x=1$};
%           \draw[line width = 1mm] (4*.6,-.95) -- (4*.6,-1.05) ;
           \draw[line width = 1mm] (3*.6,-.95) -- (3*.6,-1.05) node[below]{$x=1$};
     %      \draw[line width = 1mm] (2*.6,-.95) -- (2*.6,-1.05) ;
           \draw[line width = 1.4mm] (1*.6,-.93) -- (1*.6,-1.07) node[below]{$x_i=ih$};
           \draw[line width = 1.4mm] (-5*.586,1) -- (-5*.614,1) node[left]{$y_j=jh$};
           \draw[line width = 1.4mm] (1*.586,1) -- (1*.614,1) node[above]{$(x_i,y_j)$};
%           \draw[line width = 1mm] (0,-.95) -- (0,-1.05) ;
           \draw[line width = 1mm] (-1*.6,-.95) -- (-1*.6,-1.05) ;
%           \draw[line width = 1mm] (-2*.6,-.95) -- (-2*.6,-1.05) ;
           \draw[line width = 1mm] (-3*.6,-.95) -- (-3*.6,-1.05) ;
   %        \draw[line width = 1mm] (-4*.6,-.95) -- (-4*.6,-1.05) ;        
           \draw[line width = 1mm] (-5*.6,-.95) -- (-5*.6,-1.05) ;
  %        \draw[line width = 1mm] (5*.6,.05) -- (5*.6,-.05) ;
%           \draw[line width = 1mm] (4*.6,.05) -- (4*.6,-.05) ;
           \draw[line width = 1mm] (3*.6,.05) -- (3*.6,-.05) ;
   %        \draw[line width = 1mm] (2*.6,.05) -- (2*.6,-.05) ;
           \draw[line width = 1mm] (1*.6,.05) -- (1*.6,-.05) ;
      %     \draw[line width = 1mm] (0,.05) -- (0,-.05) ;
           \draw[line width = 1mm] (-1*.6,.05) -- (-1*.6,-.05) ;
   %        \draw[line width = 1mm] (-2*.6,.05) -- (-2*.6,-.05) ;
           \draw[line width = 1mm] (-3*.6,.05) -- (-3*.6,-.05) ;
%           \draw[line width = 1mm] (-4*.6,.05) -- (-4*.6,-.05) ;        
           \draw[line width = 1mm] (-5*.6,.05) -- (-5*.6,-.05) ;

%          \draw[line width = 1mm] (5*.6,.95) -- (5*.6,1.05) ;
  %         \draw[line width = 1mm] (4*.6,.95) -- (4*.6,1.05) ;
           \draw[line width = 1mm] (3*.6,.95) -- (3*.6,1.05) ;
%           \draw[line width = 1mm] (2*.6,.95) -- (2*.6,1.05) ;
           \draw[line width = 1.4mm] (1*.6,.93) -- (1*.6,1.07) ;
   %        \draw[line width = 1mm] (0,.95) -- (0,1.05) ;
           \draw[line width = 1mm] (-1*.6,.95) -- (-1*.6,1.05) ;
   %        \draw[line width = 1mm] (-2*.6,.95) -- (-2*.6,1.05) ;
           \draw[line width = 1mm] (-3*.6,.95) -- (-3*.6,1.05) ;
   %        \draw[line width = 1mm] (-4*.6,.95) -- (-4*.6,1.05) ;        
           \draw[line width = 1mm] (-5*.6,.95) -- (-5*.6,1.05) ;

%          \draw[line width = 1mm] (5*.6,1.95) -- (5*.6,2.05) ;
%           \draw[line width = 1mm] (4*.6,1.95) -- (4*.6,2.05) ;
           \draw[line width = 1mm] (3*.6,1.95) -- (3*.6,2.05) ;
  %         \draw[line width = 1mm] (2*.6,1.95) -- (2*.6,2.05) ;
           \draw[line width = 1mm] (1*.6,1.95) -- (1*.6,2.05) ;
  %         \draw[line width = 1mm] (0,1.95) -- (0,2.05) ;
           \draw[line width = 1mm] (-1*.6,1.95) -- (-1*.6,2.05) ;
   %        \draw[line width = 1mm] (-2*.6,1.95) -- (-2*.6,2.05) ;
           \draw[line width = 1mm] (-3*.6,1.95) -- (-3*.6,2.05) ;
   %        \draw[line width = 1mm] (-4*.6,1.95) -- (-4*.6,2.05) ;        
           \draw[line width = 1mm] (-5*.6,1.95) -- (-5*.6,2.05) ;

          \draw[line width = 1mm] (5*.6,2.95) -- (5*.6,3.05) ;
   %        \draw[line width = 1mm] (4*.6,2.95) -- (4*.6,3.05) ;
           \draw[line width = 1mm] (3*.6,2.95) -- (3*.6,3.05) ;
   %        \draw[line width = 1mm] (2*.6,2.95) -- (2*.6,3.05) ;
           \draw[line width = 1mm] (1*.6,2.95) -- (1*.6,3.05) ;
    %       \draw[line width = 1mm] (0,2.95) -- (0,3.05) ;
           \draw[line width = 1mm] (-1*.6,2.95) -- (-1*.6,3.05) ;
     %      \draw[line width = 1mm] (-2*.6,2.95) -- (-2*.6,3.05) ;
           \draw[line width = 1mm] (-3*.6,2.95) -- (-3*.6,3.05) ;
       %    \draw[line width = 1mm] (-4*.6,2.95) -- (-4*.6,3.05) ;        
           \draw[line width = 1mm] (-5*.6,2.95) -- (-5*.6,3.05) ;


            \draw[line width = 1mm] (-5*.59,0) -- (-5*.61,0) ;
            \draw[line width = 1mm] (-5*.59,2) -- (-5*.61,2) ;
            \end{tikzpicture}
  \end{center}
og bruker den kjente og kjære differanseformelen
\begin{equation*}
u_{xx}(x,t) \approx \frac{u(x+h,t)-2u(x,t)+u(x-h,t)}{h^2}
\end{equation*}
i både $x$ og $y$, slik at vi får skjemaet 
  \[
4u_{i,j}-u_{i+1,j}-u_{i-1,j}- u_{i,j+1}-u_{i,j-1}=0.
\]

Nå er $u$ kjent på sidekantene til kvadratet, 
slik at det er kun de indre gitterpunktene som gir opphav til likninger og ukjente.
La oss si at det er $n-1$ indre gitterpunkter i hver koordinatretning, 
slik at antall likninger og ukjente blir $(n-1)^2$.
Vi organiserer $u_{ij}$ i en vektor på formen 
\[
\V u=
\begin{bmatrix}
u_{11} \\
\vdots \\
u_{n-1,1} \\
u_{1,2} \\
\vdots \\
u_{n-1,2} \\
u_{1,3} \\
\vdots \\
u_{n-1,n-1}
\end{bmatrix} 
\]
og får følgelig et likningssystem på formen 
\[
A\V u=\V f
 \]
 der 
 \[
B=\begin{bmatrix}
4      & -1      &  &  &   &  \\
-1      & 4      & -1 &   &  & \\[-2pt]
    & -1      & \ddots    & \ddots  & &  \\[-2pt]
    &  &\ddots & \ddots &-1  &    \\
      &   &    & -1   & 4  & -1  \\
    &    &   &    &  -1 & 4 
\end{bmatrix} 
 \]
og
 \begin{align*}
A=&\begin{bmatrix}
B      & 0      &  &  &   &  \\
0      & B      & 0 &   &  & \\[-2pt]
    & 0      & \ddots    & \ddots  & &  \\[-2pt]
    &  &\ddots & \ddots &0  &    \\
      &   &    & 0   & B  & 0  \\
      &    &   &    &  0 & B 
\end{bmatrix} 
\\
&\qquad +
\begin{bmatrix}
0      & \cdots      & 0  & -1  &   & 0 \\
\vdots      &    0   &  &   & \ddots & \\[-2pt]
 0   &       & \ddots    & \ddots  & & -1  \\[-2pt]
  -1  &  &\ddots & \ddots &  &  0  \\
      & \ddots  &    &   &  \ddots & \vdots  \\
0      &    &  -1 &  0  & \cdots  & 0 
\end{bmatrix} 
 \end{align*}
 Den første matrisen i uttrykket for $A$ er en såkalt blokkdiagonal matrise, 
 og den andre matrisen har to sammenhengende super- og subdiagonaler med -1 i alle komponenter, 
 der superdiagonalen begynner i rad $n$, og subdiagonalen begynner i kolonne $n$. 

Under er et plot av en løsning der  $f(x)=\sin \pi x$.
 \begin{figure}[htbp]
  \begin{center}
  \vspace{-2cm}
	\includegraphics[scale=.3]{laplace.pdf}
	\label{fig:Num1}
	\end{center}
	  \vspace{-3cm}
\end{figure}


 \section*{Varmelikningen på en plate}
 
 Hvis man først har forstått hvordan laplaces likning diskretiseres, 
 er det ikke så vanskelig å forstå hvordan man diskretiserer varmelikningen i to romlige dimensjoner.
Temperaturen i platen er gitt ved funksjonen $u(x,y,t)$, 
men nå tillater vi at temperaturen varierer med tiden, og problemet vi skal løse, er  
 \begin{equation*}
u_t=u_{xx}+u_{yy}
\end{equation*}
med randkrav
\begin{align*}
u(0,y,t)=u(1,y,t)&=u(x,0,t)=0 \\ u(x,1,t)&=f(x)
\end{align*}
og initialkrav
\begin{equation*}
u(x,y,0)=g(x,y).
\end{equation*}

Gitter i $x$ og $y$ blir det samme som for laplaces likning, 
mens $k$ er tidssteget. Vi lar telleren i tid være $m$.
Vi diskretiserer med de samme differanseformlene som før, 
og får eksplisitt skjema:
\begin{align*}
\frac{u_{i,j,m+1}-u_{i,j,m}}{k} =&\\ 
\frac{u_{i+1,j,m}-2u_{i,j,m}+u_{i-1,j,m}}{h^2}& +\frac{u_{i,j+1,m}-2u_{i,j,m}+u_{i,j-1,m}}{h^2} 
\end{align*}
implisitt skjema:
\begin{align*}
\frac{u_{i,j,m+1}-u_{i,j,m}}{k} =&\\ 
\frac{u_{i+1,j,m+1}-2u_{i,j,m+1}+u_{i-1,j,m+1}}{h^2}& +\frac{u_{i,j+1,m+1}-2u_{i,j,m+1}+u_{i,j-1,m+1}}{h^2} 
\end{align*}
og Crank-Nicolson
\begin{align*}
\frac{u_{i,j,m+1}-u_{i,j,m}}{k} =&\\ 
\frac{u_{i+1,j,m}-2u_{i,j,m}+u_{i-1,j,m}}{h^2}& +\frac{u_{i,j+1,m}-2u_{i,j,m}+u_{i,j-1,m}}{h^2} 
+&\\ 
\frac{u_{i+1,j,m+1}-2u_{i,j,m+1}+u_{i-1,j,m+1}}{h^2}& +\frac{u_{i,j+1,m+1}-2u_{i,j,m+1}+u_{i,j-1,m+1}}{h^2} 
\end{align*}
Det går an å skrive disse skjemaene om til matrise-vektor-iterasjoner
på samme vis som i det endimensjonale tilfellet. 
For eksplisitt skjema får vi
\[
\V u_{m+1}= \left(I-\frac{k}{h^2}A\right) \V u_m + \frac{k}{h^2}\V f
\]
der $A$ er matrisen fra diskretiseringen av Laplaces likning, og 
\[
\V f=
\begin{bmatrix}
0 \\
\vdots \\
0 \\
f(x_1) \\
\vdots \\
f(x_{n-1})
\end{bmatrix} 
\]
Likninger for implisitt skjema og Crank-Nicolson blir henholdsvis
\[
 \left(I+\frac{k}{h^2}A\right) \V u_{m+1}=\V u_m + \frac{k}{h^2}\V f
\]
og
\[
 \left(I+\frac{k}{2h^2}A\right) \V u_{m+1}=\left( \left(I-\frac{k}{2h^2}A\right)  \V u_m + \frac{k}{2h^2}\V f \right)+ \frac{k}{2h^2}\V f.
\]



% \begin{equation*}
%A = \begin{bmatrix} 
%    1+\frac{k}{h^2} & -\frac{k}{2h^2} & 0 & 0 \\
%    -\frac{k}{2h^2} & 1+\frac{k}{h^2} & -\frac{k}{2h^2} &0  \\
%   0 &    -\frac{k}{2h^2}    &  1+\frac{k}{h^2} & -\frac{k}{2h^2} \\
%   0 &  0 &    -\frac{k}{2h^2}    &  1+\frac{k}{h^2}
%    \end{bmatrix}
%\end{equation*}
% og 
% \begin{equation*}
%B = \begin{bmatrix} 
%    1-\frac{k}{h^2} & \frac{k}{2h^2} & 0 & 0 \\
%    \frac{k}{2h^2} & 1-\frac{k}{h^2} & \frac{k}{2h^2} &0  \\
%   0 &    \frac{k}{2h^2}    &  1-\frac{k}{h^2} & \frac{k}{2h^2} \\
%   0 &  0 &    \frac{k}{2h^2}    &  1-\frac{k}{h^2}
%    \end{bmatrix}
%\end{equation*}

\clearpage


\section*{\Huge Oppgaver}

\section*{\LARGE Med LF}


\begin{oppgave}
 \end{oppgave}

\section*{\LARGE Uten LF}

\clearpage


\section*{\Huge Løsningsforslag}

\noindent 1. 


\kapittelslutt
