TMV035 Analysis and Linear Algebra B^{1}

Ordinära differentialekvationer 5

Läraren går igenom avsnitt 1.1 och 1.2. Sedan gör du övningarna med penna och papper och med Matlab.

1.1 Allmän och speciell ekvation

Vi har löst det allmänna begynnelsevärdesproblemet

\begin{matrix}u′(t) = f(t, u(t)), t ∈ [a, b], \\ u(a) = u_{a}, \end{matrix}  (1)

med Eulers metod (framlänges Euler)

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

Man kan visa att om f är Lipschitzkontinuerlig, så konvergerar U(t) mot en unik lösning u(t) till (1) då antalet delintervall N → ∞, åtminstone på något intervall [a, b]. 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}(t):

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

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

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

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

och så vidare.

Beviset går sedan enligt följande.

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

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

3. Vi visar att u löser (1), dvs u är deriverbar med u′ = f(t, u) 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(t). I allmänhet kan man inte skriva ned u med en formel (analytisk lösning) men den kan beräknas med godtycklig noggrannhet med Eulers metod.

Vi ska nu studera några speciella typer av differentialekvationer som kan lösas analytiskt, dvs med en formel. Det är:

• separabel ekvation: u′(x) = h(x)/g(u(x));

• linjär ekvation av första ordningen: u′(x) + a(x)u(x) = f(x)

• linjär ekvation av andra ordningen med konstanta koefficienter: u′′(x)+bu′(x)+cu(x) = f(x)

• system av första ordningen med konstanta koefficienter: u′(x) = Au(x) (i ALA-c).

12006-11-28 /stig

1.2 Separabel ekvation

(Adams kap. 17.1.) Vi betraktar nu ett begynnelsevärdesproblem av formen

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

Differentialekvationen är av formen u′(x) = f(x, u(x)) med funktionen f av den speciella formen

f(x, s) = h(x)/g(s).

En sådan ekvation kallas separabel för vi kan separera variablerna x och u:

g(u(x))u′(x) = h(x)

Om vi kan gissa primitiva funktioner G och H till g och h så blir ekvationen

D(G(u(x))) = g(u(x))u′(x) = h(x) = DH(x)

dvs

D(G(u(x))) = DH(x)

och därmed

G(u(x)) = H(x) + C.

Konstanten C bestäms med hjälp av begynnelsevärdet i (2):

C = G(u_{a}) − H(a)

så att

G(u(x)) = H(x) + G(u_{a}) − H(a)  (3)

Om vi har tur kan vi sedan lösa ut u(x):

u(x) = G^{−1}(H(x) + G(u_{a}) − H(a)).

Det är praktiskt att skriva räkningarna så här. Utgå från (2). Separera variablerna: g(u(x))u′(x) = h(x)

g(u)\frac{du}{dx} = h(x)

g(u) du = h(x) dx

Integrera:

∫_{u_{a}}^{u(x)} g(s) ds = ∫_{a}^{x} h(s) ds
[G(s)]_{s=u_{a}}^{s=u(x)} = [H(s)]_{s=a}^{s=x}
G(u(x)) − G(u_{a}) = H(x) − H(a)
G(u(x)) = H(x) + G(u_{a}) − H(a)

vilket är samma som (3). Sedan löser man ut u(x) ur denna ekvation (om man kan) och då är det klart.

Exempel 1.

u′(x) = −xu(x), x > 0,
u(0) = 2.

Separera variablerna och integrera:

\frac{1}{u} \frac{du}{dx} = −x

\frac{1}{u}du = −x dx

∫_{2}^{u(x)} \frac{ds}{s} = − ∫_{0}^{x} s ds

[ ln(s)]_{s=2}^{s=u(x)} = −[\frac{1}{2}s^{2}]_{s=0}^{s=x}

ln(u(x)) − ln(2) = −\frac{1}{2}x^{2}

ln(u(x)) = ln(2) − \frac{1}{2}x^{2}

Lös ut u(x):

u(x) = exp(ln(2) − \frac{1}{2}x^{2}) = exp(ln(2)) exp(−\frac{1}{2}x^{2}) = 2 exp(−\frac{1}{2}x^{2})

Med Matlab: funktionsfil funk.m

function y=funk(t,u)
y=-t*u;

Kommandorad

>> [x,U]=myode(@funk,[0,3],2,1e-3);
>> plot(x,U)

1.3 Exercises

