\input{kapittel}

\kapittel{4}{Partielle differensiallikninger I}


Fysiske situasjoner der det er behov for mer enn en uavhengig variabel, 
beskrives gjerne av partielle differensiallikninger. 
Vi skal ta for oss tre av de mest grunnleggende likningene,
nemlig bølgelikningen, varmelikningen, og Laplaces likning.

\section*{ \huge Bølgelikningen}

Bølgelikningen er en matematisk beskrivelse av en vibrerende streng, 
eller en stående luftbølge i en orgelpipe. 
Vi skal ta for oss selve likningen,
hvor den kommer fra, 
og to forskjellige løsningsteknikker - en for et intervall på $x$-aksen, 
og en for hele $x$-aksen.

\subsubsection*{Utledning}

Vi tenker at vi har en vibrerende streng som er spent opp i $x=0$ og $x=L$. La $u(x,t)$ være en funksjon som for hvert tidspunkt $t$ og hvert punkt $x$  beskriver utslaget fra likevektslinjen, som ligger langs $x$-aksen. Strengen har konstant massetetthet $\rho \, [kg/m]$.
 \begin{center}
    \begin{tikzpicture}[thick]
%      \draw[->] (-1.5,0) -- (1.5,0) node[right]{$x$};
%      \draw[->] (0,-1.5) -- (0,1.5) node[above]{$y$};
%      \filldraw[fill=blue!75!white,fill opacity=0.55] (-0.866,0.5)
%        [snake=coil,segment aspect=0,segment amplitude=1pt] --
%        (0.866,0.5) arc(30:-210:1cm);
%      \draw (0,0) circle (1cm);
%      \node[below right] at (1,0){$1$};
%      \node[below left] at (-1,0){$-1$};
%      \node[above left] at (0,1){$1$};
%      \node[below left] at (0,-1){$-1$};
%      \node[draw,fill=white,align=center] at (0,0){$75$ \%};

      \begin{scope}[xshift=5cm,yshift=-1cm]
          \draw[->] (-2.5,0) -- (2.5,0) node[right]{$x$};
          \draw[->] (-2.5,-0.25) -- (-2.5,1.3) node[above]{$y$};
          \draw (2,0) arc(55:125:3.9cm) node[right,xshift=2cm,yshift=1.2cm]{$u(x,t)$};
%          \draw (0,0) -- (2*0.866,1);
%          \draw[snake=coil,segment aspect=0,segment amplitude=1pt]
%            (0,1) -- (2*0.866,1);
%          \draw (0,0.35) arc(90:30:0.35cm)
%            node[above,xshift=-0.05cm,yshift=0.05cm]{$\theta$};
%          \draw[snake=brace] (-0.05,0) --
%            node[left,xshift=-0.05cm]{$h$} (-0.05,1);
          \node[below] at (2,0){$L$};
%          \node[left] at (0,2){$L$};
      \end{scope}
    \end{tikzpicture}
  \end{center}
%  
%\begin{center}
%    \begin{tikzpicture}
%      % Vanntank
%      \filldraw[black,fill=blue!20] (0,0) -- (1,2) arc(0:180:1cm and
%        0.25cm) -- cycle;
%      \draw (-1,2) arc(-180:0:1cm and 0.25cm);
%      \draw (0,0) -- (1.5,3) arc(0:180:1.5cm and 0.375cm) -- cycle;
%      \draw (-1.5,3) arc(-180:0:1.5cm and 0.375cm);
%
%      % Mål
%      \draw (-0.0725,0) -- (-1.2,0);
%      \draw[|<->|] (-1.2,-0.007) --
%        node[midway,fill=white,font=\scriptsize]{$600$} (-1.2,2);
%      % \draw[densely dashed] (0,2) -- (0,2.625);
%      % \draw (0,2.7) -- (0,3.75);
%      % \draw[densely dashed] (1,2) -- (1,2.725);
%      % \draw (1,2.8) -- (1,3.75);
%      % \draw[|<->|] (-0.007,3.75) --
%      %   node[midway,fill=white,font=\scriptsize]{$120$}
%      %   (1.007,3.75);
%      \draw (0,3) -- (0,4);
%      \draw[|<->|] (-0.007,4) --
%        node[midway,fill=white,font=\scriptsize]{$160$} (1.507,4);
%      \draw (1.5,4) -- (1.5,3);
%      \draw (0.0725,0) -- (1.7,0);
%      \draw[|<->|] (1.7,-0.007) --
%        node[midway,fill=white,font=\scriptsize]{$800$} (1.7,3.007);
%
%      % Pil ned
%      \draw[->] (0,-0.05) -- (0,-0.25);
%    \end{tikzpicture}
 % \end{center}


Vi tar en nærmere titt på strekkreftene på et lite stykke av strengen. Vi antar at tyngdekraften er neglisjerbar, og at strengen er helt elastisk, slik at strengestrekket, som virker parallelt med strengen, er eneste kraft. Vi antar at hvert punkt på strengen kun beveger seg loddrett, og at den horisontale komponenten av strengestrekket er konstant lik $T$. 
 \begin{center}
    \begin{tikzpicture}[thick]
 \begin{scope}[xshift=5cm,yshift=-1cm]
          \draw[->] (-1,.63) -- (.5,.9) node[right,yshift=.3cm]{$[T,T_2]$};
          \draw[-] (-1,.05) -- (-1,-.05) node[below,xshift=.4cm, yshift=.07cm]{$x+h$};
                    
          \draw[->] (-1.5,.5) -- (-3,0) node[left]{$[T,T_1]$};
          \draw[-] (-1.5,.05) -- (-1.5,-.05) node[below]{$x$};
          
          \draw (2,0) arc(55:125:3.9cm) node[right,xshift=4cm,yshift=.5cm]{$u(x,t)$};
          \draw[->] (-2.5,0) -- (2.5,0);
