TMV035 Analysis and Linear Algebra B

Ordinära differentialekvationer 2

1.1 Exponentialfunktionen

Vi erinrar oss att exponentialfunktionen uppfyller

D exp(x) = exp(x), exp(0) = 1.

Det betyder att u(x) = exp(x) är lösning till begynnelsevärdesproblemet

\begin{matrix}u′(x) = u(x), x ∈ [0, b] (b > 0) \\ u(0) = 1. \end{matrix}

Man kan visa att lösningen till detta problem är unik (vi återkommer till detta) så att u(x) = exp(x) är den enda lösningen. Vi ska använda detta begynnelsevärdesproblem för att konstruera u(x) = exp(x) för x ≥ 0. För x < 0 kan vi sedan definiera exp(x) = 1/ exp(−x).

Mer allmänt kan vi betrakta begynnelsvärdesproblemet

\begin{matrix}u′(x) = cu(x), x ∈ [0, b] \\ u(0) = u_{0}. \end{matrix}

Det är lätt kolla att lösningen är u(x) = u_{0} exp(cx). Notera att u är växande om c > 0, avtagande om c < 0 och konstant om c = 0.

Ännu mer allmänt kan vi betrakta

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

Det är lätt kolla att lösningen är u(x) = u_{a} exp(c(x − a)).

Vi konstruerar lösningen till (1) med hjälp av Eulers metod. 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}) = cu(x_{i−1}),

vilket leder till

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

Vi beräknar nu en approximativ lösning enligt

U(x_{0}) = u_{a}
U(x_{i}) = U(x_{i−1}) + hcU(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}.

Konvergens

Man kan visa att U(x) konvergerar mot en unik lösning u(x) till (1) 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 (1), dvs u är deriverbar med u′ = cu och u(a) = u_{a}.

4. Vi visar att lösningen till (1) är unik.

Detta bevis är ganska långt och lite för svårt och vi nöjer oss med att nämna dessa steg. Vi har nu konstruerat en ny funktion u(x) som kallas exp(x).

Implementering i Matlab

Algoritmen är

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}) + hcU(x_{i−1}). \end{matrix}

Övning 1. Skriv ett program myexp.m med anropet [x,U]=myexp(c,I,ua,h) som löser begynnelsevärdesproblemet (1). Du skall använda programskalet

myexp.m.

Börja med c = 1, I = [0, 1], u_{a} = 1. 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}. Obs att Matlabs exp är inte en exakt lösning, det är också en approximativ lösning, dock beräknad med hög noggrannhet.

Plotta till sist flera lösningar med olika värden på c, positiva och negativa.

1.2 Trigonometriska funktioner

Vi vet att de trigonometriska funktionerna uppfyller

D sin(t) = cos(t), D^{2} sin(t) = − sin(t), sin(0) = 0,
D cos(t) = − sin(t), D^{2} sin(t) = − cos(t), cos(0) = 1.

(Notera att vi nu går över till att skriva den oberoende variabeln som t istället för x. Anledningen till detta är att den ofta är tiden i ett begynnelsevärdesproblem.) Det betyder att u(t) = sin(t) löser begynnelsevärdesproblemet

\begin{matrix}u′′(t) = −u(t), t ∈ [0, b], (b > 0) \\ u(0) = 0, u′(0) = 1, \end{matrix}

och att u(t) = cos(t) löser begynnelsevärdesproblemet

\begin{matrix}u′′(t) = −u(t), t ∈ [0, b], (b > 0) \\ u(0) = 1, u′(0) = 0, \end{matrix}

Mer allmänt har vi

\begin{matrix}u′′(t) = −c^{2}u(t), t ∈ [0, b], (b > 0) \\ u(0) = u_{0}, u′(0) = u_{1}, \end{matrix}

med lösningen

u(t) = u_{0} cos(ct) + \frac{u_{1}}{c} sin(ct).  (2)

Verifiera detta genom att derivera och genom att sätta in t = 0!

Denna differentialekvation är av andra ordningen. För att kunna använda vår algoritm skriver vi om den som ett system av två differentialkevationer av första ordningen. Vi inför två nya variabler

