TMV035 Analysis and Linear Algebra B

Ordinära differentialekvationer 1

I en serie studioövningar ska vi lära oss hur man konstruerar och beräknar lösningar till ordinära differentialekvationer av typen

u′(x) = f(x, u(x))

Vi börjar i denna övning med det enklaste fallet då f = f(x) är en funktion som bara beror på x. I senare övningar låter vi f bero också på u och till sist låter vi u och f vara vektorfunktioner, dvs vi löser system av ordinära differentialekvationer.

1.1 Primitiv funktion

Antag att f : I → R är en funktion. Om funktionen u är deriverbar och uppfyller

u′(x) = f(x), ∀x ∈ I,  (1)

så säger vi att u är en primitiv funktion till f. Vi ska konstruera primitiva funktioner. Vi börjar med att notera att primitiv funktion är väsentligen unik. För om två funktioner u och v är primitiva funktioner till samma funktion f så gäller

(u − v)′(x) = u′(x) − v′(x) = f(x) − f(x) = 0, ∀x ∈ I.

Härav följer att

u − v = konstant

dvs

u(x) = v(x) + C.

För att bestämma konstanten kan man kräva att funktionen skall vara lika med ett bestämt värde i en bestämd punkt: u(a) = u_{a}. Vi kan nu precisera vår uppgift. Givet en kontinuerlig funktion f : [a, b] → R och en konstant u_{a} söker vi en deriverbar funktion u sådan att

\begin{matrix}u′(x) = f(x), x ∈ [a, b], \\ u(a) = u_{a}. \end{matrix}  (2)

Eftersom vi anger funktionens värde i intervallets vänstra ändpunkt så kallas problemet (2) ett begynnelsevärdesproblem.

Exempel 1. Om man känner ett par av en funktion och dess derivata så kan man lösa motsvarande begynnelsevärdesproblem. Till exempel känner vi till följande derivata:

Dx^{m} = mx^{m−1} (m ≥ 0)

Detta leder till

D \frac{x^{m}}{m} = x^{m−1} (m > 0)

dvs, med r = m − 1,

D \frac{x^{r+1}}{r + 1} = x^{r} (r > −1)

Alltså kan vi säga att problemet

u′(x) = x^{r} (r > −1)

har lösningen

u(x) = \frac{x^{r+1}}{r + 1} + C.

Vi bestämmer konstanten, till exemel, genom begynnelsevillkoret

u(0) = 0

vilket leder till C = 0 och u(x) = \frac{x^{r+1}}{r+1} . Mer allmänt: villkoret u(0) = u_{0} leder till C = u_{0} och u(x) = \frac{x^{r+1}}{r+1} + u_{0}. Alltså har begynnelsevärdesproblemet

\begin{matrix}u′(x) = x^{r}, x ∈ [0, ∞), \\ u(0) = u_{0}, \end{matrix}

den unika lösningen

u(x) = \frac{x^{r+1}}{r + 1} + u_{0}.

Denna metod fungerar bara om vi redan känner ett funktion/derivata-par. Vi ska nu beskriva hur man kan konstruera en primitiv funktion för varje kontinuerlig funktion f, dvs konstruera en unik lösning till begynnelsevärdesproblemet (2). Vi börjar med att dela in intervallet [a, b] i N stycken delintervall av längden h = (b − a)/N:

a = x_{0} < x_{1} < x_{2} < … x_{i−1} < x_{i} < ⋯ < x_{N−1} < x_{N} = b,
x_{i} = a + hi, h = (b − a)/N = x_{i} − x_{i−1}.

Sedan utgår vi från approximationen

\frac{u(x_{i}) − u(x_{i−1})}{h} ≈ u′(x_{i−1}) = f(x_{i−1})

vilket leder till

u(x_{i}) ≈ u(x_{i−1}) + hf(x_{i−1}).

Vi beräknar nu en approximativ lösning enligt

U(x_{0}) = u_{a}
U(x_{i}) = U(x_{i−1}) + hf(x_{i−1}).

Genom att förbinda punkterna (x_{i}, U(x_{i})) med räta linjer får vi en graf och funktionen U(x) blir definierad också mellan beräkningsnoderna x_{i}. Rita figur! Denna algoritm kallas Eulers metod.