%          \draw (0,0) -- (2*0.866,1);
%          \draw[snake=coil,segment aspect=0,segment amplitude=1pt]
%            (0,1) -- (2*0.866,1);
%          \draw (0,0.35) arc(90:30:0.35cm)
%            node[above,xshift=-0.05cm,yshift=0.05cm]{$\theta$};
%          \draw[snake=brace] (-0.05,0) --
%            node[left,xshift=-0.05cm]{$h$} (-0.05,1);
%          \node[below] at (2,0){$L$};
%          \node[left] at (0,2){$L$};
      \end{scope}
    \end{tikzpicture}
  \end{center}
Vi setter opp Newtons andre lov for den lille biten fra $x$ til $x+h$. Massen til en bit med lengde $h$ er $h\rho$, og akselerasjonen til strengen i punktet $x$ er $u_{tt}(x,t)$. Netto kraft på biten er gitt ved $T_2 + T_1$, slik at 
\[
h\rho u_{tt}(x,t) = T_2 + T_1,
\]
eller 
\begin{align*}
\frac{\rho}{T} u_{tt}(x,t) &=\frac{T_2/T + T_1/T}{h}\\&=\frac{u_{x}(x+h,t) - u_{x}(x,t)}{h},
\end{align*}
siden stigningstallet til tangenten til strengen er gitt ved $u_{x}$. Lar vi nå $h\to 0$, får vi bølgeligningen
\[
 u_{tt}(x,t) = c^2 u_{xx}(x,t),
\]
der $c^2=\frac{T}{\rho}$.

\subsubsection*{Oppstilling av problem}

Det er ikke nok med en differensialligning som beskriver strengens bevegelse. 
Vi må også ha informasjon om hvordan bevegelsen blir satt igang, 
og hvor strengen er spent opp. Et fullstendig oppstilt problem er:
\begin{equation}\label{bolgeligningen}
 u_{tt}(x,t) = c^2 u_{xx}(x,t),
\end{equation}
med randkrav 
\begin{equation}\label{randkrav}
u(0,t)=u(L,t)=0,
\end{equation}
og initialkrav
\begin{equation}\label{initialkrav}
u(x,0)=f(x) \quad \quad u_t(x,0)=g(x). 
\end{equation}

Selve differensialligningen (\ref{bolgeligningen}) forteller oss hva slags fysiske lover som skal tilfredsstilles (i dette tilfelle Newtons andre lov), eller hva slags oppførsel vi kan forvente av løsningen, for eksempel at det er en vibrerende streng det er snakk om. Randkravene (\ref{randkrav}) forteller oss at strengen er spent opp i $x=0$ og $x=L$, slik at løsningen står helt i ro der. Initialkravene (\ref{initialkrav}) forteller oss noen om hvordan bevegelsen settes i gang; $f$ angir strengens posisjon ved $t=0$, mens $g$ angir strengens fart ved $t=0$. Når man spiller en tone på en gitar ved å dra i strengen og slippe den, slik man vanligvis gjør, er $g=0$. Randkravene
\begin{equation}
u(0,t)=u(L,t)=0,
\end{equation}
kalles dirichletrandkrav.

\subsubsection*{Løsning på intervall - separasjon av variable}

Et stort geni har engang tenkt at løsningen på bølgeligningen kan skrives
\[
u(x,t)=F(x)G(t).
\]
Han hadde rett. Innsetting i (\ref{bolgeligningen}) gir 
\[
F(x)G''(t)=c^2F''(x)G(t).
\]
Vi deler på $c^2F(x)G(t)$ og får
\[
\frac{F''(x)}{F(x)}=\frac{G''(t)}{c^2G(t)}.
\]
Siden $x$ og $t$ skal kunne varieres uavhengig av hverandre, må vi ha at 
%\[
%\frac{F''(x)}{F(x)}=\frac{G''(t)}{c^2G(t)}=-k \quad \text{der} \; k>0,
%\]
\[
\frac{F''(x)}{F(x)}=\frac{G''(t)}{c^2G(t)}=k 
\]
der $k$ er en foreløbig ubestemt konstant. 
Vi ganger opp med $c^2F(x)G(t)$ og bytter fortegn på $k$, slik at 
\[
F''(x)+kF(x)=0
\]
og 
\[
G''(t)+kc^2G(t)=0.
\]



Vi skal først prøve å finne ut hva $k$ kan være. 
Vi kan bruke $F$ og randkravene 
\[
u(0,t)=u(L,t)=0,
\]
til dette. 
La oss kikke på
\[
F''(x)+kF(x)=0.
\]
Dersom $k=0$, får vi 
\[
F''(x)=0
\]
som gir
\[
F(x)=Ax+B.
\]
Er dette en interessant løsning?
Vel, nei. Dersom 
\[
F(0)G(t)=u(0,t)=0,
\]
må enten $F(0)=0$ eller $G(t)=0$. 
At $G(t)=0$, slik at $u(x,t)=0$, er en gyldig løsning av bølgelikningen, som også tilfredsstiller randkravene.
Men dette er åpenbart ikke en spesielt interessant løsning, 
så jeg tror vi går for $F(0)=0$, som impliserer $B=0$. 
Det andre randkravet
\[
u(L,t)=0
\]
gir likeledes at $F(L)=0$, altså at $AL=0$, 
som impliserer at $A=0$. 
Dette impliserer igjen at $u(x,t)=0$, 
vi konkluderer at $k=0$ og $F(x)=0$ ikke er en interessant løsning av problemet. 

