‹

TMA4212: Numerisk løsning av differensialligninger med differansemetoder

Tags:

$\newcommand{\C}{\mathbb{C}}$ $\newcommand{\R}{\mathbb{R}}$ $\newcommand{\O}{\mathcal{O}}$

Introduksjon

Gershgorin's teorem

La $A = (a_{ik}) \in \C^{n\times n}$ være en kvadratisk matrise. Definer $$S = \bigcup_{j} \left\{ z \in \C: \lvert z - a_{jj} \rvert \leq \sum_{k \neq j} \lvert a_{jk} \rvert \right\},$$ til å være unionen av sirklene sentrert rundt de diagonale elementene i $A$, med radius lik summen av absoluttverdien til de resterende elementene, respektivt på hver rad. Da inneholder $S$ alle egenverdiene til $A$.

Spektralradius

La $A = (a_{ik}) \in \C^{n\times n}$ være en kvadratisk matrise. Spektralradiusen til $A$ er definert som egenverdien til $A$ med størst absoluttverdi: $$\rho(A) = \max(\{\lambda_1, \lambda_2,\ldots,\lambda_n\}).$$

Vektornormer

Vektorrom kan utstyres med en norm. Uformelt kan vi si at en slik norm gir oss en metode for å måle elementer av vektorrommet. Naturlig nok, er den velkjente lengdeformelen $\mathrm{\lVert x \rVert} = \sqrt{x_1^2 + x_2^2 + \ldots + x_n^2}$ en norm over $\R^n$.

La $X$ være et vektorrom. En vektornorm er en avbildning $\lVert \cdot \rVert: X \to \C\>(\R)$, slik at følgende krav er tilfredstilt:

$\lVert x \rVert \geq 0,\>\forall x \in X,\> \lVert x \rVert = 0 \iff x = 0$

$\lVert ax \rVert = \lvert a \rvert \lVert x \rVert,\>\forall a \in \C\>(\R), \forall x \in X$

$\lVert x + y \rVert \leq \lVert x \rVert + \lVert y \rVert$

Eksempler på vektorrom, for $x = [x_1, x_2,\ldots,x_n], X = \C^n$.

$$ \begin{align} \lVert x \rVert_1 &= \sum_i \lvert x_i \rvert \\ \lVert x \rVert_2 &= \sqrt{\sum_i \lvert x_i\rvert ^2} \\ \lVert x \rVert_p &= \left(\sum_i \lvert x_i\rvert^p\right)^{\frac 1p} \\ \lVert x \rVert_{max} &= \max_{1 \leq i \leq n} \lvert x_i \rvert \end{align} $$

Matrisenormer

Matriserom kan utstyres med normer på samme måter som vektorrom.

La $K^{m \times n}$ være et matriserom. En matrisenorm er en avbildning $\lVert \cdot \rVert: K \to \C\>(\R)$, slik at følgende krav er tilfredstilt:

De tre kravene for vektornormer.

$\lVert A\cdot B \rVert \leq \lVert A \rVert \cdot \lVert B \rVert$, hvis matrisemultiplikasjon er definert for $A,B \in K^{m \times n}$.

Eksempler på matrisenormer, for $A = (a_{ik}), K = C^{m \times n}$:

$$ \begin{align} \lVert A \rVert_1 &= \max_{1 \leq k \leq n} \sum_{i=1}^n \lvert a_ik \rvert \\ \lVert A \rVert_2 &= \sqrt{\rho(A^TA)} \\ \lVert A \rVert_{\infty} &= \max_{1 \leq i \leq n} \sum_{k=1}^n \lvert a_ik \rvert \\ \lVert A \rVert_{\mathrm{frobenius}} &= \left(\sum_{i=1}^n \sum_{j=1}^n \lvert a_ij\rvert^2\right)^{\frac 12} \end{align} $$

Stor-$\O$ notasjon

