TMV035 Analysis and Linear Algebra B^{1}

Ordinära differentialekvationer 3

1.1 Allmänt system av ODE

Vi har nu löst differentialekvationer

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

i några specialfall där f bara beror av x (myprim.m):

u′(x) = f(x)

och bara beror av u (myexp.m):

u′(x) = u(x)

och även system av denna senare typ (mytrig.m):

u′_{1}(x) = u_{2}(x),
u′_{2}(x) = −u_{1}(x).

Vi ska nu behandla differentialekvationer där f(x, u) är en allmän funktion av både x och u. Till exempel,

u′(x) = −xu(x), x ∈ [0, 3]; u(0) = 1,

med analytisk lösning u(x) = exp(−x^{2}/2). Verifiera detta! Här ges alltså f : R^{2} → R av formeln f(x, s) = −xs. När vi sätter in i högerledet av differentialekvationen ersätter vi den andra variabeln s med s = u(x). Vi ska även studera system av sådana ekvationer.

Vi kommer då till ett allmänt system av n stycken differentialekvationer:

u′_{1}(t) = f_{1}(t, u_{1}(t), …, u_{n}(t)), t ∈ [a, b],
⋮
u′_{n}(t) = f_{n}(t, u_{1}(t), …, u_{n}(t)), t ∈ [a, b],

tillsammans med lika många begynnelsevillkor:

u_{1}(a) = u_{a1}
⋮
u_{n}(a) = u_{an}.

Vi skriver dessa på matrisform:

u(t) = [u_{1}(t) ; ⋮ ; u_{n}(t)] , u_{a} = [u_{a1} ; ⋮ ; u_{an}] , f(t, u) = [f_{1}(t, u_{1}, …, u_{n}) ; ⋮ ; f_{n}(t, u_{1}, …, u_{n})] ,

och får följande begynnelsevärdesproblem

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

12006-11-21 /stig

För att lösa detta använder vi samma algoritm som tidigare (Eulers metod). Vi börjar med att dela in intervallet [a, b] i N stycken delintervall av längden h = (b − a)/N:

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

Vi beräknar nu en approximativ lösning enligt

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

Genom att förbinda punkterna (t_{i}, U(t_{i})) med räta linjer får vi en graf och funktionen U(t) blir definierad också mellan beräkningsnoderna t_{i}.

Implementering i Matlab

Algoritmen är

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

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

Övning 1. Skriv ett program myode.m med anropet [t,U]=myode(f,I,ua,h) som löser begynnelsevärdesproblemet (1). Du skall använda programskalet myode.m. Programmet blir en enkel modifikation av ditt tidigare program mytrig.m.

Prova programmet på följande begynnelsevärdesproblem. För varje exempel måste du skriva en funktionsfil av typen

function y=funk(t,u)

Plotta den numeriska lösningen tillsammans med den analytiska lösningen.

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

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

(c) \{ \begin{matrix}[\begin{matrix}u′_{1}(t) \\ u′_{2}(t) \end{matrix}] = [\begin{matrix}0 & 1 \\ −1 & 0 \end{matrix}] [\begin{matrix}u_{1}(t) \\ u_{2}(t) \end{matrix}] , t ∈ [0, 10], \\ u(0) = [\begin{matrix}0 \\ 1 \end{matrix}] . \end{matrix}

(d) \{ \begin{matrix}u′(x) = −xu(x), x ∈ [0, 3], \\ u(0) = 1. \end{matrix}

analytisk lösning: u(x) = exp(−x^{2}/2)

Facit

function y=funka(t,u)
y=t^2;
>> [x,U]=myode(@funka,[1,3],1,1e-3)
function y=funkb(t,u)
y=u;
>> [t,U]=myode(@funkb,[0,2],1,1e-3)
function y=funkc(t,u)
A=[0 1;-1 0];
y=A*u;
>> [t,U]=myode(@funkc,[0,3],[0;1],1e-3)
function y=funkd(t,u)
y=-t*u;
>> [x,U]=myode(@funkd,[0,3],1,1e-3)

1.2 Populationsdynamik

Vi betraktar population av bytesdjur (kaniner) som lever tillsammans med en population rovdjur (rävar). Låt u_{1}(t) respektive u_{2}(t) beteckna antalet kaniner respektive rävar vid tiden t. En enkel matematisk modell för populationernas utveckling ges av Volterra-Lotka-ekvationerna:

\begin{aligned}u′_{1}(t) & = au_{1}(t) − bu_{1}(t)u_{2}(t) \\ u′_{2}(t) & = −cu_{2}(t) + du_{1}(t)u_{2}(t) \end{aligned}  (2)

Koefficienterna a, b, c, d är positiva. Termen au_{1}(t) representerar netto-födelse-dödstalet i en ensam kaninpopulation. Termen −cu_{2}(t) är motsvarande för rävarna. Termen −bu_{1}(t)u_{2} är antalet kaniner som blir uppätna per tidsenhet. Termen du_{1}(t)u_{2}(t) är antalet rävar per tidsenhet som överlever på grund av tillgång på föda.

Observera teckenkombinationen i ekvationerna. Vad blir lösningen om populationerna är ensamma (b = d = 0)?

Lös Volterra-Lotka-ekvationerna med ditt program myode.m. Lämpliga värden är a = .5, b = 1, c = .2, d = 1.

Facit

Med b = d = 0 får vi u_{1}(t) = u_{01} exp(at), u_{2}(t) = u_{02} exp(−ct), dvs kaninpopulationen exploderar och rävarna dör ut. Matlab funktionsfil:

   function y=volterra(t,u)
   a=.5; b=1; c=.2; d=1;
   y=zeros(2,1);
   y(1)= a*u(1)-b*u(1)*u(2);
   y(2)=-c*u(2)+d*u(1)*u(2);

Matlab kommandon:

>> [t,U]=myode(@volterra,[0 50],[.5;.3],1e-2);
>> figure(1)
>> plot(t,U)
>> figure(2)
>> plot(U(:,1),U(:,2))