w_{1} = u, w_{2} = u′.

Vi deriverar

w′_{1} = u′ = w_{2},
w′_{2} = u′′ = −c^{2}u = −c^{2}w_{1}.

Begynnelsevillkoren blir w_{1}(0) = u_{0}, w_{2}(0) = u_{1}. Vi ser att w_{1} och w_{2} löser begynnelsevärdesproblemet

\{ \begin{matrix}w′_{1} = w_{2}, \\ w′_{2} = −c^{2}w_{1}, \end{matrix} \{ \begin{matrix}w_{1}(0) = u_{0}, \\ w_{2}(0) = u_{1}. \end{matrix}

Vi skriver detta på matrisform:

[\begin{matrix}w_{1} \\ w_{2} \end{matrix}]′ = [\begin{matrix}w′_{1} \\ w′_{2} \end{matrix}] = [\begin{matrix}w_{2} \\ −c^{2}w_{1} \end{matrix}] = [\begin{matrix}0 & 1 \\ −c^{2} & 0 \end{matrix}] [\begin{matrix}w_{1} \\ w_{2} \end{matrix}] , [\begin{matrix}w_{1}(0) \\ w_{2}(0) \end{matrix}] = [\begin{matrix}u_{0} \\ u_{1} \end{matrix}] .

Med beteckningarna

A = [\begin{matrix}0 & 1 \\ −c^{2} & 0 \end{matrix}] , w = [\begin{matrix}w_{1} \\ w_{2} \end{matrix}] , w_{0} = [\begin{matrix}u_{0} \\ u_{1} \end{matrix}] ,

får vi slutligen

\begin{matrix}w′(t) = Aw(t), t ∈ [0, b], \\ w(0) = w_{0}. \end{matrix}  (3)

Vår algoritm (Eulers metod) blir nu

initiera: \{ \begin{matrix}t_{0} = 0 \\ W (t_{0}) = w_{0} \end{matrix}

uppdatera: \{ \begin{matrix}while t_{i} < b \\ t_{i} = t_{i−1} + h \\ W (t_{i}) = W (t_{i−1}) + hAW (t_{i−1}). \end{matrix}

Övning 2. Skriv ett program mytrig.m med anropet [t,W]=mytrig(c,I,w0,h) som löser begynnelsevärdesproblemet (3). Du skall använda programskalet

mytrig.m.

I Matlab sparas W (t_{i}) som i:te kolonnen W(:,i) i matrisen W, den blir av typ 2 × (N + 1). Matrisen skapas kolonnvis. Till exempel blir initieringen W(:,1)=w0. När beräkningen är klar bör du transponera matriserna t och W bland annat för att plottning ska fungera smidigt.

Tag först c = 1 och w_{0} = [\begin{matrix}0 \\ 1 \end{matrix}].

Det finns nu två sätt att plotta lösningen.

1. Vi kan plotta de två lösningskurvorna y = w_{1}(t), y = w_{2}(t) med kom-

     mandot plot(t,W).

2. Vi kan plotta w_{2} mot w_{1} med kommandot plot(W(:,1),W(:,2)). Sådana kurvor med olika begynnelsepunkter bildar ett fasporträtt för differentialekvationssystemet.

1.3 Hyperboliska funktioner

Vi vet att de hyperboliska funktionerna uppfyller

D sinh(t) = cosh(t), D^{2} sin(t) = sinh(t), sinh(0) = 0,
D cosh(t) = sinh(t), D^{2} sinh(t) = cosh(t), cosh(0) = 1.

Det betyder att både u(t) = sinh(t) och u(t) = cosh(t) uppfyller differentialekvationen

u′′(t) = u(t).

Övning 3. Visa hur denna ekvation skrivs om till ett system av två differentialekvationer av första ordningen och hur vi får ett begynnelsevärdesproblem av typen (3) fast med matrisen

A = [\begin{matrix}0 & 1 \\ 1 & 0 \end{matrix}] .

Gör motsvarande ändring i programmet mytrig.m och använd detta till att plotta de två kurvorna y = sinh(t), y = cosh(t).

2006-11-16 /stig