Stor-$\O$ notasjon brukes til å beskrive en øvre grense for hvordan en funksjon oppfører seg nær en grense. Ta for eksempel $f(x) = x^2$. Vi sier at $f(x) = \O(x^3),\> x \to \infty$, siden $x^3$ vokser fortere enn $x^2$ når $x$ øker mot uendelig. Merk at dette ikke betyr at $x^2 < x^3\>\forall x$. Siden vi her snakker om øvre grenser, er det klart at $f(x)$ har mange stor-$\O$ ekvivalenser: $f(x) = \O(x^4), f(x) = \O(2^x), f(x) = \O(x^2)$.

Tilsvarende kan vi se på stor-$\O$ notasjon ved grensen $0$. Ta igjen $f(x) = \O(x^2)$. Her er $f(x) = \O(x), f(x) = \O(x^2)$, men ikke $f(x) = \O(x^3)$, da $x^2$ er (mye) større enn $x^3$ når $x$ nærmer seg $0$ (ovenifra). Denne måten å beskrive hvor fort en funksjon "nærmer seg null" ved hjelp av polynomene $x^n,\>n \in \mathbb{Z}^+$ er hyppig brukt i numerisk matematikk.

La $f, g$ være to funksjoner $f, g:C \to C$. Da er $f(x) = \O(g(x))$ i grensen mot $x_0 \in R \cup \infty$ hvis det finnes $M, \delta > 0$ slik at

$$ \lvert f(x) \rvert \leq M \lvert g(x) \rvert, \forall x; \lvert x - x_o \rvert < \delta$$

Definisjonen sier med andre ord at vi kan alltids finne en konstant $M$ slik at funksjonen $f$ vil være øvrig bundet av $g$ ganget med $M$, tilstrekkelig nærme grensen vår.

La $f$ være en funksjon $f:C\to C$. La $f(x) = \O(x^p)$ mot en grensen $x_0 \in R \cup \infty$(typisk 0). Vi sier da at $f$ er av orden p.

Taylorrekker

La $f(x)$ være en kontinuerlig $n+1$-deriverbar funksjon på et intervall $I$. Da har vi for $x+h \in I$:

$$ f(x+h) = f(x) + h f'(x) + h^2 \frac{f''(x)}{2!} + h^3 \frac{f'''(x)}{3!} + \dots + h^n \frac{f^{(n)}(x)}{n!} + R_n $$

hvor $R_n = h^{n+1}\frac{f^{(n+1)}(x + \theta h)}{(n+1)!},\>\theta \in [0, 1]$. Eventuelt kan vi si at $R_n = \O(h^{n+1})$.

Diskretisering og differanseoperatorer

Grunnleggende metode

Vi ønsker å løse (partielle) differensiallikninger. Det lar seg ikke alltid gjøre så enkelt analytisk. Istedet for den analytiske tilnærmingen, ønsker vi dermed å anvende numeriske metoder for å estimere løsningen.

Vi begynner med definisjonen av den deriverte:

$$ f'(x) = \lim_{h\to0} \frac{f(x+h) - f(x)}{h} $$

Det er tydelig at vi kan estimere $f'(x)$ ved å sette inn en verdi for $h$. Estimatet vil forbedres ved å gjøre $h$ så liten som nødvendig. Så hvor stor er egentlig feilen her? Ser vi på taylorrekken rundt $f(x+h)$ får vi

$$ f(x+h) = f(x) + h f'(x) + \O(h^2) \implies f'(x) = \frac{f(x+h) - f(x)}{h} + \O(h) $$

Vi har med andre ord en tilnærming i første orden av h til den deriverte. Vi har nå en formel vi kan bruke til å iterativt regne ut differensiallikninger av første orden. Ta eksempelvis:

$$ u'(x) - 3 = 5u(x),\>u(0) = 10,\> x \in [0, 3] $$

Vi diskretiserer domenet, ved å la $h=0.1$. Vi kan nå regne ut $u(x)$ ved å "hoppe" steg med lengde $h$:

$$ \begin{align} u'(x) - 3 =& 5u(x) \\ \frac{u(x+h) - u(x)}{h} - 3 =& 5u(x) \\ u(x+h) - u(x) =& h(5u(x) + 3) \\ u(x+h) =& u(x) + h(5u(x) + 3) \\ \end{align} $$

