\input{preamble}

\title{ 48 - NUMERISK LØSNING}

\begin{document}
\pagenumbering{gobble}

\maketitle








Når man skal løse sin første partielle differensiallikning numerisk,
er varmelikningen grei å begynne med, 
for dette er den enkleste å få til.



Den deriverte til en funksjon $f$ er
\[
f'(x)=\lim_{h \to 0} \frac{f(x+h)-f(x)}{h}
\]
dersom grensen eksisterer.
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. 


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}$.

\begin{oppgave}{1}
Gjenta eksperimentet med  $h=0.01$, $h=0.001$ osv,
og se hva som skjer. 
Hvor liten $h$ kan du ta før det går åt skogen?
\end{oppgave}


Feilen i forrige  er tydelig proporsjonal med $h$ - 
deler du $h$ på 10, deler du feilen på 10. 
En god illustrasjon av lineær feil. 
Men dette gjelder bare inntil et visst punkt, 
blir $h$ liten nok blir presisjonen dårligere igjen. 
Dette er fordi datamaskinen din kun regner med 16 desimaler. 
Skal vi få til noe bedre, 
må vi bruke en annen formel.


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$.


\begin{oppgave}{2}
Gjenta eksperimentet i forrige oppgave, 
men nå med formelen
\[
\frac{f(x+h)-f(x-h)}{2h}=f'(x)
\] 
og se hva som skjer. 
Hvor liten $h$ kan du nå ta før det går åt skogen?
Bruk taylorrekker til å forklare oppførselen.
\end{oppgave}


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.


\begin{oppgave}{3}
Hvis du virkelig vil slå på stortrommen, kan du bruke formelen
\[
f'(x)= \frac{f(x-2h)-8f(x-h)+8f(x+h)-f(x+2h)}{12h}
\]
Prøv.
\end{oppgave}


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)+\frac{f''''(x)}{12}h^2+\cdots
\]
Det finnes bedre og mer kompliserte formler, 
men vi skal ikke plages med dem.
Vi skal nå lage en numerisk metode 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*}
Det er relativt lett å konstruere liknende metoder for andre likninger
men varmelikningen er som sagt den enkleste å begynne med.


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.

\begin{tcolorbox}
Vi skal finne approksimasjoner til løsningen $u$ i alle gitterpunktene.
\end{tcolorbox}

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}


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*}
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_{i}(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_{i}(t)+u_{i-1}(t)}{h^2},
\end{equation*}
Siden $u_i$ kun er en funksjon av $t$, 
passer det å skrive alt dette som et koblet system av differensiallikninger: 
\[
u_i'(t) = \frac{u_{i+1}(t)-2u_{ij}(t)+u_{i-1}(t)}{h^2}
\]
Systemet har $n-1$ likninger, 
for $u_0=0$ og $u_n=0$ kjenner vi; 
disse er gitt av randbetingelsene.


Nå kan vi løse systemet med en ønsket numerisk metode. 
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:
\[
\frac{u_{i,j+1}-u_{ij}}{k} = \frac{u_{i+1,j}-2u_{i,j}  + u_{i-1,j}}{2h^2} +\frac{u_{i+1,j+1}-2u_{i,j+1}+u_{i-1,j+1}}{2h^2}
\]

Metodene har litt forskjellige egenskaper. 
Vi får se hva vi rekker av analyse.

\begin{oppgave}{4}
Løs varmelikningen med eksplisitt skjema. 
Prøv forskjellige kombinasjoner av $h$ og $k$.\\
(Metoden funker ikke med mindre $k/h^2\geq 1/2$.)
\end{oppgave}

\begin{oppgave}{5}
Kravet $k/h^2\geq 1/2$ er ganske restriktivt,
og det er derfor vi har implisitte metoder. 
De er litt vanskeligere å implementere, 
og litt mer beregningstunge,
men slipper unna slike hårreisende betingelser på $h$ og $k$.
Implementer implisitt.
\end{oppgave}


\begin{oppgave}{6}
Men hvorfor har vi Crank-Nicolson?
Implementer Crank-Nicolson, 
og sammenlikne de tre metodene for samme $h$ og $k$ med den 
analytiske løsningen til varmelikningen med for eksempel initialkrav $f(x)=\sin x$.
\end{oppgave}
\begin{center}
\includegraphics[scale=.26]{../figurer/IMG_0319}
\end{center}



Hvis man har orket å komme seg gjennom alt dette greiene her, 
kommer belønningen i form av evnen til å modellere
problemer som er umulig å regne særlig på med penn og papir. 


\begin{oppgave}{7}
Gå helt tilbake til start, 
og utled at varmelikningen for anisotrope media blir: 
\[
u_t=\nabla\cdot \left(\kappa\nabla u\right) 
\]
Denne går ikke å løse analytisk. 
Modifiser kodene dine til å takle dette i en romlig dimensjon.
\end{oppgave}

\begin{oppgave}{8}
Utvid til løse varmelikningen i to romlige koordinater:
\[
\rho u_t=\Delta u
\]
Eksplisitt skjema er ikke så vanskelig, 
men for implisitt og C-N er det litt knot å sette opp korrekte matriser.
\end{oppgave}


\begin{oppgave}{9}
Løs varmelikningen 
og løs varmelikningen numerisk på rektangulær boks $\Omega$ med 
initialkrav 
\[
u(\mathbf x,0)=f(\mathbf x)
\] 
og randkrav 
\[
u(\mathbf x,t)=0
\] 
på $\partial \Omega$.
\end{oppgave}
\includegraphics[scale=.3]{../figurer/M1004345}

\begin{oppgave}{10}
Du kan også løse bølgelikningen slik, 
men da må du knote litt mer med den doble tidsderiverte.
\end{oppgave}

Og du kan diskretisere laplaces likning slik:
\[
0 =\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^2}  + \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h^2}  
\]
Man får en likning for hvert indre gitterpunkt,
og hvis man organiserer de ukjente i en kolonnevektor, 
får man et likningssystem på formen 
\[
A\mathbf u = f
\]
der $f$ er en vektor med randverdier.
 

\begin{center}
\hspace*{-3em}
\begin{tikzpicture}[thick]
\draw[->] (-1.2,-1) -- (4,-1) node[right]{$x$};
\draw[-] (-1.2,2) -- (4*.6,2) ;
\draw[->] (-2*.6,-1) -- (-2*.6,3) node[right]{$y$};
\draw[-] (4*.6,-1) -- (4*.6,2) ;


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

\draw[line width = 1mm] (2*.6,-.95) -- (2*.6,-1.05) node[below]{$x_i=ih$};
\draw[line width = 1mm] (-2*.586,) -- (-2*.614,1) node[left]{$y_j=jh$};

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

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

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

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


\end{tikzpicture}
\end{center}



\begin{oppgave}{11}
Finn matrisen $A$, 
\end{oppgave}

\begin{oppgave}{12}
Løs problemet
\[
u_{xx}+u_{yy}=0
\]
på rektangelet $[0,1] \times [0,1]$ med randkrav 
\[
u(x,0)=u(0,y)=u(1,y)=0
\]
og
\[
u(x,1)=f(x).
\]
\end{oppgave}



\end{document}