La oss prøve $k<0$. 
I så fall løses 
\[
F''(x)+kF(x)=0
\]
 av
 \[
 F(x)=Ae^{\sqrt{-k}x}+Be^{-\sqrt{-k}x}
\]
Dersom vi nå bruker randkravene, får vi de to likningene
\begin{align*}
A+B=0 \\
Ae^{\sqrt{-k}L}+Be^{-\sqrt{-k}L}=0
\end{align*}
Dersom $k \neq 0$ (som vi jo allerede vet), 
er determinanten til dette systemet gitt ved 
\[
e^{\sqrt{-k}L}+e^{-\sqrt{-k}L} \neq 0
\]
og vi konkluderer med at $A=B=0$, slik at $F(x)=0$. 
Altså er heller ikke $k<0$ en interessant løsning. 


Dersom $k>0$, går alt så meget bedre, og vi får
\[
F(x)=A\cos \sqrt{k}x + B \sin \sqrt{k}x.
\]
Bruker vi  randkravet $u(x,t)=0$, får vi 
\[
A=0,
\]
og krever vi
\[
F(L)=B \sin (\sqrt{k}L) G(t)=0,
\]
kan dette oppnås ved å sette $B=0$, som er uinteressant siden da blir $u(x,t)=0$, eller ved å kreve
\[
\sqrt{k}L = n\pi.
\]
slik at 
\[
k=\left(\frac{n\pi}{L}\right)^2.
\]
%Nå ser vi hvorfor $k>0$. Dersom $k<0$, blir $F(x)=A e^{\sqrt{k}x} + B e^{-\sqrt{k}x}$, og da gir randkravene $F(x)=0$. Dersom $k=0$, blir $F(x)=Ax+B$ og da gir randkravene også $F(x)=0$. 
Vi ser også at $n>0$, for dersom $n<0$ byttes bare fortegnet på $B$, som ennå er ubestemt. Ved å ta en titt på det endelige løsningen av problemet nedenfor, 
ser man at $B$ kommer til  å bli overflødig, så vi velger $B=1$. 

Ligningene
\[
G''(t)+c^2\left(\frac{n\pi}{L}\right)^2G(t)=0.
\]
løses av
\[
G_n(t)=A_n\cos c\frac{n\pi}{L} t+B_n\sin c\frac{n\pi}{L} t,
\]
så de generelle løsningene til bølgeligningen med randkrav (\ref{randkrav}) blir
\begin{align*}
u_n(x,t)&=F(x)G_n(t)\\&=\left(A_n\cos c\frac{n\pi}{L} t+B_n\sin c\frac{n\pi}{L} t\right) \sin \frac{n\pi}{L}x.
\end{align*}
Vi har ennå ikke tatt stilling til initialkravene (\ref{initialkrav}). Det kan vi klare ved å skrive 
\begin{align*}
u(x,t)&=\sum_{n=1}^{\infty}u_n(x,t)\\&=\sum_{n=1}^{\infty}\left(A_n\cos c\frac{n\pi}{L} t+B_n\sin c\frac{n\pi}{L} t\right) \sin \frac{n\pi}{L}x.
\end{align*}

Dersom vi nå krever
\[
f(x)=u(x,0)=\sum_{n=1}^{\infty}u_n(x,0)=\sum_{n=1}^{\infty}A_n \sin \frac{n\pi}{L}x, 
\]
ser vi at summen til høyre bør være fourierrekken til den odde utvidelsen til $f$, og følgelig bør 
\[
A_n=\frac{1}{L} \int_{-L}^L f(x) \sin \frac{n\pi}{L}x \; dx=\frac{2}{L} \int_{0}^L f(x) \sin \frac{n\pi}{L}x \; dx.
\]
På samme vis, dersom
\[
g(x)=u_t(x,0)=\sum_{n=1}^{\infty}B_nc\frac{n\pi}{L} \sin \frac{n\pi}{L}x, 
\]
bør summen til høyre være fourierekken til den odde utvidelsen til $g$, og følgelig må
\[
B_nc\frac{n\pi}{L}=\frac{2}{L} \int_{0}^L g(x) \sin \frac{n\pi}{L}x \; dx.
\]
slik at 
\[
B_n=\frac{2}{cn\pi} \int_{0}^L g(x) \sin \frac{n\pi}{L}x \; dx.
\]

Når vi skal oppsummere, 
må vi være litt forsiktige, 
for vi har ikke sagt noe om deriverbarhet. 
En løsning av bølgelikningen bør helst være to ganger kontinuerlig deriverbar, 
men dersom du drar i en gitarstreng med fingeren og slipper strengen, 
vil $f$ har form som en sagtann, 
og ikke være deriverbar i punktet der fingeren din holder strengen. 
(Ihvertfall hvis strengen er elastisk, og fingen din er uendelig tynn.)
Av disse grunnene bør man kanskje ikke kalle oppsummeringen for et teorem, 
men det får få en rød ramme siden det ligger et visst resonnement bak.

\begin{tcolorbox}[colback=red!10, title=Løsningsformel for dirichletrandkrav]
Bølgeligningen
\[
 u_{tt}(x,t) = c^2 u_{xx}(x,t),
\]
med randkrav 
\[
u(0,t)=u(L,t)=0,
\]
og initialkrav
\[
u(x,0)=f(x) \quad \quad u_t(x,0)=g(x), 
\]
løses av 
\begin{align*}
u(x,t)=\\ \sum_{n=1}^{\infty}&\left(A_n\cos c\frac{n\pi}{L} t+B_n\sin c\frac{n\pi}{L} t\right) \sin \frac{n\pi}{L}x,
\end{align*}
der 
\[
A_n=\frac{2}{L} \int_{0}^L f(x) \sin \frac{n\pi}{L}x \; dx 
\]
og 
\[
 B_n=\frac{2}{cn\pi} \int_{0}^L g(x) \sin \frac{n\pi}{L}x \; dx.
 \]

\end{tcolorbox}

\subsubsection*{Løsning for fløyte}

