Gå till innehållet

SF1547 Numeriska metoder, grundkurs 6,0 hp

Tillfälleskod Termin(er) Period(er) Föreläsare
60945 VT2025 3 Katarina Gustavsson

Sammanfattning av föreläsningarna VT2025 (period 3).

Föreläsning 1

Minsta kvadratmetoden

Mål: Hitta de modellparametrar som minimerar summan av kvadratavvikelserna (modell kontra data), dvs. som minimerar \(||\vec r||_2^2=||A\vec c-\vec y||_2^2\).

Exempel: Anpassa modellen \(f(t)=c_0+c_1t+c_2t^2+c_3\sin(kt)+c_4\cos(kt)\) efter data (antag \(k=2\pi\)).

Lösning

Vi vill lösa ekvationssystemet

\[\begin{bmatrix}1&t_1&t_1^2&\sin(kt_1)&\cos(kt_1)\\1&t_2&t_2^2&\sin(kt_2)&\cos(kt_2)\\\vdots&\vdots&\vdots&\vdots&\vdots\\1&t_n&t_n^2&\sin(kt_n)&\cos(kt_n)\end{bmatrix}\begin{bmatrix}c_0\\c_1\\c_2\\c_3\\c_4\end{bmatrix}=\begin{bmatrix}y_1\\y_2\\\vdots\\y_n\end{bmatrix}\]

Systemet är överbestämt. Minstakvadratlösningen fås genom att istället lösa

\[A^TA\vec c=A^T\vec y\]

Linjära ersättningsmodeller

Anta att vi vill anpassa data till en icke-linjär modell \(p=ae^{bt}\). Genom att logaritmera båda led får vi en linjär ersättningsmodell:

\[y=m+bt\]

där \(y=\ln p\) och \(m=\ln a\).

Minstakvadratlösningen ger oss värden på \(m\) och \(b\). Dessa kan sedan användas för att direkt beräkna \(a=e^m\).

Föreläsning 2

Polynominterpolation

Mål: Hitta ett polynom som går genom givna datapunkter.

Antag att vi har \(n\) datapunkter. Vi vill hitta ett polynom (av gradtal \(n-1\); kan möjligen visa sig bli lägre) som passerar genom dessa. Ett sådant polynom är alltid unikt.

Naiv ansats

\[p_{n-1}(x)=c_0+c_1x+c_2x^2+\cdots+c_{n-1}x^{n-1}\]

För att hitta modellparametrarna \(c_0,\dots,c_{n-1}\) löses systemet

\[A\vec c=y\iff\begin{bmatrix}1&x_1&x_1^2&\cdots&x_1^{n-1}\\1&x_2&x_2^2&\cdots&x_2^{n-1}\\\vdots&\vdots&\vdots&\ddots&\vdots\\1&x_n&x_n^2&\cdots&x_n^{n-1}\end{bmatrix}\begin{bmatrix}c_0\\c_1\\c_2\\\vdots\\c_{n-1}\end{bmatrix}=\begin{bmatrix}y_1\\y_2\\\vdots\\y_n\end{bmatrix}\]

Matrisen \(A\) kallas för Vandermondematrisen och kan vara illa-konditionerad.

Newtons ansats

Antag att vi har \(n\) stycket datapunkter, alla med olika \(x\)-värden.

Vi ställer nu upp ett polynom på formen

\[p_{n-1}(x)=a_0+a_1(x-x_1)+a_2(x-x_1)(x-x_2)+\cdots+a_{n-1}(x-x_1)(x-x_2)\cdots(x-x_{n-1})\]

Obs! Båda dessa ansatser ger samma polynom, även om koefficienterna normalt inte är desamma.

Interpolationsfel

Maxfelet definieras som \(E_n=\max_{a\leq x\leq b}{|f(x)-p_n(x)|}\).

Går \(E_n\to0\) när \(n\to\infty\)?

Konstant funktion

En konstant funktion (\(n=0\)) för vilken \(p(x)=f(x_0)\) har felet