Eksplisitt for de 3 første stegene:

$$ \begin{align} u(0.1) =& u(0) + 0.1\cdot(5\cdot u(0) + 3) = 15.3 \\ u(0.2) =& u(0.1) + 0.1\cdot(5\cdot u(0.1) + 3) = 23.25 \\ u(0.3) =& u(0.2) + 0.1\cdot(5\cdot u(0.2) + 3) = 35.175 \\ \end{align} $$

De tilsvarende analytiske verdiene for systemet er $16.88$, $28.21$, og $46.91$.

Differanseoperatorer

Det er flere måter å approksimere den deriverte på. Og da vi etterhvert har lyst til å arbeide med deriverte av høyere grad, trenger vi å introdusere flere differanseoperatorer:

Deriverte av første grad

Foroverdifferens	$\Delta u(x) = u(x+h) - u(x)$
Bakoverdifferens	$\nabla u(x) = u(x) - u(x-h)$
Sentraldifferens	$\delta u(x) = u(x+h/2) - u(x - h/2)$
Middelverdi	$\mu u(x) = \frac{1}{2}(u(x + h/2) + u(x - h/2))$

Disse kan da kombineres på interessante måter til å approksimere deriverte. Eksempelvis er det vanlig å approksimere andrederiverte ved å anvende sentraldifferense to ganger:

Deriverte av andre grad

$$ \delta^2 u(x) = \delta(u(x+h/2) - u(x - h/2)) = u(x+h) - u(x) - (u(x) - u(x-h)) = u(x+h) - 2u(x) - u(x-h)$$

Sentraldifferense i annen	$u(x+h) - 2u(x) + u(x-h)$
Middelverdi i annen	$\frac 14 (u(x+h) + 2u(x) + u(x-h)$

Som regel (i dette kurset) er det sentraldifferansen i annen som blir anvendt til approksimering av annen-deriverte.

Notasjon

Vi bruker fremover følgende notasjon for deriverte og differanser:

$$ \begin{align} \frac{\partial u}{\partial x} &= u_x \\ \frac{\partial u}{\partial t} &= u_t \\ u(mh) &= u_m \\ u(mh+h) &= u_{m+1} \\ u(mh, nk) &= u_m^n \\ u(mh+h, nk+k) &= u_{m+1}^{n+1} \\ \end{align} $$

For en partiell differensiallikning $Pu = f$, hvor $P$ er en lineær kombinasjon av n'te-deriverte; denoterer vi en bestemt differanseapproksimasjon til $u$ ved $P_k u$(én variabel) eller $P_{kh} u$(to variable).

Eksplisitte og implisitte metoder

Eksplisitte metoder

Eksplisitte metoder er metoder som er avhengig bare av de forige tidsstegene for å regne ut neste tidssteg. Dette betyr at utregningen av det neste tidssteget kan gjøres eksplisitt med en lineær kombinasjon av de tidligere tidsstegene. Ta for eksempel liknigen $u_t - u_x = 0$. Vi diskretiserer ved

$$ \begin{align} \frac{u_m^{n+1} - u_m^n}{k} - \frac{u_m^n - u_{m-1}^n}{h} = 0 \Longleftrightarrow \\ u_m^{n+1} = (1 + \lambda)u_m^n + \lambda u_{m-1}^n,\quad \lambda = k/h \end{align} $$

Vi kan her altså regne ut neste tidssteg direkte ved å bruke resultater fra de tidligere tidsstegene.

Implisitte metoder

Vi kan derimot diskretisere "rundt" et tidssteg $(n+1)k$, som gir oss en implisitt metode:

$$ \begin{align} \frac{u_m^{n+1} - u_m^n}{k} - \frac{u_m^{n+1} - u_{m-1}^{n+1}}{h} = 0 \Longleftrightarrow \\ (1 - \lambda)u_m^{n+1} - \lambda u_{m-1}^{n+1} = \frac{1}{h}u_m^n \\ \end{align} $$

Her har vi et nytt problem, i at vi på venstresiden har en differanse i tidssteg $n+1$, som vi altså ønsker å finne. Dette betyr at vi på hvert tidsteg må løse $m$ likninger. Implisitte metoder har ofte (i konteksten av problemene vi skal løse) "greiere" stabilitetskrav enn eksplisitte metoder.

Avbruddsfeil (Truncation error)

Avbruddsfeil defineres som feilen introdusert ved å "brekke av" en uendelig sum ved et endelig steg. Ta for eksempel taylorrekken ovenfor

$$ f(x+h) = f(x) + h f'(x) + h^2 \frac{f''(x)}{2!} + h^3 \frac{f'''(x)}{3!} + \dots + h^n \frac{f^{(n)}(x)}{n!} + R_n. $$