En stående trykkbølge inne i fløyte, 
beskrives av problemet
\[
 u_{tt}(x,t) = c^2 u_{xx}(x,t),
\]
med randkrav 
\[
u_x(0,t)=u_x(L,t)=0,
\]
og initialkrav
\[
u(x,0)=f(x) \quad \quad u_t(x,0)=g(x). 
\]
Konstanten $c$ avhenger nå av trykket og lydhastigheten i mediet der lydbølgene produseres. 
(Dersom du spiller på fløyten inni en gassballong full av helium, blir $c$ høyere enn i luft.) 
Randkravene kalles \defterm{von-Neumann-randkrav}, og $L$ er lengden på fløyten.

Vi løser problemet på samme måte som for den vibrerende strengen. 
Den eneste forskjellen blir i forbindelse med løsning av 
\[
F''(x)+kF(x)=0.
\]
På samme vis må $k>0$, slik at 
\[
F(x)=A\cos \sqrt{k}x + B \sin \sqrt{k}x.
\]
Siden 
\[
F'(x)=-A\sqrt{k}\sin \sqrt{k}x + B \sqrt{k}\cos \sqrt{k}x
\]
impliserer 
\[
u_x(0,t)=F'(0)=0,
\]
at $B=0$, mens 
\[
u_x(L,t)=F'(x)=-A\sqrt{k}\sin \sqrt{k}L=0,
\]
gir som før at 
\[
\sqrt{k}L=n\pi.
\]
Merk at $n=0$ er en gyldig løsning her siden en konstant funksjon løser både bølgelikningen og von-Neumann-randkravene.
Resten blir som før, men vi må bruke cosinusrekkene til $f$ og $g$ istedet for sinusrekkene. Vi tar ikke resten av regningen.

\begin{tcolorbox}[colback=red!10, title=Løsningsformel for fløyte]

Bølgeligningen
\[
 u_{tt}(x,t) = c^2 u_{xx}(x,t),
\]
med randkrav 
\[
u_x(0,t)=u_x(L,t)=0,
\]
og initialkrav
\[
u(x,0)=f(x) \quad \quad u_t(x,0)=g(x), 
\]
løses av 
\begin{align*}
u(x,t)=\\ A+ \sum_{n=1}^{\infty}&\left(A_n\cos c\frac{n\pi}{L} t+B_n\sin c\frac{n\pi}{L} t\right) \cos \frac{n\pi}{L}x,
\end{align*}
der $A$ er en vilkårlig konstant,
\[
A_n=\frac{2}{L} \int_{0}^L f(x) \cos \frac{n\pi}{L}x \; dx 
\]
og 
\[
 B_n=\frac{2}{cn\pi} \int_{0}^L g(x) \cos \frac{n\pi}{L}x \; dx.
 \]

\end{tcolorbox}