Konvergens

Man kan visa att U(x) konvergerar mot en unik lösning u(x) till (2) då antalet delintervall N → ∞. Beviset är upplagt enligt samma princip som för bisektionsalgoritmen och fixpunktsiterationen.

Vi delar intervallet [a, b] med hjälp av halvering. Vi låter n beteckna antalet halveringar, antalet delintervall blir N = 2^{n}, och för varje n får vi en approximativ lösning U_{n}(x):

n = 1, N = 2^{1} = 2, U_{1}(x)

n = 2, N = 2^{2} = 4, U_{2}(x)

n = 3, N = 2^{3} = 8, U_{3}(x)

n = 4, N = 2^{4} = 16, U_{4}(x) och så vidare.

Beviset går sedan enligt följande.

1. Algoritmen Eulers metod ger en följd av approximativa lösningar {U_{n}(x)}_{n=1}^{∞}.

2. Vi visar att U_{n}(x) är en Cauchy-följd, dvs U_{n}(x) − U_{m}(x) → 0 då m, n → ∞. Vi får alltså ett reellt tal (decimalutveckling) u(x) = lim_{n∞} U_{n}(x).

3. Vi visar att u löser (2), dvs u är deriverbar med u′ = f och u(a) = u_{a}.

4. Vi visar att lösningen till (2) är unik. (Det har vi redan gjort ovan.)

Detta bevis är ganska långt och lite för svårt och vi nöjer oss med att nämna dessa steg.

Implementering i Matlab

Detta är lätt programmera i Matlab. Algoritmen är följande. Eftersom vi behöver både x_{i} och U(x_{i}) för att till exempel plotta funktionen så måste algoritmen räkna ut även noderna x_{i}.

initiera: \{ \begin{matrix}x_{0} = a \\ U(x_{0}) = u_{a} \end{matrix}

uppdatera: \{ \begin{matrix}while x_{i} < b \\ x_{i} = x_{i−1} + h \\ U(x_{i}) = U(x_{i−1}) + hf(x_{i−1}). \end{matrix}

I Matlab kan index inte börja med 0, så att till exempel initieringen av x skrivs x(1)=a. Likaså måste U(x_{i}) representeras av en vektor U med elementen U(i). Man kan inte skriva U(x(i)).

Övning 1. Skriv ett program myprim.m med anropet [x,U]=myprim(f,I,ua,h) som löser begynnelsevärdesproblemet (2). Du skall använda programskalet myprim.m. In- och ut-variablerna förklaras i programskalet.

Testa programmet på följande. Lös först begynnelsevärdesproblemet analytiskt (dvs som en formel med penna och papper). Plotta både den analytiska lösningen u och den approximativa lösningen U i samma figur. Använd stort steg h, till exempel, h = 10^{−1}, när du först skriver programmet så blir det lättare att se vad som händer. Tag sedan litet h, till exempel, h = 10^{−3}.

(a) \begin{matrix}u′(x) = x^{2}, x ∈ [0, 2], \\ u(0) = 3, \end{matrix}

(b) \begin{matrix}u′(x) = x^{2}, x ∈ [1, 4], \\ u(1) = 3, \end{matrix}

(c) \begin{matrix}u′(x) = cos(2x), x ∈ [0, π], \\ u(0) = 0, \end{matrix}

(d) \begin{matrix}u′(x) = 1/x, x ∈ [1, 10], \\ u(1) = 0, \end{matrix}

Obs att naturliga logarimen ln(x) heter log(x) i Matlab.

Vi får en noggrannare approximation om vi utgår från

\frac{u(x_{i}) − u(x_{i−1})}{h} ≈ u′(x̂_{i}) = f(x̂_{i}), x̂_{i} = (x_{i} + x_{i−1})/2,

vilket leder till

u(x_{i}) ≈ u(x_{i−1}) + hf(x̂_{i−1}).

Vi beräknar nu en approximativ lösning enligt U(x_{0}) = u_{a} U(x_{i}) = U(x_{i−1}) + hf(x̂_{i}).

Denna algoritm kallas mittpunktsmetoden. Övning 2. Modifiera programmet myprim.m så att det använder mittpunktsmetoden istället. Prova med samma exempel som ovan.

2006-11-14 /stig