Om vi her approksimerer $f(x+h)$ved å bruke de $n$ første leddene, får vi en avbruddsfeil $\tau = R_n = \O(x^{n+1})$.

Vi får som regel en avbruddsfeil ved å bruke differansemetoder på en differensiallikning, som er avhengig av hvilke differenaseoperatorer vi bruker til å estimere de deriverte. Ta for eksempel, som før:

$$ u_x - 3 - 5u(x) = 0 $$

Vi estimerer $u_x$ ved foroverdifferens som tidligere, og får:

$$ \begin{align} \frac{1}{h}(u_{n+1} - u_{n}) - 3 - 5u_n &= 0 \\ \end{align} $$

Vi ekspanderer $u_{n+1}$ i en taylorrekke, og får

$$ \begin{align} u_{n+1} &= u_n + h u_x + \frac 12 h^2 u_{xx} + \O(h^3) \implies \\ P_{kh} &= \frac{1}{h}(u_{n+1} - u_{n}) - 3 - 5u_n = 0 \implies \\ P_{kh} &= u_x + \frac{h}{2} u_{xx} + \O(h^2) - 3 - 5u(x) = 0 \end{align} $$

merk her at $u_n = u(x)$.

Så vi har nå funnet et utrykk som likner på den originale differensiallikningen, med en avbruddsfeil, nemlig $\frac{h}{2}u_{xx} + \O(h^2)$.

Bruker vi notasjonen $P_{kh}$ for differanseapproksimasjoner, kan vi definere avbruddsfeilen $$\tau_n^m = Pu - P_{kh}$$

Konsistens

En differanseapproksimasjon er konsistent hvis

$\tau_n^m \to 0$ for alle $n, m$ når $k, h \to 0$.

Konsistens er ett av kravene for

Stabilitet

Uformelt beskriver stabilitet hvorvidt en differansemetode vil holde seg innenfor et begrenset område, i løpet av estimeringsprosessen for et gitt system. Kontrapositivt vil en ustabil metode kunne "eksplodere", oppføre seg sporadisk og aldri tilnærme seg en løsning. Eventuelt kan vi si at stabilitet betyr at en differansemetode ikke vil øke feil fra steg til steg i approksimasjonen.

La $\lVert \cdot \rVert_h$ være en skalert variant av den vanlige $\lVert \cdot \rVert_2$ normen for vektorrom. Da har vi følgende:

En differanseapproksimasjon $P_{kh}u$ er stabil i en region $\Lambda$ hvis det finnes et heltall $J$ slik at for enhver positiv tid $T$, så er det en konstant $C_T$ slik at

$$\lVert u^n \rVert_k \leq C_T \sum_{j=1}^J \lVert u^j \rVert_k$$ for enhver $0 \leq nk \leq T$, $(k, h) \in \Lambda$.

Definisjonen forteller oss at vi ved et punkt $J$, burde bare være en konstant unna å endre på approksimasjonen vår så mye som vi kommer til å gjøre i løpet av ethvert tidsløp $T$.

Definisjonen er dessverre vanskelig å anvende, og vi bruker som regel andre verktøy for å vise stabilitet.

En-stegs metoder