\subsubsection*{Løsning på hele $x$-aksen - D'alembert}

Hvis man ikke ønsker å bruke bøglelikningen til å beskrive en oppspent streng, 
men heller bølgene fra et steinkast på et endimensjonalt og uendelig langt hav,
må man studere bølgelikningen
\[
 u_{tt}(x,t) = c^2 u_{xx}(x,t),
\]
med initialkrav
\[
u(x,0)=f(x) \quad \quad u_t(x,0)=g(x),
\]
og ingen randkrav.
Her skal løsningen gjelde for alle $x$, 
ikke bare på intervallet $[0,L]$

I noen tilfeller kan dette gjøre alt vanskeligere, 
men akkurat i tilfellet bølgelikningen, 
blir alt fryktelig enkelt.
Det er enkelt å sjekke at funksjonen
\[
u(x,t)=\phi(x+ct) + \psi(x-ct)
\]
passer i likningen 
\[
 u_{tt}(x,t) = c^2 u_{xx}(x,t).
\]
uansett hva $\phi$ og $\psi$ er, så lenge de er to ganger kontinuerlig deriverbare.
Spørsmålet blir bare da hvordan vi skal klare å innpasse initialkravene.

Dersom man bruker
\[
u(x,0)=f(x) 
\]
får man 
\[
\phi(x) + \psi(x)=f(x) 
\]
og dersom man bruker
\[
u_t(x,0)=g(x),
\]
får man 
\[
c\phi'(x) - c\psi'(x)=g(x),
\]
eller 
\[
\phi(x) - \psi(x)=\frac{1}{c}\int^x g(t)\; dt +C
\]
Vi har nå et lineært $2 \times 2$-likningssystem for $\phi$ og $\psi$.
Legger vi likningene sammen, 
får vi
\begin{align*}
2\phi(x)=f(x)+\frac{1}{c}\int^x g(t)\; dt +C
\end{align*}
og trekker vi dem fra hverandre, 
får vi
får vi
\begin{align*}
2\psi(x)=f(x)-\frac{1}{c}\int^x g(t)\; dt -C
\end{align*}
Vi setter nå alt sammen igjen
\begin{align*}
u(x,t)=&\phi(x+ct)+\psi(x-ct) \\
=&\frac{1}{2}\left(f(x+ct)+f(x-ct)\right)\\[3pt]
&+\frac{1}{2c}\left(\int^{x+ct} g(t)\; dt -\int^{x-ct} g(t)\; dt\right) \\[3pt]
=&\frac{1}{2}\left(f(x+ct)+f(x-ct)\right)+\frac{1}{2c}\int^{x+ct}_{x-ct} g(t)\; dt,
\end{align*}
og oppsummerer.

\begin{tcolorbox}[colback=red!10, title=D'Alemberts løsningsformel]

Bølgeligningen
\[
 u_{tt}(x,t) = c^2 u_{xx}(x,t),
\]
på hele $x$-aksen, med initialkrav
\[
u(x,0)=f(x) \quad \quad u_t(x,0)=g(x), 
\]
løses av 
\begin{align*}
u(&x,t)= \\
&\frac{1}{2}\left(f(x+ct)+f(x-ct)\right)+\frac{1}{2c}\int^{x+ct}_{x-ct} g(t)\; dt.
\end{align*}

\vspace{.3cm}
\end{tcolorbox}


\section*{ \huge Varmelikningen}


Varmelikningen beskriver kort og godt varmeflyt i materialer. 
Vi skal ta for oss den enkleste varianten, 
der varmeledningsevnen er identisk overalt og i alle retninger i materialet.

\begin{tcolorbox}[title=Lesestoff]
\url{https://www.feynmanlectures.caltech.edu/II_02.html}\\
\url{https://www.feynmanlectures.caltech.edu/II_03.html}
\end{tcolorbox}

\subsubsection*{Utledning}
I en ideell gass er det et lineært forhold mellom temperatur og opplagret bevegelsesenergi, 
gitt ved 
\[
\frac{1}{2}m v^2=\frac{3}{2} k T
\]
der $k=1.38064852(79)\times 10^{-23}$ J/K er Boltzmanns konstant, 
$v$ er gjennomsnittsfarten til gassmolekylene, 
og $T$ er temperaturen målt i Kelvin.
Vi kan altså tenke på temperatur som et mål 
på opplagret kinetisk energi.
I andre materialer enn ideelle gasser kan den kinetiske energien også lagres opp som rotasjon og vibrasjon.



La oss tenke at temperaturen er gitt ved funksjonen $u(\mathbf x,t)$,  
der $\mathbf x= (x_1,x_2,x_3)$ er de romlige koordinatene og $t$ er tidspunktet.
La $\Omega$ være et område i $\mathbb R^3$.
Den totale opplagrede varmeenergien i $\Omega$ er gitt ved
\[
\rho c 
\iiint_{\Omega} u(\mathbf x,t)\; d\mathbf x,
\]
der $\rho$ er massetettheten og $c$ er den spesifikke varmekapasiteten.
Dersom vi antar at $\rho$ og $c$ er konstante og legger noen restriksjoner på glattheten til $u$ (for kompliserte for oss),
er endringen i denne energien over tid gitt ved
\[
\rho c 
\frac{\partial }{\partial t}\iiint_{\Omega} u(\mathbf x,t)\; d\mathbf x = 
\rho c 
\iiint_{\Omega} u_t(\mathbf x,t)\; d\mathbf x 
\]

Transport av varmeenergi er for mange materialer proporsjonal med temperaturgradienten.
Dette kalles Fouriers lov, og skrives
\[
\V q=-\kappa \nabla u
\]
der $\kappa$ er den termiske konduktiviteten, og $\mathbf q$ er varmefluksen.  
Endringen i varmeenergi i $\Omega$ skyldes enten varme som produseres inne i $\Omega$ 
(tenk at det står et stearinlys og brenner inne i $\Omega$),  
eller varme som slipper inn og ut gjennom randen $\partial \Omega$.
Dersom vi antar at $\kappa$ er konstant og at det ikke produseres varme i $\Omega$, må 
\begin{align*}
-\rho c \iiint_{\Omega} u_t\; d\mathbf x&=\iint_{\partial \Omega} \V q \cdot dS \\ 
&=\iiint_{\Omega} \nabla \cdot \V q \; d\mathbf x\\
&=-\kappa \iiint_{\Omega} \nabla \cdot \nabla u \; d\mathbf x=-\kappa \iiint_{\Omega} \Delta u \; d\mathbf x
\end{align*}
eller 
\[
\rho c \iiint_{\Omega} u_t\; d\mathbf x=
\kappa \iiint_{\Omega} \Delta u \; d\mathbf x.
\]
Dersom  dette  gjelder for et tilfeldig valgt område $\Omega$, 
bør integrandene være like:
\[
\rho c u_t=\kappa \Delta u
\]
Det er vanlig å samle den spesifikke varmekapasiteten, massetettheten og den termiske konduktiviteten i en konstant
\[
\alpha= \frac{\kappa}{\rho c}
\]
som kalles den termiske diffusiviteten. 
Likningen 
\[
u_t=\alpha \Delta u
\]
kalles 
{varmelikningen}
eller {diffusjonslikningen}.
Det er også vanlig å sette $\alpha=1$, 
for dersom $u$ løser 
$
u_t=\alpha \Delta u
$
må $v(\mathbf x, t)=u(\mathbf x, t/\alpha)
$
løse 
$
v_t= \Delta v.
$

\subsubsection*{Løsning på begrenset intervall}


Nå skal vi tenke oss en isolert stang der temperaturen er gitt av $u(x,t)$. Varme kan slippe ut eller tilføres gjennom endepunktene, men ingen andre steder. 
Vi skal løse varmeligningen 
\[
u_t= u_{xx}
\]
på stang med lengde $L$. Randkravene
\begin{equation}
u(0,t)=u(L,t)=0,
\end{equation}
sier at temperaturen holdes konstant lik 0 i endepunktene, og initialkravet
\begin{equation}
u(x,0)=f(x)
\end{equation}
sier at temperaturfordelingen er gitt ved $f(x)$ når tiden begynner å gå. 

Vi prøver separasjon av variable
\[
u(x,t)=F(x)G(t).
\]
På samme vis som for bølgeligningen, får vi 
\[
F''(x)+kF(x)=0
\]
og 
\[
G'(t)+kG(t)=0.
\]
Å finne $F$ og $k$ blir eksakt repetisjon av argumentet for bølgeligningen, mens 
\[
G_n(t)=A_n e^{-\left(\frac{n\pi}{L}\right)^2t}
\]
slik at de generelle løsningene til bølgeligningen med randkrav $u(0,t)=u(L,t)=0$ blir
\[
u_n(x,t)=F(x)G_n(t)=A_n e^{-\left(\frac{n\pi}{L}\right)^2t} \sin \frac{n\pi}{L}x.
\]



Likeledes kan initialkravet $u(x,0)=f(x)$ inkorporeres ved å legge sammen alle løsninger 
\[
u(x,t)=\sum_{n=1}^{\infty}u_n(x,t)=\sum_{n=1}^{\infty} A_n e^{-\left(\frac{n\pi}{L}\right)^2t}\sin \frac{n\pi}{L}x,
\]
sette $t=0$
\[
f(x)=u(x,0)=\sum_{n=1}^{\infty} A_n \sin \frac{n\pi}{L}x,
\]
og observere at vi får dette til dersom $A_n$ er fourierkoeffisientene til den odde utvidelsen til $f$:
\[
A_n=\frac{2}{L} \int_{0}^L f(x) \sin \frac{n\pi}{L}x \; dx
\]




\begin{tcolorbox}[colback=red!10, title=Løsningsformel for varmeledning på gittarstreng]

Varmelikningen
\[
 u_{t}(x,t) =  u_{xx}(x,t),
\]
med randkrav 
\[
u(0,t)=u(L,t)=0,
\]
og initialkrav
\[
u(x,0)=f(x) 
\]
løses av 
\begin{align*}
u(x,t)=\sum_{n=1}^{\infty} A_n e^{-\left(\frac{n\pi}{L}\right)^2t}\sin \frac{n\pi}{L}x,
\end{align*}
der 
\[
A_n=\frac{2}{L} \int_{0}^L f(x) \sin \frac{n\pi}{L}x \; dx.
\]

\end{tcolorbox}


\subsubsection*{Løsning på hele $x$-aksen}

Vi skal nå løse varmeligningen
\[
u_t= u_{xx}
\]
på en uendelig lang stang. Randkravene skal være
\begin{equation}
\lim_{x \to \pm \infty}u(x,t)=0,
\end{equation}
mens initialkravet er som før
\begin{equation}
u(x,0)=f(x),
\end{equation}
der vi antar
\begin{equation}
\lim_{x \to \pm \infty}f(x)=0.
\end{equation}


Vi skal løse dette problemet ved å fouriertransformere varmeligningen med hensyn på $x$. 
Siden 
$\lim_{x \to  \infty}u(x,t)=\lim_{x \to -\infty}u(x,t)=0$, 
satser vi på at det går bra, og vi får
\[
\mathcal{F}(u_t)=\mathcal{F}\left(u_{xx}\right).
\]
Hvis vi antar at $u$ er kontinuerlig deriverbar i $t$, og fouriertransformen er med hensyn på $x$, kan vi skrive
\begin{align*}
\mathcal{F}(u_t)&=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}u_t(x,t)e^{-i\omega x}\; dx\\
&=\frac{\partial}{\partial t}\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}u(x,t)e^{-i\omega x}\; dx=\hat{u}_t(w,t),
\end{align*}
og 
\[
\mathcal{F}(u_{xx})=-\omega ^2\mathcal{F}(u)=-\omega ^2\hat{u}(w,t),
\]
slik at 
\[
\hat{u}_t=-\omega^2 \hat{u}.
\]

Dette er en ordinær differensialligning for $\hat{u}$, med løsning
\[
\hat{u}(\omega,t)=A(\omega)e^{-\omega^2 t}.
\]
Fouriertransformerer vi initialkavet, ser vi at 
\[
A(\omega )=\hat{u}(\omega,0)=\mathcal{F}\left(u(x,0)\right)=\mathcal{F}(f(x))=\hat{f}(\omega),
\]
og da er det egentlig bare å inversfouriertransformere, og så har vi løsningen
\[
u(x,t)=\mathcal{F}^{-1}\left(\hat{u}\right)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} \hat{f}(\omega) e^{-\omega^2t} e^{i\omega x}\;d\omega.
\]

Vel og bra, men vi kan komme oss et knepp til. Husk at 
\[
f \ast g= \mathcal{F}^{-1}\left(\sqrt{2\pi} \; \hat{f}\hat{g}\right),
\]
og observer at 
\[
\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} \hat{f}(\omega) e^{-\omega^2t} e^{i\omega x}\;d\omega=\mathcal{F}^{-1}\left( \hat{f}e^{-\omega^2t}\right).
\]
Med andre ord, hvis vi kunne inverstransformere $e^{-\omega^2t}$ hadde det vært bra. Husk at
\[
\mathcal{F}\left(e^{-ax^2}\right)=\frac{1}{\sqrt{2a}}e^{\frac{-\omega^2}{4a}}.
\]
Setter vi 
\[
a=\frac{1}{4 t},
\]
og ganger og deler litt med konstanter, kan vi si at 
\[
\mathcal{F}\left(\frac{1}{\sqrt{2 t}} e^{\frac{-x^2}{4 t}}\right)=e^{-\omega^2t},
\]
som gir
\begin{align*}
u(x,t)=&\mathcal{F}^{-1}(\hat{f}e^{-\omega^2t})\\[3pt]
=&\frac{1}{\sqrt{2\pi}}\left(f \ast \frac{1}{\sqrt{2 t}} e^{\frac{-x^2}{4 t}}\right)\\[3pt]
=&\frac{1}{2\sqrt{\pi t}}\int_{-\infty}^{\infty} f(v)e^{\frac{-(x-v)^2}{4t}}\;dv
\end{align*}
Onde tunger vil kanskje hevde at $u(x,t)$ ikke er definert for $t=0$. 
Dette er ikke bra, for vi har jo prøvd å få til at $u(x,0)=f(x)$, 
som åpenbart ikke kan være helt riktig. 
Nå kan det heldigvis vises at 
\[
\lim_{t \to 0} u(x,t)=f(x)
\]
men dette er for vanskelig for oss.