The only kind of nonlinear differential equations that we solve analytically is the so-called separable differential equations (including autonomous equations). Solution method: separate the variables and integrate.

Problem 1.1. Find analytical solution formulas for the following initial value problems. In each case sketch the graphs of the solutions and determine the half-life. See: P. Atkins and L. Jones, Chemical Principles. The Quest for Insight. Freeman, New York, second edition, 2002, pp. 698– 706. Then solve the equation numerically with your Matlab program myode and plot the solution.

(a) First order decay rate:

\begin{matrix}u′ = −ku, t > 0, (k > 0, u_{0} > 0) \\ u(0) = u_{0}. \end{matrix}

(b) Second order decay rate:

\begin{matrix}u′ = −ku^{2}, t > 0, (k > 0, u_{0} > 0) \\ u(0) = u_{0}. \end{matrix}

(c) Third order decay rate:

\begin{matrix}u′ = −ku^{3}, t > 0, (k > 0, u_{0} > 0) \\ u(0) = u_{0}. \end{matrix}

Problem 1.2. Find analytical solution formulas for the following initial value problems. Sketch the graphs.

(a) First order increase rate:

\begin{matrix}u′ = ku, t > 0, (k > 0, u_{0} > 0) \\ u(0) = u_{0}. \end{matrix}

(b) Second order increase rate:

\begin{matrix}u′ = ku^{2}, t > 0, (k > 0, u_{0} > 0) \\ u(0) = u_{0}. \end{matrix}

(c) Third order increase rate:

\begin{matrix}u′ = ku^{3}, t > 0, (k > 0, u_{0} > 0) \\ u(0) = u_{0}. \end{matrix}

(d) The logistic equation. The equation u′ = ku leads to exponential growth, u(t) = u_{0} exp(kt). In practice the exponential growth stops when u becomes large. One way of modeling this is to replace the rate ku by the rate ku(1 − u/M) which behaves like ku when u is small, but which is small when u approaches the large number M. We get the so-called “logistic equation” (with M = 1 for simplicity):

\begin{matrix}u′ = ku(1 − u), t > 0, (k > 0, u_{0} > 0) \\ u(0) = u_{0}. \end{matrix}

(e)

\begin{matrix}u′ = −tu^{2}, t > 0, (u_{0} > 0) \\ u(0) = u_{0}. \end{matrix}

1.4 Answers and solutions

1.1. Reaction of order 1 (decay rate of order 1):

\{ \begin{matrix}u′ = −ku \\ u(0) = u_{0} \end{matrix}

u(t) = u_{0}e^{−kt}

The half-life T_{1/2} is given by

u(T_{1/2}) = u_{0}e^{−kT_{1/2}} = \frac{1}{2}u_{0},

which leads to

T_{1/2} = \frac{log(2)}{k}.

Reaction of order n > 1 (decay rate of order n > 1):

\{ \begin{matrix}u′ = −ku^{n} \\ u(0) = u_{0} \end{matrix}

\frac{du}{u^{n}} = −k dt

∫_{u_{0}}^{u(T)} u^{−n} du = − ∫_{0}^{T} k dt

[\frac{u^{−n+1}}{−n + 1}]_{u_{0}}^{u(T)} = −kT

u(T)^{−n+1} − u_{0}^{−n+1} = (n − 1)kT

\frac{1}{u(T)^{n−1}} = \frac{1}{u_{0}^{n−1}} + (n − 1)kT = \frac{1 + (n − 1)u_{0}^{n−1}kT}{u_{0}^{n−1}}

u(T) = \frac{u_{0}}{(1 + (n − 1)u_{0}^{n−1}kT )^{1/(n−1)}}

The half-life T_{1/2} is given by

u(T_{1/2}) = \frac{u_{0}}{(1 + (n − 1)u_{0}^{n−1}kT_{1/2})^{1/(n−1)}} = \frac{1}{2}u_{0}

which implies

T_{1/2} = \frac{2^{n−1} − 1}{(n − 1)u_{0}^{n−1}k}

1.2.

(a) u(t) = u_{0}e^{kt}

(b–c) order n > 1: u(t) = \frac{u_{0}}{(1 − (n − 1)u_{0}^{n−1}kt)^{1/(n−1)}} , 0 ≤ t < \frac{1}{(n − 1)u_{0}^{n−1}k}

(d) u(t) = \frac{u_{0}}{u_{0} + (1 − u_{0})e^{−kt}}
(e) u(t) = \frac{u_{0}}{1 + t^{2}u_{0}/2}