La oss se på en generell en-stegs differanseapproksimasjon $AU^{n+1} = BU^n + c^n$, som kan skrives $$U^{n+1} = CU^n + q^n,\>C=A^{-1}B,\>q^n = A^{-1}c^n.$$ Vi har for dette systemet en enklere definisjon for stabilitet;

La $U^{n+1} = CU^n + q^n$ være en differanseapproksimasjon for en PDE i $[0,1] \times [0, T]$. Denne er stabil hvis det finnes en konstant $L$ uavhengig av $l$ og $k$ slik at for enhver vektor $w^0$ og sekvensen $w^{n+1} = Cw^n$ så har vi

$$ \lVert w^n \rVert \leq L \lVert w^0 \rVert,\> 0 \leq nk \leq T$$

Om vi lar $T = \infty$ så kalles stabiliteten for $F$-stabilitet. Ekvivalent, så har vi

La $U^{n+1} = CU^n + q^n$ være en differanseapproksimasjon for en PDE i $[0,1] \times [0, T]$. Denne er stabil hvis det finnes en $L$ slik at

$$\lVert C^n \rVert \leq L,0 \leq nk \leq T$$

Vi kan formilde kravet for stabilitet noe:

La $U^{n+1} = CU^n + q^n$ være en differanseapproksimasjon for en PDE i $[0,1] \times [0, T]$. Denne er stabil hvis det finnes en $\mu$ uavhengig av $k$ og $h$ slik at

$$\lVert C \rVert \leq 1 + \mu k$$

Til slutt er følgende et krav for stabilitet (men er ikke alltid tilstrekkelig for å garantere stabilitet):

La $U^{n+1} = CU^n + q^n$ være en differanseapproksimasjon for en PDE i $[0,1] \times [0, T]$. Om den er stabil, finnes det en $\nu \geq 0$ uavhengig av $k$ og $h$ slik at

$$\rho(C) \leq 1 + \nu k$$

Da mange av metodene vi jobber med er en-stegs metoder, er definisjonene og kravene ovenfor nyttige.

Von-Neumann stabilitet

Metoden som ofte brukes til å (av)vise stabilitet, er Von-Neumann analyse. Metoden kan kort beskrives ved følgende steg:

Erstatt alle $u_n^m$ i differanseapproksimasjonen med $g^n e^{im\theta}$.
Manipuler til et utrykk for $g$.
Bruk en av de følgende teoremene for å (av)vise stabilitet.

Funksjonen $g(\theta)$ kalles for amplifikasjonsfaktoren til differanseapproksimasjonen.

En ett-stegs differanseapproksimasjon er stabil i en region $\Lambda$ hvis og bare hvis det finnes en konstant $\mu$ (uavhengig av $\theta$, $k$ og $h$, slik at

$$\lvert g(\theta, k, h) \rvert \leq 1 + \mu k$$

Om $g$ også er uavhengig av $k$ og $h$, får vi følgende:

En ett-stegs differanseapproksimasjon er stabil i en region $\Lambda$ hvis og bare hvis det finnes en konstant $\mu$ (uavhengig av $\theta$, $k$ og $h$, slik at

$$\lvert g(\theta) \rvert \leq 1$$

Konvergens

Det store spørsmålet vi er ute etter når vi anvender en numerisk metode for å løse en gitt PDE, er om metoden vår konvergerer. Med andre ord er vi interessert i å vite om den approksimerte løsningen vår tilnærmer seg den analytiske løsnigen.

Etter en runde med en numeriske metode, har vi en approksimert løsning $U^n$. Den analytiske løsnigen er $u^n$. Vi definerer feilen til løsningen $U^n$ ved $e^n = U^n - u$. Vi sier så at $U \to u$ i normen $\lVert \cdot \rVert$, hvis

$$ \max_{0 \leq nk \leq T} \lVert e^n \rVert \to 0 $$

når $k\to 0, h\to 0$.

Lax' ekvivalensteorem

En differanseapproksimasjon er konvergent hvis og bare hvis den er både stabil og konsistent.

Wikipendium