\begin{tcolorbox}[colback=red!10, title=Kjempe lang stang på Vang]

Varmelikningen
\[
 u_{t} = u_{xx},
\]
på hele $x$-aksen med
initialkrav
\[
u(x,0)=f(x) 
\]
løses av 
\begin{align*}
u(x,t)=\frac{1}{2\sqrt{\pi t}}\int_{-\infty}^{\infty} f(v)e^{\frac{-(x-v)^2}{4t}}\;dv,
\end{align*}
og vi har 
\[
\lim_{t \to 0} u(x,t)=f(x).
\]

\vspace{.3cm}
\end{tcolorbox}


\section*{ \huge Laplaces ligning}

For å få en følelse for Laplaces ligning kommer fra, må man betrakte varmeligningen i to romlige dimensjoner
\[
u_t=u_{xx}+u_{yy}.
\]
Det fysiske bildet du kan ha, er en tynn varmeisolert plate der varmeenergi kan forsvinne ut eller inn gjennom sidekantene. 
Løsningen til varmeligningen vil da gi temperaturen $u(x,y,t)$ i platen gitt fornuftige rand- og initialkrav. 
Vi skal ikke løse varmeligningen i to romlige dimensjoner, men derimot betrakte spesialtilfellet der $t \to \infty$, slik at varmestrømmen er konstant, eller steady-state flow som vi sier på norsk. 
Dersom varmestrømmen er konstant, er $u_t(x,y,t)=0$, slik at varmeligningen blir
\[
u_{xx}+u_{yy}=0.
\]
Dette kalles Laplaces ligning, og den kan vi løse med teknikker vi har lært i kurset. Her er en figur som beskriver en typisk situasjon. 
\begin{figure}[htbp]
  \begin{center}
	\includegraphics[scale=.2]{../figurer/fig_heat.pdf}
	\label{fig:Num1}
	\end{center}