\[E_0=|f(x)-p(x)|=|f(x)-f(x_0)|=|f'(\xi)|(x-x_0)\]

för något \(\xi\in(x_0,x)\).

Linjär funktion (linjär interpolation)

Det interpolerade polynomet ger rätt värde i två punkter \(x_0\) och \(x_1\). Felet ges av

\[E_1=|f(x)-p(x)|=\left|\frac{f''(\xi)}2\right|(x-x_0)(x-x_1)\]

för något \(\xi\in(x_0,x_1)\).

Allmänt fel

Om \(p_n(x)\) är ett interpolerande polynom av gradtal \(n\) till funktionen \(f(x)\) på intervallet \(x=[a,b]\), och avståndet mellan punkterna är ekvidistant lika med \(h\), gäller att

\[E_n=\max_{a\leq x\leq b}\left|f(x)-p_n(x)\right|\leq\max_{a\leq\xi\leq b}\frac{|f^{n+1}(\xi)|}{4(n+1)}h^{n+1}\]

där \(f^{n+1}\) är originalfunktionens derivata av ordning \(n+1\).

Mer om interpolationsfel

Tumregel: \(E_1=O(h^2),E_2=O(h^3)\)

Runges fenomen:

Oscillationer uppträder mot intervallets ändar när \(n\to\infty\).

Lösning: Använd styckvis interpolation (dela upp intervallet och interpolera med polynom av lägre grad).

Föreläsning 3

Lösningar till olinjära ekvationer

Olinjära ekvationer \(f(x)=0\) kan inte lösas allmänt på analytisk väg. Istället kan en numerisk metod användas.

Följande metoder är alla iterativa och består av

  • en startgissning
  • ett stoppvillkor
  • en iterationsformel

Hur snabbt metoden hittar lösningen ges av den förstnämndes konvergenshastighet.

Intervallhalvering

Enkel metod. Beskrivs i boken.

Fixpunktsmetoden

Vi skriver om ekvationen på fixpunktsform (\(g(x)=x\)). Givet en startgissning \(x_0\) ges iterarationsformeln som:

\[x_{n+1}=g(x_n),n=0,1,2,\dots\]

Stoppvillkor: \(|x_{n+1}-x_n|\lt\tau\) där \(\tau\) är önskad tolerans.

För att fixpunktsiterationen ska konvergera mot fixpunkten måste det gälla att \(|g'(x)|\lt1\) i en omgivning till fixpunkten.

Fixpunktsmetoden konvergerar linjärt med hastigheten \(S=\lim_{n\to\infty}\frac{e_{n+1}}{e_n}=|g'(x^*)|\), där \(x^*\) är fixpunkten.

Härledning: av konvergensvillkoret ovan

(se bild #30 i föreläsning 3)

Newtons metod

Anta att vi har en ekvation \(f(x)=0\).

Iterationsformeln ges av

\[x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)},n=0,1,2,\dots\]

Nollstället till \(x_n\)-tangenten väljs som ny gissning \(x_{n+1}\).

Stoppvillkor: \(|x_{n+1}-x_n|\lt\tau\)

Newtons metod konvergerar kvadratiskt och lokalt om \(f\) är två gånger kontinuerligt deriverbar samt \(f'(x^*)\neq0\).

Härledning av konvergens... (se bild #14 i F4)

Obs! Satsen innebär enbart att Newtons metod är lokalt konvergent. Metoden kan dock misslyckas om \(f'(x_n)=0\) för någon iterand \(x_n\). I vissa fall kan metoden "fastna" och alternera mellan två eller flera värden.

Om \(f(x)\) har ett nollställe med högre multiplicitet än 1 (t.ex. dubbelrot) konvergerar Newtons metod endast linjärt.

Föreläsning 4

Lösningar till olinjära ekvationer (fortsättning)

Sekantmetoden

Liknar Newtons metod, men använder en approximation av derivatan.

Iterationsformeln är

\[x_{n+1}=x_n-f(x_n)\cdot\frac{x_n-x_{n-1}}{f(x_n)-f(x_{n-1})},n=1,2,3,\dots\]

Sekantmetoden kräver två startvärden: \(x_0\) och \(x_1\)

Sekantmetoden ger superlinjär konvergens (mellan linjär och kvadratisk).

Sammanfattning konvergenshastighet

Generellt kan konvergens hos en iterativ numerisk metod beskrivas av sambandet

\[\lim_{n\to\infty}\frac{e_{n+1}}{e_n^p}=S\]

där \(S\) kallas för den asymptotiska felkonstanten (fixpunktsmetoden: konvergenshastighet).

Givet att \(S\) är ett ändligt tal kallas \(p\) för konvergensordning.

  • Fixpunktsmetoden: \(p=1,S\lt1\) (linjär konvergens)
  • Sekantmetoden: \(1\lt p\lt2\) (superlinjär konvergens)
  • Newtons metod: \(p=2\) (kvadratisk konvergens)

Newtons metod för system

Vi skriver problemet på formen

\[\vec f(\vec x)=0\iff\begin{bmatrix}f_1(x_1,x_2,\dots,x_k)\\f_2(x_1,x_2,\dots,x_k)\\\vdots\\f_k(x_1,x_2,\dots,x_k)\end{bmatrix}=\vec0\]

och söker alltså vektorn \(\vec x\) s.a. \(\vec f(\vec x)=0\).

Iterationsformeln blir

\[\vec x_{n+1}=\vec x_n-J^{-1}(\vec x_n)\cdot\vec f(\vec x_n),n=0,1,2,\dots\]

vilket ofta skrivs om på formeln

\[J(\vec x_n)\left(\vec x_{n+1}-\vec x_n\right)=-\vec f(\vec x_n)\]

där jakobianen \(J(\vec x_n)=\begin{bmatrix}\frac{\partial f_1(\vec x)}{\partial x_1}&\frac{\partial f_1(\vec x)}{\partial x_2}&\cdots&\frac{\partial f_1(\vec x)}{\partial x_k}\\\frac{\partial f_2(\vec x)}{\partial x_1}&\frac{\partial f_2(\vec x)}{\partial x_2}&\cdots&\frac{\partial f_2(\vec x)}{\partial x_k}\\\vdots&\vdots&\ddots&\vdots\\\frac{\partial f_k(\vec x)}{\partial x_1}&\frac{\partial f_k(\vec x)}{\partial x_2}&\cdots&\frac{\partial f_k(\vec x)}{\partial x_k}\end{bmatrix}\).

Detta är ett linjärt ekvationssystem vars lösning ger oss

\[\vec h_n=\vec x_{n+1}-\vec x_n\\\implies\vec x_{n+1}=\vec x_n+\vec h_n\]

Obs! Jakobianen måste vara inverterbar för alla \(n\) för att systemet skall ha en lösning.

Stoppvillkor: \(||h_n||_2\lt\tau\)

Likt för Newtons metod för skalära problem konvergerar metoden kvadratiskt, dvs. \(||\vec h_n||_2\approx K||\vec h_{n-1}||_2^2\).

Föreläsning 5

Lösning av (stora) linjära ekvationssystem

Övertriangulära matriser

\[U=\begin{bmatrix}a_{11}&a_{12}&a_{13}&a_{14}\\0&a_{22}&a_{23}&a_{24}\\0&0&a_{33}&a_{34}\\0&0&0&a_{44}\end{bmatrix}\]

En kvadratisk matris där alla element under huvuddiagonalen är noll.

System på formen \(U\vec x=\vec b\) löses enkelt med bakåtsubstitution (lös sista raden och sätt in på raden ovan, osv.).

Undertriangulära matrisen

\[L=\begin{bmatrix}a_{11}&0&0&0\\a_{21}&a_{22}&0&0\\a_{31}&a_{32}&a_{33}&0\\a_{41}&a_{42}&a_{43}&a_{44}\end{bmatrix}\]

En kvadratisk matris där alla element ovanför huvuddiagonalen är noll.

System på formen \(L\vec x=\vec b\) löses enkelt med framåtsubstitution (lös första raden och sätt in på nästa rad, osv.).

LU-faktorisering

I många fall kan man faktorisera en matris så att \(A=LU\).

Att lösa \(LU\vec x=\vec b\) görs sedan i två steg:

  1. Låt \(U\vec x=\vec y\) och lös systemet \(L\vec y=\vec b\) med framåtsubstitution.
  2. Med \(\vec y\) från steg 1, lös systemet \(U\vec x=\vec y\) med bakåtsubstitution.

Matrisen \(U\) fås genom att reducera \(A\) till trappstegsform (Gausselimination).

\[\begin{align*}A=&\left.\begin{bmatrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\end{bmatrix}\right|\begin{bmatrix}1&0&0\\0&1&0\\0&0&1\end{bmatrix}\\&\left.\begin{bmatrix}a_{11}&a_{12}&a_{13}\\\red0&\tilde a_{22}&\tilde a_{23}\\\red0&\tilde a_{32}&\tilde a_{33}\end{bmatrix}\right|\begin{bmatrix}1&0&0\\\red{-\frac{a_{21}}{a_{11}}}&1&0\\\red{-\frac{a_{31}}{a_{11}}}&0&1\end{bmatrix}\\U=&\left.\begin{bmatrix}a_{11}&a_{12}&a_{13}\\0&\tilde a_{22}&\tilde a_{23}\\0&\red0&\hat a_{33}\end{bmatrix}\right|\begin{bmatrix}1&0&0\\-\frac{a_{21}}{a_{11}}&1&0\\-\frac{a_{31}}{a_{11}}\red{+\frac{\tilde a_{32}}{\tilde a_{22}}\frac{a_{21}}{a_{11}}}&\red{-\frac{\tilde a_{32}}{\tilde a_{22}}}&1\end{bmatrix}\end{align*}\]

Matrisen till höger fås genom att utföra samma radoperationer på enhetsmatrisen. Inversen av denna matris blir \(L\):

\[L=\begin{bmatrix}1&0&0\\\frac{a_{21}}{a_{11}}&1&0\\\frac{a_{31}}{a_{11}}&\frac{\tilde a_{32}}{\tilde a_{22}}&1\end{bmatrix}\]

Om radbyten krävs för att lösa systemet \(A\vec x=\vec b\) kan dessa beskrivas med en permutationsmatris \(P\). På så sätt kan det radbytta systemet skrivas som:

\[PA\vec x=P\vec b\]

Vi LU-faktoriserar nu den radbytta matrisen \(PA\). Sammantaget får vi

\[A\vec x=\vec b\iff PA\vec x=P\vec b\iff LU\vec x=P\vec b\]
Användningsområden för LU/PLU-faktorisering

Antag att vi vill lösa många system \(A\vec x_i=\vec b_i,i=1,2,\dots\).

Vi kan antingen lösa dessa ett i taget med Gausselimination, eller LU-faktorisera \(A\) & därefter lösa varje system med framåt/bakåt-substitution.

Komplexitet

Exempel:

För att beräkna \(A\vec x\), där \(A\in\mathbb R^{n\times n}\) och \(\vec x\in\mathbb R^n\), krävs \(2n^2-n\) operationer.

Exempel 2: För att lösa \(A\vec x=\vec b\) med Gausselimination, där \(A\in\mathbb R^{n\times n}\) och \(\vec b\in\mathbb R^n\), krävs totalt

\[N=\frac23n^3+\frac32n^2-\frac76n\]

flyttalsoperationer.

Jämförelse av algoritmer
Algoritm Komplexitet
Matris-vektor-multiplikation \(\mathcal O(n^2)\)
Matris-matris-multiplikation \(\mathcal O(n^3)\)
Gausselimination \(\mathcal O(n^3)\)
LU-faktorisering \(\mathcal O(n^3)\)
Framåt/bakåt-substitution \(\mathcal O(n^2)\)
Lösa \(A\vec x=\vec b\) där \(A\) är en gles bandmatris \(\mathcal O(n)\)

I komplexitetsanalys är man oftast bara intresserad av den dominerande termen, eftersom \(n\) är ett stort tal.

Glesa bandmatriser

En matris är gles om nästan alla element \(=0\). Ett sådant ekvationssystem kan lösas betydligt snabbare än ett system med fylld matris.

En gles bandmatris har nollskilda element endast i ett diagonalt band. Ett sådant system kan lösas med komplexitet \(n\).

LU-faktorisering (forts.)

Proceduren består av två steg:

  1. LU-faktorisering (Gausselimination): \(T_{LU}\propto n^3\)
  2. Lösa \(LU\vec x=\vec b\) med framåt- och bakåt-substitution: \(T_L\propto n^2,T_U\propto n^2\)

Steg 1 behöver bara utföras en gång. Därefter kan alla system

\[A\vec x_i=\vec b_i,i=1,2,\dots\]

lösas genom att helt sonika upprepa steg 2.

I Matlab kan LU-faktorisering utföras med funktionen lu: 

[L, U] = lu(A);
% eller
[L, U, P] = lu(A);

Konditionstal

Vid lösning av \(A\vec x=\vec b\), där \(\vec b\) innehåller fel (t.ex. mätfel, avrundningsfel), kommer detta fel att påerka lösningen \(\vec x\).

Om felet i \(\vec b\) är litet, men felet i \(\vec x\) blir stort, är problemet illakonditionerat.

Maxnorm

\[||A||_\infty=\max_i\sum_{j=1}^n|a_{ij}|="\text{maxsumma av absoluta radelement}"\\||\vec x||_\infty=\max_i|x_i|\]

För matrisers maxnorm gäller att

\[||A\vec x||_\infty\leq||A||_\infty||\vec x||_\infty\]

Konditionstal för matriser

Låt \(\vec x\) vara lösningen till \(A\vec x=\vec b\), och \(\vec x_\epsilon\) vara lösningen till \(A\vec x_\epsilon =\vec b_\epsilon\), där \(\vec b_\epsilon\) innehåller ett litet fel.

Relativt fel in: \(e_\text{in}=\frac{||\vec b-\vec b_\epsilon||_\infty}{||\vec b||_\infty}\)

Relativt fel ut: \(e_\text{ut}=\frac{||\vec x-\vec x_\epsilon||_\infty}{||\vec x||_\infty}\)

Matrisers konditionstal: (definition)

\[||A||_\infty||A^{-1}||_\infty\]

Konditionstalet anger värsta möjliga förstoringsfaktor för felet:

\[e_\text{ut}\leq\textcolor{#f32}{||A||_\infty||A^{-1}||_\infty}e_\text{in}\]

Vid lösning av ett system \(A\vec x=\vec b\), där \(\vec b\) innehåller en osäkerhet \(\Delta\vec b\), ges den fortplantade osäkerheten (maximala felet; gäller för samtliga komponenter i \(\vec x\)) av

\[||\Delta\vec x||_\infty\leq||A||_\infty||A^{-1}||_\infty\cdot\frac{||\Delta\vec b||_\infty}{||\vec b||_\infty}||\vec x||_\infty\]

Föreläsning 6

Numerisk integration

Målet: Approximera \(I=\int_a^bf(x)dx\)

Trapetsregeln

Approximerar funktionen som en rät linje \((a,f(a))\to(b,f(b))\), och integrerar under denna:

\[\int_a^bf(x)dx\approx(b-a)\cdot\frac{f(a)+f(b)}2\]
Sammansatta trapetsregeln

För att få en bättra approximation delas intervallet upp i delintervall. Med \(n\) delintervall fås \(n+1\) st integrationspunkter \(x_i=a+i\cdot h\), där steglängden \(h=\frac{b-a}n\). Regeln kan formuleras som:

\[\boxed{\int_a^bf(x)dx\approx h\cdot\left[\frac{f(x_0)}2+\left(\sum_{i=1}^{n-1}f(x_i)\right)+\frac{f(x_n)}2\right]}\]
Approximationsfel

Om \(I\) är integralens exakta värde och \(T_h(f)\) är approximationen med trapetsregeln (steglängd \(h\)), gäller för approximationsfelet (även kallat diskretiseringsfelet) \(e_h\)

\[e_h=|I-T_h(f)|\leq Ch^2\]

där \(C=\left|\frac{b-a}{12}f''(c)\right|\) och \(c\in(a,b)\).

Felet är av storleksordning \(h^2\) -- trapetsregeln har noggrannhetsordning 2.

Simpsons regel

Använder kvadratisk interpolation istället.

Basintervallet är \(2h\), eftersom det krävs 3 punkter för varje kvadratiska interpolation. Observera att antalet intervall (\(n\)) måsta vara jämnt.

\[\boxed{\int_a^bf(x)dx\approx\frac h3\left[f(x_0)+4\left(\sum_{i=1}^{n/2}f(x_{2i-1})\right)+2\left(\sum_{i=1}^{n/2-1}f(x_{2i})\right)+f(x_n)\right]}\]
Approximationsfel

Om \(I\) är integralens exakta värde och \(S_h(f)\) är approximationen med Simpsons regel (steglängd \(h\)), gäller för approximationsfelet \(e_h\) att

\[e_h=I-S_h(f)\leq Ch^4\]

där \(C=\left|\frac{b-a}{180}f^{(4)}(c)\right|\) och \(c\in(a,b)\).

Felet är av storleksordning \(h^4\) -- Simpsons regel har noggrannhetsordning 4.

Noggrannhetsordning - allmänt

Om en numerisk metod med steglängd \(h\) producerar en lösning med approximationsfelet

\[e_h\leq Kh^p\]

gäller då att metoden har noggrannhetsordning \(p\)**.

Om vi halverar steglängden får vi approximationsfelet \(e_{h/2}\). Observera följande förhållande:

\[\frac{e_h}{e_{h/2}}\approx\frac{Kh^p}{K\left(\frac h2\right)^p}=2^p\]

Kontroll av noggrannhetsordning (Trapetsregeln)

Då vi inte känner till den exakta lösningen kan vi inte beräkna felet. Istället kontrollerar vi noggrannhetsordningen experimentellt med hjälp av steglängdshalvering: \(T_h\),\(T_{h/2},T_{h/4},\dots\) beräknas.

När \(h\) successivt halveras skall kvoten

\[\frac{e_h}{e_{h/2}}\approx\frac{\left|T_{h/2}-T_h\right|}{\left|T_{h/4}-T_{h/2}\right|}\]

gå mot \(2^p=4\) (för Trapetsregeln).

Notera att \(\left|T_{h/2}-T_h\right|\) används som en approximation av felet \(e_h\)

Motsvarande metod kan användas för att kontrollera noggrannhetsordningen för Simpsons regel, men kvoten skall då gå mot \(2^4=16\).

Richardson-extrapolering

\[T^*(f)=\frac{4T_{h/2}(f)-T_h(f)}3\]

Felet i \(T^*\) är i storleksordningen \(h^4\), jämfört med \(T_{h/2}\) där felet är i stl.ordning \(h^2\).

Föreläsning 7

Numerisk derivering

Framåtdifferens

\[f'(x)\approx\frac{f(x+h)-f(x)}h=D_+f\]

Central differens

\[f'(x)\approx\frac{f(x+h)-f(x-h)}{2h}=D_cf\]

Bakåtdifferens

\[f'(x)\approx\frac{f(x)-f(x-h)}h=D_-f\]

Observera att \(D_c=\frac{D_++D_-}2\).

Approx. \(f''(x)\) (central differens)

\[f''(x)\approx\frac{f(x+h)-2f(x)+f(x-h)}{h^2}=D_{+-}f\]

Approximationsfel och noggrannhetsordning

Se slide 13-15 (F7 m. anteckningar)!!!

Numerisk lösning av ordinära diff.ekvationer (ODE:er)

Vi vill, givet en funktion \(f(t,y)\), samt \(t_0\) och \(y_0\) (begynnelsevärde), hitta en funktion \(y(t)\) som uppfyller

\[\begin{equation}\frac{dy}{dt}=f(t,y),\quad y(t_0)=y_0\end{equation}\]

En numerisk lösning ger oss värden \(y_i\approx y(t_i)\) som uppfyller (1). Tillsammans ger \((t_i,y_i)\) oss punkter som approximerar funktionen \(y=y(t)\).

Eulers metod (framåt och bakåt)

Eulers framåt

Vi approximerar \(y'(t_i)\approx\frac{y(t_{i+1})-y(t_i)}h\), där vi använder våra approximerande värden \(y_i\approx y(t_i),\ y_{i+1}\approx y(t_{i+1})\). Punkten \((t_0,y_0)\) och funktionen \(f(t,y)\) är givna.

Detta ger

\[\frac{y_{i+1}-y_i}h=f(t_i,y_i)\implies y_{i+1}=y_i+hf(t_i,y_i)\]

där \(h\) är steglängden och \(i=0,1,\dots,n\).

Eulers framåt är en explicit metod, vilket innebär att varje värde följer direkt av det föregående.

Eulers bakåt
\[y_{i+1}=y_i+hf(t_{i+1},y_{i+1}),\ i=0,1,\dots,n\]

Euler bakåt är en implicit metod, eftersom \(y_{i+1}\) även förekommer i HL. För att hitta \(y_{i+1}\) krävs att vi löser ekvationen.

Noggrannhetsordning

Felet för Eulers metoder är

\[e_h=\left|y(t_i)-y_i\right|\approx Kh\]

Metoderna har noggrannhetsordning 1, och felet minskar alltså linjärt i takt med att \(h\) minskar.

Felet kommer från approximationen av derivatan:

\[\left|y'(t_i)-\frac1h\left[y(t_i+h)-y(t_i)\right]\right|\approx\{\text{Taylorutveckla }y(t_i+h)\}\\\approx\left|y'(t_i)-\frac1h\left[\green{y(t_i)}+y'(t_i)\cdot(\red{t_i}+h\red{-t_i})+\frac{y''(t_i)}{2!}\cdot(\red{t_i}+h\red{-t_i})^2\green{-y(t_i)}\right]\right|\\\approx\left|\red{y'(t_i)}-\frac1{\blue h}\left[\red{y'(t_i)}\cdot \blue h+\frac{y''(t_i)}2\cdot \blue h^2\right]\right|\\\approx\frac h2\left|y''(t_i)\right|\]
Kontroll av noggrannhetsordning:

Görs experimentellt med steglängdshalvering. Kvoten

\[\frac{\left|y_{h/2}-y_h\right|}{\left|y_{h/4}-y_{h/2}\right|}\to2^p,\quad p=1\]

där \(y_h,y_{h/2},y_{h/4}\) är y-värden, beräknade med olika steglängder, för något \(t=T\).

Stabilitet

Det finns krav på \(h\) för att den numeriska lösningen skall vara stabil.

Implicita metoder (t.ex. Euler bakåt) är ofta mer stabila.

För att avgöra för vilka problem och steglängder en metod är stabil används följande testproblem:

\[y'=-\lambda y,\ \lambda\gt0,\ y(0)=1\]

Detta problem har den analytiska lösningen \(y(t)=e^{-\lambda t}\), vilket innebär att \(y\to0\) när \(t\to\infty\).

Eulers framåt: \(y_{i+1}=(1-\lambda h)y_i\implies y_i=(1-\lambda h)^i\)

  • För att \(\lim_{i\to\infty}y_i=0\) måste \(h\lt\frac2\lambda\).

Exempel:

\(\frac{dy}{dt}=\sin(3t)-2y,\ y(0)=1.2,\ t\in[0,8]\)

har lösningen \(y(t)=\frac{93}{65}e^{-2t}-\frac3{13}\cos(3t)+\frac2{13}\sin(3t)\)

I detta fallet är \(\lambda=2\), varför kravet på steglängd är \(h\lt1\).

Eulers bakåt: \(y_{i+1}=\frac{y_i}{1+\lambda h}\implies y_i=\frac1{(1+\lambda h)^i}\)

  • Metoden är stabil för alla \(h\).

Val av steglängd

Steglängdshalvering kan användas: \(h\) halveras tills \(e_h=\left|y_{h/2}-y_h\right|\) är tillräckligt liten för önskad noggrannhet i svaret.

Föreläsning 8

Runge-Kutta 4 (RK4)

\[\boxed{\begin{aligned}y_{i+1}&=y_i+\frac h6\left(k_1+2k_2+2k_3+k_4\right)\\k_1&=f(t_i,\ y_i)\\k_2&=f(t_i+\frac h2,\ y_i+\frac h2k_1)\\k_3&=f(t_i+\frac h2,\ y_i+\frac h2k_2)\\k_4&=f(t_{i+1},\ y_i+hk_3)\\t_{i+1}&=t_i+h\\i&=0,1,2,\dots\end{aligned}}\]

Noggrannhetsordning 4, dvs. felet \(e\propto h^4\). Stabilitetsanalys på testproblemet ger \(h\lt\frac{2.78}\lambda\).

Funktioner i Matlab

I Matlab finns flera olika ODE-lösare inbyggda, däribland ode45:

% fun har två inparametrar och motsvarar högerledet f(t, y).
% [t0 T] anger start- och sluttid
% y0 är begynnelsevillkoret
% [t, y] är två kolumnvektorer med t- respektive y-värden
% 
% Obs! ode45 använder adaptiv steglängd och punkterna är
%      därför inte ekvidistant placerade längs t-axeln.
[t, y] = ode45(fun, [t0 T], y0);

System av ODE:er

Standardform:

\[\frac{d\vec u}{dt}=\vec u'=\vec f(t,\ \vec u),\quad\vec u(t_0)=\vec c\]

Exempel: SIR

\[\vec u(t)=\begin{bmatrix}S(t)\\I(t)\\R(t)\end{bmatrix},\quad\vec f(t,\ \vec u)=\begin{bmatrix}-\beta SI-vS\\\beta SI-\gamma I\\\gamma I+vS\end{bmatrix},\quad\vec c=\begin{bmatrix}1-I_0\\I_0\\0\end{bmatrix}\]

Samma numeriska metoder som för skalära ODE:er kan användas för system -- endast notation behöver bytas (t.ex. \(y\) ist. f. \(\vec u\)).

Högre ordningens ODE

\[\begin{equation}y''=f(t,\ y,\ y'),\quad y(t_0)=c_1,\quad y'(t_0)=c_2\end{equation}\]

där \(f(t,\ y,\ y'),\ t_0,\ c_1,\ c_2\) är givna. Notera att två begynnelsevillkor krävs.

Ett begynnelsevärdeproblem på denna form kan skrivas om som ett system på standardform:

Vi inför två nya funktioner:

\[a(t)=y(t),\ b(t)=y'(t)\]

Derivering av dessa ger

\[\begin{aligned}\green{a'(t)}&=y'(t)\green{=b(t)}\\\blue{b'(t)}&=y''(t)=\{\text{ekv. (2)}\}=f(t,\ y,\ y')\blue{=f(t,\ a,\ b)}\end{aligned}\]

Det kompletta systemet ser ut som följer:

\[\vec u(t)=\begin{bmatrix}a(t)\\b(t)\end{bmatrix},\quad \vec f(t)=\begin{bmatrix}\green{b(t)}\\\blue{f(t,\ a,\ b)}\end{bmatrix},\quad \vec u(0)=\vec c=\begin{bmatrix}c_1\\c_2\end{bmatrix}\]

Föreläsning 9

Begynnelsevärdesproblem, fortsättning

Styva problem

I vissa system av ODE:er varierar någon komponent snabbt, medan andra varierar långsamt men blir dominerande när den snabba komponenten avtagit.

Explicita metoder är ineffektiva vid styva problem, eftersom den snabba komponenten ställer krav på kort steglängd för stabilitet.

Randvärdesproblem (RVP)

Andra (eller högre) ordningens ODE där lösningens värde är givet i två (rand)punkter.

\[y''+g(x)y'+h(x)y=f(x),\quad y(a)=\alpha,\quad y(b)=\beta,\quad x\in[a,b]\]

Finita differensmetoden

Låt \(y=y(x)\) och \(y_i\approx y(x_i)\). Då kan vi (likt tidigare) approximera derivatorna med centraldifferenser:

\[y'(x_i)\approx\frac{y_{i+1}-y_{i-1}}{2h}\\y''(x_i)\approx\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}\]

\(y''\)-approximationen har noggrannhetsordning 2.

Derivatan i problemet ersätts med approximationen ovan. Därefter samlas termer framför \(\red{y_{i-1}},\green{y_i},\blue{y_{i+1}}\); övriga termer flyttas till HL. Slutligen kan de tre samlade koefficienterna placeras ut längs diagonalerna i en tridiagonal matris, där huvuddiagonalen består av \(\green{y_i}\)-koefficienten.

Exempel:

\[\begin{bmatrix}\green{1-\frac2{h^2}}&\blue{\frac1{h^2}}&0\\\red{\frac1{h^2}}&\green{1-\frac2{h^2}}&\blue{\frac1{h^2}}\\0&\red{\frac1{h^2}}&\green{1-\frac2{h^2}}\end{bmatrix}\begin{bmatrix}y_1\\y_2\\y_3\end{bmatrix}=\begin{bmatrix}\textcolor{#fc0}{f(x_1)-\frac\alpha{h^2}}\\\textcolor{#fc0}{f(x_2)}\\\textcolor{#fc0}{f(x_3)-\frac\beta{h^2}}\end{bmatrix}\]
Noggrannhet

Noggrannheten i lösningen ges av noggrannheten hos derivataapproximationen. I fallet ovan: 2:a ordningens central differens ⇒ noggrannhetsordning 2.

Föreläsning 10

Randvillkor

Det finns två sorters randvillkor:

  • Dirichletvillkor: \(y(b)=\beta\)
  • Neumannvillkor: \(y'(a)=\alpha\)
    • Derivatan approximeras med en (enkelsidig) differensapproximation.

Neumannvillkor

Antag att vi har

\[y''(x)+g(x)\cdot y'(x)+y(x)=f(x),\quad y'(a)=\alpha,\quad y(b)=\beta\]

Detta kan lösas på två sätt:

  1. Låt \(y_0\) finnas med bland de obekanta -- den extra ekvationen/raden blir \(-y_0+y_1=h\alpha\) (approximation av randvillkoret).
  2. I ekvationen för första inre punkten, ersätt \(y_0\) med \(y_1-h\alpha\). Detta ger ingen extra obekant, men kräver att man i slutet beräknar \(y_0=y_1-h\alpha\).

Noggrannhet

Noggrannheten avgörs av noggrannheten i derivata-approximationerna som används.

Dirichletvillkor

Om endast centrala differenser används: fel i stl.ordning \(h^2\). Hela metoden har då noggr.ordn. 2.

Neumannvillkor

Om även 1:a ordningens differenser används: fel i stl.ordning \(h\). Hela metoden har då noggr.ordn. 1.

Kontroll av noggrannhetsordning

Studera lösningen som fås i en fix punkt \(x=\hat x\). Bilda kvoterna

\[\frac{y_{h/2}-y_h}{y_{h/4}-y_{h/2}},\frac{y_{h/4}-y_{h/2}}{y_{h/8}-y_{h/4}},\dots\to2^p\]

Olinjära randvärdesproblem

Exempel:

\[\frac{d^2y}{dx^2}-y+y^2=0,\quad u(a)=\alpha,\quad u(b)=\beta,\quad x\in[a,b]\]

Vi tillämpar finita differensmetoden och får ett system av olinjära ekvationer. Detta system kan lösas med Newtons metod.

Felfortplantning och störningsanalys

Notation

Notation Betydelse Notation Betydelse
\(\tilde x=x+e_x\) Värde med fel \(\tilde y=y+e_y\) Utdata med fel
\(x\) Exakt värde (okänt) \(y\) Verkligt värde (okänt)
\(e_x\) Felet i indata (okänt) \(e_y\) Felet i utdata (okänt)
\(E_x\) Absolut felgräns (känd) \(E_y\) Absolut felgräns (okänd)
\(R_x=E_x/x\) Relativa felet i \(\tilde x\) \(R_y=E_y/\tilde y\) Relativa felet i \(y\)

Vi vet att \(|e_x|\leq E_x\).

Att uppskatta felgränsen \(E_y\)

Två sätt:

  1. Om vi vet explicit hur \(F\) beror av indata (dvs. \(y=F(x)\)), med felfortplantningsformeln (\(E_y\leq E_x|F'(\tilde x)|\)).
  2. Om vi inte vet explicithur \(F\) beror av indata, med experimentell störningsanalys.
Felfortplantningsformeln
\[\begin{aligned}e_y&=\tilde y-y\approx F(\tilde x)-y^*=F(\tilde x)-(F(\tilde x)+(x-\tilde x)F'(\tilde x))\\&=(\tilde x-x)F'(\tilde x)=e_xF'(\tilde x)\implies\{|e_x|\leq E_x\}\implies\boxed{E_y\leq E_x\left|F'(\tilde x)\right|}\end{aligned}\]

Om \(y=F(x_1,\ x_2)\):

\[E_y\leq E_{x_1}\left|\frac{\partial F}{\partial x_1}\right|+E_{x_2}\left|\frac{\partial F}{\partial x_2}\right|\]

\(E_y\) är gränsen för absolutfelet i utdata -- dvs. \(y=\tilde y\pm E_y\).

Experimentell störningsanalys

Exempel: Antag att vi vill lösa ett problem \(\tilde y=F(\tilde x;\ \tilde\alpha,\ \tilde\beta)\).

Stör en variabel i taget (uppåt) och uppskatta \(E_y\) som

\[E_y\approx|\tilde y_1-\tilde y|+|\tilde y_2-\tilde y|+|\tilde y_3-\tilde y|\]

där

\[\begin{aligned}\tilde y_1&=F(\tilde x\red{+E_x};\ \tilde\alpha,\ \tilde\beta)\\\tilde y_2&=F(\tilde x;\ \tilde\alpha\red{+E_\alpha},\ \tilde\beta)\\\tilde y_3&=F(\tilde x;\ \tilde\alpha,\ \tilde\beta\red{+E_\beta})\end{aligned}\]

Om störningarna t.ex. är stora, eller om \(F\) varierar kraftigt, kan det vara intressant att beräkna \(\tilde y\) för alla kombinationer av störningar, dvs. \(3^n\) stycken, där \(n\) är antalet variabler.