\end{figure}
Dette er temperaturen i en tynn plate der temperaturen holdes konstant lik null på tre sider, og konstant lik fem eller en eller noe på den siste siden. 

\subsubsection*{Løsning på begrenset intervall - separasjon av variable}

Vi setter opp problemet 
\[
u_{xx}+u_{yy}=0,
\]
på rektangelet $[0,a] \times [0,b]$ med randkrav 
\[
u(x,0)=u(0,y)=u(a,y)=0
\]
og
\[
u(x,b)=f(x).
\]
Vi separerer i vei
\[
u(x,y)=F(x)G(y).
\]
På samme vis som for bølgeligningen og varmeligningen, får vi 
\[
F''(x)+kF(x)=0
\]
og 
\[
G''(y)-kG(y)=0.
\]

Å finne $F$ og $k$ blir eksakt repetisjon av argumentet for bølgeligningen, mens 
\[
G_n(y)=A_n e^{-\frac{n\pi }{a}y} +B_n e^{\frac{n\pi }{a}y}
\]
slik at de generelle løsningene som tilfredsstiller $u(0,y)=u(b,y)=0$ blir
\begin{align*}
u_n(x,y)&=F(x)G_n(y)\\
&=\left(A_n e^{-\frac{n\pi }{a}y} +B_n e^{\frac{n\pi }{a}y}\right) \sin \frac{n\pi}{a}x.
\end{align*}
Krever vi $u(x,0)=0$, får vi
\[
A_n=-B_n,
\]
slik at 
\begin{align*}
u_n(x,y)&=A_n\left( e^{-\frac{n\pi }{a}y} -e^{\frac{n\pi }{a}y}\right) \sin \frac{n\pi}{a}x\\[3pt]
&=A_n\sinh \frac{n\pi }{a}y \sin \frac{n\pi}{a}x,
\end{align*}
og krever vi $u(x,b)=f(x)$, får vi det til ved å legge sammen alle løsninger
\[
u(x,y)=\sum_{n=1}^{\infty} u_n(x,y)=\sum_{n=1}^{\infty} A_n\sinh \frac{n\pi }{a}y \sin \frac{n\pi}{a}x,
\]
og kreve 
\[
f(x)=u(x,b)=\sum_{n=1}^{\infty} A_n\sinh \frac{n\pi b}{a} \sin \frac{n\pi}{a}x,
\]
slik at 
\[
A_n=\frac{2}{a\sinh \frac{n\pi b}{a}} \int_{0}^a f(x) \sin \frac{n\pi x}{a}\; dx.
\]


\subsubsection*{Løsning i halvplanet}

Vi skal se på Laplaces ligning
\[
u_{xx}+u_{yy}=0,
\]
i halvplanet. Randkrav skal være 
\[
\lim_{x \to \pm \infty} u(x,y)=\lim_{y \to \infty} u(x,y)=0
\]
og 
\[
u(x,0)=f(x),
\]
der 
\[
\lim_{x \to \pm \infty} f(x)=0.
\]

Siden $\lim_{x \to \pm \infty} u(x,y)=0$ satser vi på at $u$ er absolutt integrerbar i $x$, slik at vi kan fouriertransformere med hensyn på $x$. Da får vi
\[
\mathcal{F}(u_{xx})+\mathcal{F}(u_{yy})=\mathcal{F}(0)
\]
eller 
\[
-w^2\hat{u}(w,y)+\hat{u}_{yy}(w,y)=0.
\]

Akkurat som med varmeligningen, later vi som om dette er en ordinær differensialligning i $y$, og får 
\[
\hat{u} (w,y)= A(w) e^{-|w|y}+B(w) e^{|w[y}.
\]
Nå må vi fouriertransformere randkravet $\lim_{y \to \infty} u(x,y)=0$, og da får vi
\[
\lim_{y \to \infty} \hat{u}(w,y)=0.
\]
Bruker vi dette, ser vi at $B(w)=0$, og følgelig er 
\[
\hat{u} (w,y)= A(w) e^{-|w|y}.
\]
Transformerer vi $u(x,0)=f(x)$, får vi 
\[
\hat{u}(w,0)=\hat{f}(w),
\]
og bruker vi denne, får vi
\[
\hat{u} (w,y)= \hat{f}(w) e^{-|w|y}.
\]


Nå kan vi inverstransformere og få
\begin{align*}
u(x,y)&= \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\hat{f}(w) e^{-|w|y} e^{iwx}\; dw\\ &=\mathcal{F}^{-1}(\hat{f}(w) e^{-|w|y}).
\end{align*}
Nå bruker vi samme trikset som for varmeligningen. Siden dette er inversfouriertransformen til produktet mellom $\hat{f}(w)$ og $e^{-|w|y}$, og
\[
\mathcal{F}\left(\sqrt{\frac{2}{\pi}}\frac{y}{x^2+y^2}\right)=e^{-|w|y},
\]
ser vi at 
\begin{align*}
u(x,y)&=\frac{1}{\sqrt{2\pi}} f\ast \left(\sqrt{\frac{2}{\pi}}\frac{y}{x^2+y^2}\right)\\&=\frac{1}{\pi} \int_{-\infty}^{\infty} f(v) \frac{y}{(x-v)^2+y^2}\; dv.
\end{align*}

\subsubsection*{Noen artige egenskaper}
En funksjon som tilfredsstiller 
\[
u_{xx}+u_{yy}=0
\]
kalles gjerne en \defterm{harmonisk funksjon}.
Harmoniske funksjoner har morsomme egenskaper.


\begin{tcolorbox}[colback=red!10, title=Middelverdisats for harmoniske funksjoner]
\begin{thm}
La $u$ være harmonisk på $\Omega$, 
og la $\mathcal{C} \in \Omega$ være en sirkelskive med sentrum i $(x_0,y_0)$ og radius $a$. 
Da gjelder at
\[
u(x_0,y_0)=\frac{1}{2\pi a}\int_{\partial \mathcal{C}} u \; ds=\frac{1}{\pi a^2}\iint_{ \mathcal{C}} u \; dA.
\]
\end{thm}
\vspace{.3cm}
\end{tcolorbox}

Vi lar denne stå ubevist inntil videre. 
Beviset kommer ikke på eksamen i år.
Ta kontakt om du lurer.

\begin{tcolorbox}[colback=red!10, title=Maksverdisats for harmoniske funksjoner]
\begin{thm}
La $u$ være harmonisk på $\Omega$, 
og la $\mathcal{C} \in \Omega$ være en sirkelskive med sentrum i $(x_0,y_0)$ og radius $a$. 
Da ligger 
\[
\max_{(x,y) \in \mathcal{C}} u(x,y)
\]
og 
\[
\min_{(x,y) \in \mathcal{C}} u(x,y)
\]
på $\partial \mathcal{C}$.
\end{thm}
\vspace{.3cm}
\end{tcolorbox}
Vi beviser ikke denne heller.

\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
\[
\frac{f(x+h)-f(x)}{h}=f'(x)+\frac{f''(x)}{2}h+\cdots
\]
og vi ser at dette  stigninstallet består av den eksakte verdien
for $f'(x)$ pluss resten av taylorrekken til $f$
\[
\frac{f''(x)}{2}h+\frac{f'''(x)}{6}h^2+\cdots
\]
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
\[
\frac{f''(x)}{2}h+\frac{f'''(x)}{6}h^2+\cdots=O(h)
\]
for å signalisere at feilen er proporsjonal med $h$.



Vi kan relativt lett forbedre den lineære tilnærmingen fra forrige avsnitt ved å skrive
\[
f(x-h)=f(x)-f'(x)h+\frac{f''(x)}{2}h^2+\cdots
\]
 og sette opp tilnærmingen
\[
\frac{f(x+h)-f(x-h)}{2h}=f'(x)+\frac{f'''(x)}{3!}h^2+\frac{f^5(x)}{5!}h^4+\cdots
\] 
Feilen i denne tilnærmingen er 
\[
\frac{f'''(x)}{3!}h^2+\frac{f^5(x)}{5!}h^4+\cdots
\]
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
\[
\frac{f(x+h)-2f(x)+f(x-h)}{h^2}=f''(x)+O(h^2)%\frac{f''''(x)}{12}h^2+\cdots
\]
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
\[
u_t=u_{xx}
\]
med randkrav
\[
u(0,t)=u(1,t)=0
\]
og 
\[
u(x,0)=f(x).
\]
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$:
\[
u_{xx}(x,t) \approx \frac{u(x+h,t)-2u(x,t)+u(x-h,t)}{h^2}
\]
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:
\[
u_{xx}(x_i,t) \approx \frac{u_{i+1}(t)-2u_{ij}(t)+u_{i-1}(t)}{h^2},
\]
Nå setter vi inn dette uttrykket i varmelikningen: 
\[
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},
\]
Siden $u_i$ kun er en funksjon av $t$, passer det å skrive 
\[
u_i'(t) = \frac{u_{i+1}(t)-2u_{ij}(t)+u_{i-1}(t)}{h^2},
\]
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
  \[
\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}.
\]
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
\[
 (I+\frac{k}{h^2} A) \mathbf{u}_{j+1}=\mathbf{u}_{j}
\]
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 
%\[
%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}
%\]
%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 
\[
 (2I+\frac{k}{h^2} A)\mathbf{u}_{j+1}= (2I-\frac{k}{h^2} A)\mathbf{u}_{j}
\]
å 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
\[
u_{xx}+u_{yy}=0,
\]
på kvadratet $[0,1] \times [0,1]$ med randkrav 
\[
u(x,0)=u(0,y)=u(1,y)=0 \qquad \text{og} \qquad u(x,1)=f(x).
\]
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
\[
u_{xx}(x,t) \approx \frac{u(x+h,t)-2u(x,t)+u(x-h,t)}{h^2}
\]
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]{../figurer/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  
 \[
u_t=u_{xx}+u_{yy}
\]
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
\[
u(x,y,0)=g(x,y).
\]

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.
\]



% \[
%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}
%\]
% og 
% \[
%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}
%\]



\kapittelslutt
