DYNAMISKT SYSTEM. STABILITET

STIG LARSSON

Matematik, Chalmers tekniska högskola

Innehåll

1. Dynamiskt system. Linjärisering. Stabilitet 1 1.1. Dynamiskt system 1 1.2. Stationär lösning 2 1.3. Linjärisering 3 1.4. Homogent system 5 1.5. Matrisexponentialfunktionen 6 1.6. Symmetriskt system 6 1.7. Osymmetriskt system 7 1.8. Inhomogent system 8 1.9. Styva differentialekvationer 9 2. Övningar 10 3. Svar och lösningar till vissa övningar 11

1. Dynamiskt system. Linjärisering. Stabilitet

1.1. Dynamiskt system. Vi börjar med ett enkelt exempel. Exempel. Rörelseekvationerna för en dämpad, driven pendel är

(1)  mlθ′′(t) = −mg sin θ(t) − blθ′(t) + M(t)/l, [N=kgms^{−2}]

med beteckningar enlig Figur 1, där m [kg] är pendelns massa, l [m] dess längd, g [9.81 ms^{−2}] tyngdaccelerationen, θ(t) utslagsvinkeln vid tiden t [s], b [Ns/m] är en dämpningskoefficient (“luftmotståndet” är proportionellt mot farten) och M(t) [Nm] är ett vridmoment (“en elmotor på axeln”). Symbolerna inom klammer [ ] anger SI-enheter. Vi säger t ex att l bär dimensionen längd och mäts i SI-enheten m och att M(t) bär dimensionen kraft·längd och mäts i Nm. Vinkeln mäts i radianer och är dimensionslös. Det är till stor hjälp att göra sådan dimensionsanalys när man ställer upp matematiska modeller för tekniska system.

PSfrag replacements

M(t) m l mg θ

Figur 1. Pendeln.

Date: 25 November 2001 /stig.

Med beteckningarna X_{1} = θ, X_{2} = θ′, U = M, ω = \sqrt{g/l} [s^{−1}], α = b/m [s^{−1}], β = 1/(ml^{2}) får vi

(2)  \begin{matrix}X′_{1}(t) = X_{2}(t), & [s^{−1}] \\ X′_{2}(t) = −ω^{2} sin X_{1}(t) − αX_{2}(t) + βU(t). & [s^{−2}] \end{matrix}

Detta kan skrivas på matrisform,

X′(t) = F (X(t), U(t)), där X = [\begin{matrix}X_{1} \\ X_{2} \end{matrix}], F = [\begin{matrix}F_{1} \\ F_{2} \end{matrix}],

och F_{1}(X, U) = X_{2}, F_{2}(X, U) = −ω^{2} sin X_{1}(t) − αX_{2}(t) + βU(t). Till dessa ekvationer kommer ett begynnelsevillkor

X(0) = [\begin{matrix}θ_{0} \\ ν_{0} \end{matrix}] ,

där θ_{0} och ν_{0} är pendelns utslagsvinkel respektive vinkelhastighet vid tiden t = 0. □

Mer allmänt betraktar vi ett begynnelsevärdesproblem för ett system av icke-linjära ordinära differentialekvationer av typen

(5)  \begin{matrix}X′(t) = F (X(t), U(t)), t > 0, \\ X(0) = X_{0}, \end{matrix}

där X(t) ∈ R^{n} är den obekanta lösningen, F : R^{n} ×R^{m} → R^{n} är en given vektorvärd funktion av X ∈ R^{n} och U ∈ R^{m}, U(t) ∈ R^{m} är en given funktion och X_{0} ∈ R^{n} är ett givet begynnelsevärde. I exemplet är n = 2 och m = 1.

Vi betraktar detta som ett dynamiskt system, där vi ger indata X_{0} och U(t) och beräknar utdata X(t), se Figur 2. I detta sammanhang kallas X tillståndsvariabler och U styrvariabler. Tillståndsvariablerna kan t ex representera temperaturer och koncentrationer av reagerande ämnen i en kemisk reaktor och styrvariablerna kan representera flöden av tillförda ämnen och värmeflöden i värmeväxlare för kylning och uppvärmning. Ofta kan man inte mäta (observera) tillståndsvariablerna X direkt och inför då observerbara variabler Y ∈ R^{k}, som beror på X enligt Y (t) = G(X(t)), vilket representerar mätprocessen. Information från Y kan sedan adderas till U i en återkopplingsslinga (feedback loop) (t ex om det blir för varmt ökar man kylningen). Vi bortser från detta här och arbetar bara med X, U, enligt Figur 2. Observera att systemet är autonomt om det inte finns några styrvariabler: X′ = F (X).

Vi är inte intresserade av en enskild lösning X(t) utan av systemets beteende då vi varierar indata. Differentialekvationen är icke-linjär och kan i allmänhet endast lösas numeriskt. Viss information om systemets beteende kan man dock få genom att linjärisera, vilket vi nu skall beskriva.

PSfrag replacements X(t)

X_{0}
U(t)

Figur 2. Dynamiskt system.

1.2. Stationär lösning. Antag att systemet (5) har ett stationärt (tidsoberoende) tillstånd X̄, Ū. Eftersom X̄′ = 0 ger (5) att

(6)  F (X̄, Ū) = 0.

I tillämpningar kan X̄ representera en arbetspunkt (eller jämviktsläge), som man vill uppnå genom att ställa in styrvariabeln Ū på lämpligt sätt.

Exempel. För pendeln fås de stationära punkterna genom att sätta X′_{1}(t) = X′_{2}(t) = 0 i (2), vilket leder till ekvationssystemet

(7)  \begin{aligned}X̄_{2} & = 0, \\ −ω^{2} sin X̄_{1} − αX̄_{2} + βŪ & = 0. \end{aligned}

Om vi väljer t ex Ū = 0 (dvs inget yttre vridmoment) får vi lösningarna X̄_{2} = 0, X̄_{1} = kπ, k = 0, ±1, ±2, . . .. På grund av periodicitet räcker det att betrakta de två fallen X̄_{1} = 0 (nedre jämviktsläget) och X̄_{1} = π (övre jämviktsläget). Med Ū = 0 har vi alltså två stationära punkter

X̄ = [\begin{matrix}0 \\ 0 \end{matrix}] och X̄ = [\begin{matrix}π \\ 0 \end{matrix}].

Å andra sidan, om vi vill att pendeln ska ställa in sig på vinkeln X̄_{1} = ϕ (0 ≤ ϕ < π/2), väljer vi Ū så att

−ω^{2} sin ϕ + βŪ = 0,

dvs Ū = \frac{ω^{2}}{β} sin ϕ. Med detta värde på Ū ges de stationära lösningarna av X̄_{2} = 0 och ekvationen

−ω^{2} sin X̄_{1} + ω^{2} sin ϕ = 0,

med lösningarna

X̄_{1} = \begin{cases}ϕ + k2π \\ π − ϕ + k2π \end{cases} , k = 0, ±1, ±2, ….

Med Ū = \frac{ω^{2}}{β} sin ϕ får vi alltså (bortsett från periodiska upprepningar) två stationära punkter

X̄ = [\begin{matrix}ϕ \\ 0 \end{matrix}] och X̄ = [\begin{matrix}π − ϕ \\ 0 \end{matrix}]. □

Vi kan alltså ta två attityder: (i) med styrvariabeln Ū given, bestäm vilka jämviktstillstånd X̄ som är möjliga; (ii) bestäm styrvariabeln Ū så att ett givet jämviktstillstånd X̄ uppnås.

1.3. Linjärisering. Låt nu X(t) med data X_{0}, U(t) vara en lösning till (5) som är nära X̄. Med

(11)  X(t) = X̄ + ΔX(t), U(t) = Ū + ΔU(t), X_{0} = X̄ + ΔX_{0},

kan vi betrakta ΔX(t) som en störning i tillståndet X(t) orsakad av störningarna ΔX_{0} och ΔU(t) i data. Vi söker en ekvation för ΔX(t). Vi har

ΔX′(t) = X′(t) = F (X(t), U(t)) = F (X̄ + ΔX(t), Ū + ΔU(t)).

Definitionen av derivatan av F i punkten X̄, Ū ger

F (X̄ + ΔX(t), Ū + ΔU(t)) = F (X̄, Ū) + F′_{X}(X̄, Ū)ΔX(t) + F′_{U}(X̄, Ū)ΔU(t) + R(t),

där resttermen R(t) = E_{F} (X(t), U(t), X̄, Ū) uppfyller ∥R(t)∥ ≤ K_{F} (∥ΔX(t)∥^{2} + ∥ΔU(t)∥^{2}). Utskrivet i komponentform blir detta

F_{i}(X̄ + ΔX(t), Ū + ΔU(t)) = F_{i}(X̄, Ū) + [\frac{∂F_{i}}{∂X_{1}} (X̄, Ū), …, \frac{∂F_{i}}{∂X_{n}} (X̄, Ū)] [ΔX_{1}(t) ; ⋮ ; ΔX_{n}(t)]

+ [\frac{∂F_{i}}{∂U_{1}} (X̄, Ū), …, \frac{∂F_{i}}{∂U_{m}} (X̄, Ū)] [ ΔU_{1}(t) ; ⋮ ; ΔU_{m}(t)] + R_{i}(t)

= F_{i}(X̄, Ū) + ∑_{j=1}^{n} \frac{∂F_{i}}{∂X_{j}} (X̄, Ū)ΔX_{j}(t)

+ ∑_{j=1}^{m} \frac{∂F_{i}}{∂U_{j}} (X̄, Ū)ΔU_{j}(t) + R_{i}(t).

(14)

Eftersom F (X̄, Ū) = 0, får vi (på matrisform)

F (X̄ + ΔX(t), Ū + ΔU(t)) = AΔX(t) + BΔU(t) + R(t),

där matriserna

(16)  \begin{aligned}A & = F′_{X}(X̄, Ū) = [\begin{matrix}\frac{∂F_{1}}{∂X_{1}} (X̄, Ū) & … & \frac{∂F_{1}}{∂X_{n}} (X̄, Ū) \\ ⋮ & & ⋮ \\ \frac{∂F_{n}}{∂X_{1}} (X̄, Ū) & … & \frac{∂F_{n}}{∂X_{n}} (X̄, Ū) \end{matrix}] , \\ B & = F′_{U}(X̄, Ū) = [\begin{matrix}\frac{∂F_{1}}{∂U_{1}} (X̄, Ū) & … & \frac{∂F_{1}}{∂U_{m}} (X̄, Ū) \\ ⋮ & & ⋮ \\ \frac{∂F_{n}}{∂U_{1}} (X̄, Ū) & … & \frac{∂F_{n}}{∂U_{m}} (X̄, Ū) \end{matrix}] , \end{aligned}

kallas Jacobimatriser till F i punkten X̄, Ū. Den sökta ekvationen för ΔX(t) är alltså

(17)  \begin{matrix}ΔX′(t) = AΔX(t) + BΔU(t) + R(t), t > 0, \\ ΔX(0) = ΔX_{0}. \end{matrix}

Observera att denna endast är en omskrivning av (5). Om ΔX(t) och ΔU(t) är små, så kan vi förenkla (17) genom att försumma resttermen, och vi får det linjäriserade systemet

(18)  \begin{matrix}x′(t) = Ax(t) + Bu(t), t > 0, \\ x(0) = x_{0}, \end{matrix}

för den approximativa störningen x(t) ≈ ΔX(t) orsakad av datastörningarna x_{0} = ΔX_{0} och u(t) = ΔU(t). Notera sambanden X(t) ≈ X̄ + x(t), X_{0} = X̄ + x_{0}, U(t) = Ū + u(t), som gäller så länge som störningarna är små.

Det linjära systemet (18) sägs vara stabilt om lösningen x(t) förblir liten för alla t ≥ 0, om man tar tillräckligt små indata x_{0}, u(t). Man kan visa att i så fall förblir X(t) nära X̄, om man väljer X_{0} tillräckligt nära X̄ och U(t) tillräckligt nära Ū. Vi säger då att X̄, Ū är en stabil stationär punkt till det icke-linjära systemet (5). Observera att triviala lösningen x(t) = 0, u(t) = 0 till (18) motsvarar stationära lösningen X(t) = X̄, U(t) = Ū till (5). Exempel. För pendeln blir Jacobimatriserna

(19) A = [\begin{matrix}\frac{∂F_{1}}{∂X_{1}} (X̄, Ū) & \frac{∂F_{1}}{∂X_{2}} (X̄, Ū) \\ \frac{∂F_{2}}{∂X_{1}} (X̄, Ū) & \frac{∂F_{2}}{∂X_{2}} (X̄, Ū) \end{matrix}] = [\begin{matrix}0 & 1 \\ −ω^{2} cos X̄_{1} & −α \end{matrix}] , B = [\begin{matrix}\frac{∂F_{1}}{∂U} (X̄, Ū) \\ \frac{∂F_{2}}{∂U} (X̄, Ū) \end{matrix}] = [\begin{matrix}0 \\ β \end{matrix}] .

I punkten X̄ = [\begin{matrix}0 \\ 0 \end{matrix}], Ū = 0 får vi det linjäriserade systemet

(20)  x′(t) = [\begin{matrix}0 & 1 \\ −ω^{2} & −α \end{matrix}] x(t) + [\begin{matrix}0 \\ β \end{matrix}] u(t),

vilket motsvarar den linjära, dämpade, drivna svängningsekvationen

(21)  θ′′(t) + ω^{2}θ(t) + αθ′(t) = βu(t).

I punkten X̄ = [\begin{matrix}π \\ 0 \end{matrix}], Ū = 0 får vi det linjäriserade systemet

(22)  x′(t) = [\begin{matrix}0 & 1 \\ ω^{2} & −α \end{matrix}] x(t) + [\begin{matrix}0 \\ β \end{matrix}] u(t),

vilket motsvarar

(23)  θ′′(t) − ω^{2}θ(t) + αθ′(t) = βu(t),

vars lösningar inte är svängningar. Se övning 7.

Man kan enkelt göra en numerisk lösning med Matlab. Gör detta! Följande m-fil pendel.m definierar högerledet i ekvation (20) med β = 0.

  function y=pendel(t,x)
  omega=5;  alpha=1;
  A=[0 1; -omega^2 -alpha];
  y=A*x;

Följande kommandon löser ekvationen och plottar sedan lösningskurvorna, dvs x_{1}(t) och x_{2}(t) som funktion av t.

  T=10; x0=[1;1];
  [t,x]=ode23('pendel', [0 T], x0]);
  plot(t,x);

Men du använder givetvis din egen ode-lösare my ode istället för ode23!

Kommandot plot(x(:,1),x(:,2)) plottar x_{2} mot x_{1}. Detta kallas fasporträtt, se Figur 3. □

PSfrag replacements X(t)

X_{0}

U(t)

Figur 3. Den linjäriserade pendeln. Lösningskurvor och fasporträtt.

Vi ska nu undersöka stabiliteten hos det linjäriserade systemet (18). Vi börjar med den homogena ekvationen, dvs fallet då Bu(t) ≡ 0.

1.4. Homogent system. Vi betraktar ett system av linjära differentialekvationer med konstanta koefficienter

(24)  x′(t) = Ax(t), t > 0; x(0) = b.

Antag att koefficientmatrisen A är diagonaliserbar, dvs att egenvektormatrisen P = [g_{1}, …, g_{n}] är inverterbar, så att P^{−1}AP = D, D = diag{λ_{1}, …, λ_{n}}. Då utgör egenvektorerna {g_{k}}_{k=1}^{n} en bas för R^{n} (eller C^{n} om g_{k} är komplexa). Varje vektor kan då skrivas på ett entydigt sätt som en linjärkombination av {g_{k}}_{k=1}^{n}. Speciellt gör vi ansatsen x(t) = ∑_{k=1}^{n} y_{k}(t)g_{k}. Koefficienterna y_{k}(t) bestäms genom insättning i differentialekvationen. Om vi använder relationerna Ag_{k} = λ_{k}g_{k} får vi

0 = x′(t) − Ax(t) = ∑_{k=1}^{n} y′_{k}(t)g_{k} − ∑_{k=1}^{n} y_{k}(t)Ag_{k} = ∑_{k=1}^{n} (y′_{k}(t) − λ_{k}y_{k}(t))g_{k}.

Identifikation av koefficienterna för g_{k} ger ett okopplat system

(26)  y′_{k}(t) = λ_{k}y_{k}(t), k = 1, …, n,

med lösningen y_{k}(t) = c_{k}e^{λ_{k}t}, vilket ger

(27)  x(t) = ∑_{k=1}^{n} c_{k}e^{λ_{k}t}g_{k} = c_{1}e^{λ_{1}t}g_{1} + ⋯ + c_{n}e^{λ_{n}t}g_{n}.

Koefficienterna c_{k} ges av begynnelsevillkoret x(0) = b,

b = ∑_{k=1}^{n} c_{k}g_{k},

dvs c_{k} är komponenterna för vektorn b med avseende på basen {g_{k}}_{k=1}^{n}.

Likheten (27) framställer lösningen x(t) som en linjärkombination av egenmoder e^{λ_{k}t}g_{k}. Karaktären hos dessa bestäms av egenvärdena λ_{k}. Om vi skriver λ_{k} = α_{k} + iω_{k}, där α_{k} och ω_{k} ar real- respektive imaginärdelarna av λ_{k}, kan vi urskilja två fall:

• λ_{k} = α_{k} är reellt. Då är e^{λ_{k}t} = e^{α_{k}t} exponentiellt växande, exponentiellt avtagande eller konstant beroende på om α_{k} > 0, α_{k} < 0 eller α_{k} = 0.

• λ_{k} = α_{k} + iω_{k} är inte reellt (ω_{k} ≠ 0). Då är

e^{λ_{k}t} = e^{α_{k}t+iω_{k}t} = e^{α_{k}t}e^{iω_{k}t} = e^{α_{k}t}(cos ω_{k}t + i sin ω_{k}t)

en harmonisk svängning med amplituden |e^{λ_{k}t}| = e^{α_{k}t} och vinkelfrekvensen ω_{k}. Amplituden är exponentiellt växande, exponentiellt avtagande eller konstant beroende på om α_{k} > 0, α_{k} < 0 eller α_{k} = 0.

Vi kan nu dra vissa slutsatser om lösningens storlek. Vi säger att det linjäriserade systemet (24) är

• stabilt om det gäller att om b är liten, så är x(t) liten för alla t ≥ 0,

• instabilt annars. Dessutom säger vi att det linjäriserade systemet (24) är

• asymptotiskt stabilt om x(t) → 0 då t → ∞ för alla begynnelsevärden b. Kom ihåg att lösningen till det linjäriserade systemet (24) spelar rollen av approximativ störning

till det icke-linjära systemet (5), x(t) ≈ ΔX(t), X(t) ≈ X̄ + x(t). Observera att ett asymptotiskt stabilt system kan vara instabilt: störningen x(t) skulle kunna växa sig stor i början innan den till slut går mot noll, se avsnitt 1.7.

Om A är diagonaliserbar så kan vi dra vissa slutsatser genom att undersöka stabiliteten för varje egenmod separat. Av resonemanget ovan ser vi att stabiliteten för en egenmod e^{λ_{k}t}g_{k} beror på tecknet på realdelen Re λ_{k} = α_{k}. Man kan visa att om A ej är diagonaliserbar (dvs om A har färre än n linjärt oberoende egenvektorer), så är egenmoderna av typen p_{k}(t)e^{λ_{k}t}g_{k}, där p_{k} är ett polynom. Stabiliteten av en sådan mod bestäms också av tecknet på realdelen av λ_{k}. Vi formulerar detta som en sats. Sats 1. Det linjära systemet x′ = Ax är asymptotiskt stabilt om och endast om Re λ_{k} < 0 för alla egenvärden λ_{k} till A, och instabilt om Re λ_{k} > 0 för något egenvärde λ_{k}.

Däremot kan man inte alltid dra slutsatsen att systemet är stabilt om varje egenmod är stabil, vilket vi strax ska se i avsnitt 1.7.

1.5. Matrisexponentialfunktionen. Om matrisen A är diagonaliserbar, dvs om A = P DP^{−1}, så definierar vi matrisen

e^{tA} = P e^{tD}P^{−1}, där e^{tD} = diag{e^{tλ_{1}}, …, e^{tλ_{n}}}.

Genom transformationen x = P y, b = P c, övergår ekvation (24) i

y′(t) = Dy(t), t > 0; y(0) = c,

vilken är identisk med (26), med lösningen y(t) = e^{tD}c. Genom att transformera tillbaka får vi

x(t) = P y(t) = P e^{tD}c = P e^{tD}P^{−1}b = e^{tA}b,

dvs

(32)  x(t) = e^{tA}b,

vilket är ett annat sätt att skriva (27). Matrisen e^{tA} kan ses som en lösningsoperator till det homogena systemet (24). Observera att

(33)  \frac{d}{dt}e^{tA} = Ae^{tA} = e^{tA}A.

1.6. Symmetriskt system. Sats 1 är ett kvalitativt resultat: i det asymptotiskt stabila fallet får vi t ex veta att störningarna går mot noll då t → ∞, men inte hur stora de kan bli innan de dör ut. Om A är symmetrisk så kan vi även få kvantitativ information. Denna kan uttryckas med hjälp av stabilitetsfaktorn

(34)  S(T) = max_{b} \frac{max_{0≤t≤T}∥x(t)∥}{∥b∥} ,

som anger hur mycket störningarna kan växa (i värsta fall) på intervallet [0, T].

Antag alltså att A är symmetrisk. Enligt spektralsatsen är λ_{k} reella och vi kan välja g_{k} ortonormerade. Vi har då

\begin{matrix}b = ∑_{k=1}^{n} c_{k}g_{k}, c_{k} = (b, g_{k}) = b^{t}g_{k}, ∥b∥^{2} = ∑_{k=1}^{n} c_{k}^{2}, \\ x(t) = ∑_{k=1}^{n} c_{k}e^{λ_{k}t}g_{k}, ∥x(t)∥^{2} = ∑_{k=1}^{n}(c_{k}e^{λ_{k}t})^{2} = ∑_{k=1}^{n} c_{k}^{2}e^{2λ_{k}t}. \end{matrix}

Vi ordnar egenvärdena λ_{n} ≤ … ≤ λ_{k} ≤ … ≤ λ_{1}, så att e^{λ_{k}t} ≤ e^{λ_{1}t} och

(36)  ∥x(t)∥^{2} = ∑_{k=1}^{n} c_{k}^{2}e^{2λ_{k}t} ≤ e^{2λ_{1}t} ∑_{k=1}^{n} c_{k}^{2} = e^{2λ_{1}t}∥b∥^{2},

dvs

∥x(t)∥ ≤ e^{λ_{1}t}∥b∥, t ≥ 0,

med likhet om b = g_{1}. Vi ser att en störning x(t) inte kan bli större än begynnelsestörningen x(0) = b, om λ_{1} ≤ 0; stabilt system. Om λ_{1} < 0, så avtar alla störningar exponentiellt, då t → ∞; asymptotiskt stabilt system. Om λ_{1} > 0, så finns det störningar som växer exponentiellt; instabilt system.

Stabilitetsfaktorn blir

S(T) = max_{b} \frac{max_{0≤t≤T}∥x(t)∥}{∥b∥} ≤ \frac{max_{0≤t≤T}(e^{λ_{1}t}∥b∥)}{∥b∥} = max_{0≤t≤T} e^{λ_{1}t} = \begin{cases}e^{λ_{1}T} , & om λ_{1} > 0, \\ 1, & om λ_{1} ≤ 0. \end{cases}

Eftersom vi har likhet då b = g_{1}, drar vi slutsatsen att

S(T) = \begin{cases}e^{λ_{1}T} , & om λ_{1} > 0, \\ 1, & om λ_{1} ≤ 0. \end{cases}

Sats 2. Om matrisen A är symmetrisk med (de reella) egenvärdena λ_{k}, så är systemet x′ = Ax stabilt om och endast om alla λ_{k} ≤ 0, asymptotiskt stabilt om och endast om alla λ_{k} < 0, och instabilt om och endast om något λ_{k} > 0.

1.7. Osymmetriskt system. För en symmetrisk matris är egenvektorerna ortogonala och det är anledningen till att vi kunde göra den uttömmande stabilitetsundersökningen i det fallet. En osymmetrisk matris kan ha nästan parallella egenvektorer eller till och med parallella, dvs sammanfallande, egenvektorer. Vi nöjer oss med att undersöka följande typiska fall:

x′ = Ax, A = [\begin{matrix}a & 1 \\ ϵ^{2} & a \end{matrix}] ,

där a är ett reellt eller komplext tal och ϵ är ett litet positivt tal. Här blir egenvärdes- och egenvektormatriserna

D = [\begin{matrix}a + ϵ & 0 \\ 0 & a − ϵ \end{matrix}] , P = [\begin{matrix}1 & 1 \\ ϵ & −ϵ \end{matrix}] , P^{−1} = \frac{1}{2} [\begin{matrix}1 & ϵ^{−1} \\ 1 & −ϵ^{−1} \end{matrix}] .

Vi ser att egenvärdena nästan sammanfaller och att egenvektorerna är nästan parallella. Då ϵ = 0 är a ett dubbelt egenvärde med bara en egenvektor och matrisen A är inte diagonaliserbar.

En enkel räkning, se (32), visar att

x(t) = e^{tA}b = P e^{tD}P^{−1}b

= \frac{1}{2} [\begin{matrix}e^{(a+ϵ)t} + e^{(a−ϵ)t} & ϵ^{−1}(e^{(a+ϵ)t} − e^{(a−ϵ)t}) \\ ϵ(e^{(a+ϵ)t} − e^{(a−ϵ)t}) & e^{(a+ϵ)t} + e^{(a−ϵ)t} \end{matrix}] b

= e^{at} [\begin{matrix}cosh(ϵt) & ϵ^{−1} sinh(ϵt) \\ ϵ sinh(ϵt) & cosh(ϵt) \end{matrix}] b,

(42) där cosh(ϵt) = \frac{1}{2}(e^{ϵt} + e^{−ϵt}), sinh(ϵt) = \frac{1}{2}(e^{ϵt} − e^{−ϵt}). I Figur 4 ses resultatet av en numerisk beräkning gjord med Matlab. (Gör denna beräkning själv!) Systemet är asymptotiskt stabilt enligt Sats 1, eftersom båda egenvärdena är negativa. Vi noterar att komponenten x_{2}(t) avtar exponentiellt, medan x_{1}(t) först växer kraftigt innan det exponentiella avtagandet sätter in. Systemet är asymptotiskt stabilt, eftersom alla störningar går mot 0, men det kan knappast kallas stabilt, eftersom störningen kan bli mycket större än begynnelsestörningen.

PSfrag replacements X(t)

X_{0}
U(t)
Figur 4. ϵ = 0.01, a = −0.1, x(0) = [0, 1]^{t}.

Vi har nu sett att ett osymmetriskt system kan vara asymptotiskt stabilt men ändå instabilt. Man kan visa att om matrisens egenvektorer inte är “nästan parallella eller parallella”, så kan man dra slutsats om stabilitet genom att undersöka varje egenmod för sig. Annars kan man bara dra slutsats om asymptotisk stabilitet och om instabilitet.

1.8. Inhomogent system. Om Bu(t) ≢ 0, så ges lösningen till (18) av

(43)  x(t) = e^{tA}x_{0} + ∫_{0}^{t} e^{(t−s)A}Bu(s) ds,

där matrisen e^{tA} är lösningsoperatorn till det homogena systemet (24), se övning 1. Om nu det homogena systemet är stabilt, dvs om operatorn e^{tA} är stabil, så kan man visa att termen ∫_{0}^{t} e^{(t−s)A}Bu(s) ds också är stabil, dvs den blir liten om u är liten. Exempel. Vi återvänder nu till pendeln. Jacobimatrisen för det övre jämviktsläget är, se (22),

A = [\begin{matrix}0 & 1 \\ ω^{2} & −α \end{matrix}] ,

med egenvärdena

λ_{±} = −\frac{α}{2} ± \sqrt{\frac{α^{2}}{4} + ω^{2}}.

Eftersom λ_{+} > 0 och λ_{−} < 0 drar vi slutsatsen att det övre jämviktsläget är instabilt. I det undre jämviktsläget, se (20), blir Jacobimatrisen

A = [\begin{matrix}0 & 1 \\ −ω^{2} & −α \end{matrix}] ,

med egenvärdena

λ_{±} = −\frac{α}{2} ± \sqrt{\frac{α^{2}}{4} − ω^{2}}.

Om α > 2ω (stark dämpning) får vi två reella egenvärden λ_{±} < 0. Om α = 2ω (kritisk dämpning) får vi ett dubbelt reellt egenvärde λ_{±} = −α/2 < 0 (Jacobimatrisen är ej diagonaliserbar i detta fall). Om 0 < α < 2ω (svag dämpning) får vi två komplexa egenvärden med Re λ_{±} = −α/2 < 0. Om α = 0 (ingen dämpning) får vi två imaginära egenvärden λ_{±} = ±iω med Re λ_{±} = 0. Vi drar slutsatsen att det nedre jämviktsläget är asymptotiskt stabilt för α > 0. Eftersom Re λ_{k} ≤ 0 för α ≥ 0, misstänker vi att det nedre jämviktsläget är stabilt. Men vi vet att man måste vara försiktig med den slutsatsen då systemet är osymmetriskt. Ytterligare undersökning, t ex numeriska experiment, visar att systemet faktiskt är stabilt för α ≥ 0. □

1.9. Styva differentialekvationer. Antag för enkelhets skull att A är symmetrisk så att den har reella egenvärden λ_{k} och en ON-bas g_{k} av egenvektorer. Vi antar dessutom att alla egenvärdena är negativa, så att det linjära systemet

(48)  x′(t) = Ax(t), 0 < t < T;  x(0) = b,

är stabilt och asymptotiskt stabilt, se Sats 2. Med numreringen λ_{n} ≤ ⋯ ≤ λ_{1} < 0 har vi

x(t) = c_{1}e^{λ_{1}t}g_{1} + ⋯ + c_{n}e^{λ_{n}t}g_{n},

där den sista egenmoden alltså avtar snabbare än den första. Kvoten

(49)  κ(A) = \frac{max |λ_{k}|}{min |λ_{k}|}

kallas konditionstalet för den symmetriska matrisen A. I vårt fall blir κ(A) = λ_{n}/λ_{1}. Notera att κ(A) ≥ 1. Systemet sägs vara styvt om κ(A) >> 1, dvs om lösningarna innehåller både snabba och långsamma egenmoder. (För osymmetriska matriser definieras konditionstalet på ett lite mera komplicerat sätt.)

Den explicita Eulers metod med konstant steglängd h för approximativ lösning av (48) är

\frac{y_{j} − y_{j−1}}{h} = Ay_{j−1}; y_{0} = b,

dvs vi har rekursionsekvationen

y_{j} = (I + hA)y_{j−1}; y_{0} = b.

Vi får

y_{j} = (I + hA)^{j}b.

Matrisen B = I + hA är symmetrisk med egenvärdena μ_{k} = 1 + hλ_{k} och samma egenvektorer g_{k} som A. Detta ger att

y_{j} = ∑_{k=1}^{n} c_{k}μ_{k}^{j}g_{k}; där b = ∑_{k=1}^{n} c_{k}g_{k},

och vi ser att Eulers metod är stabil om och endast om alla egenvärdena har absolutbelopp mindre än eller lika med 1, dvs |μ_{k}| = |1 + hλ_{k}| ≤ 1, dvs −1 ≤ 1 + hλ_{k} ≤ 1, dvs −2 ≤ hλ_{k} ≤ 0 för alla λ_{k}. Eftersom λ_{n} ≤ … ≤ λ_{k} ≤ … ≤ λ_{1} < 0 är detta uppfyllt om och endast om

(50)  h ≤ \frac{2}{|λ_{n}|} = \frac{2}{max |λ_{k}|},

dvs steglängden måste vara liten om max |λ_{k}| är stort. Det intressanta tidsintervallet beror å andra sidan på den långsammaste egenmoden: T ≥ 1/|λ_{1}| = 1/min |λ_{k}|. Antalet tidssteg blir således

N = \frac{T}{h} ≥ \frac{1}{2} \frac{max |λ_{k}|}{min |λ_{k}|} = \frac{1}{2}κ(A).

Om systemet är styvt måste vi ta många steg för att få med hela förloppet: små steg, långt tidsintervall.

Den implicita Eulers metod är

\frac{y_{j} − y_{j−1}}{h} = Ay_{j}; y_{0} = b,

dvs

(I − hA)y_{j} = y_{j−1}; y_{0} = b.

Metoden kallas implicit eftersom detta är ett ekvationssystem som måste lösas i varje tidssteg. När detta är gjort har vi en rekursion

y_{j} = (I − hA)^{−1}y_{j−1}; y_{0} = b.

Matrisen B = (I − hA)^{−1} har egenvärdena μ_{k} = (1 − hλ_{k})^{−1} och metoden är stabil om och endast om alla egenvärdena har absolutbelopp mindre än eller lika med 1, dvs

−1 ≤ \frac{1}{1 − hλ_{k}} ≤ 1 för alla λ_{k}.

Eftersom λ_{k} ≤ 0 är detta villkor uppfyllt för alla steglängder h. Den implicita Eulers metod är alltså lämpligare för styva problem än den explicita metoden. Man kan välja steglängden enbart efter den önskade noggrannheten utan att ta extra hänsyn till ett stabilitetsvillkor som i (50). I vårt exempel (48) kan man t ex ta små steg i början tills de snabba moderna dött ut, sedan stora steg under det långa intervall då de långsamma moderna pågår. Med den explicita metoden måste man ta små steg hela tiden för att gardera sig mot instabilitet.

Styva system av differentialekvationer uppträder ofta i kemisk reaktionsteknik där olika delreaktioner i en sammansatt reaktion kan ha mycket olika reaktionshastigheter. Styva system uppträder också i samband med numerisk lösning av partiella differentialekvationer, t ex värmeledningsekvationen. Sammanfattningsvis kan sägas att man behöver ta hänsyn till styvheten när man väljer numerisk metod för differentialekvationsproblem. Lämpliga metoder är ofta implicita. Matlab har programmen ode23 och ode45 för icke-styva system och ode15s för styva system. Nu kan du göra övning 8.

2. Övningar

Svar och lösningar till vissa övningar finns i nästa avsnitt. Övning 1. Bevisa (43) genom att multiplicera ekvation (18) med den integrerande faktorn e^{−tA} (se (33)) och integrera. Övning 2. Avgör om systemet x′ = Ax är stabilt när A är

(a) [\begin{matrix}−1001 & 999 \\ 999 & −1001 \end{matrix}] (b) [\begin{matrix}−1 & 3 & 5 \\ 0 & −4 & 6 \\ 0 & 0 & 1 \end{matrix}] (c) [\begin{matrix}0 & 1 \\ −1 & −4 \end{matrix}]

Övning 3. Bestäm de stationära lösningarna till systemet

X′_{1}(t) = X_{2}(t)(1 − X_{1}(t)^{2}),
X′_{2}(t) = 2 − X_{1}(t)X_{2}(t).

Linjärisera kring de stationära lösningarna och undersök de linjäriserade systemens stabilitet. Övning 4. Mineas ekvation. Bestäm vilka stationära lösningar man kan få för olika värden på Ū och δ > 0. Skriv för varje stationär lösning ned den linjäriserade ekvationen och undersök dess stabilitet.

X′_{1}(t) = −X_{1}(t) − δ(X_{2}(t)^{2} + X_{3}(t)^{2}) + U(t),
X′_{2}(t) = −X_{2}(t) − δX_{1}(t)X_{2}(t),
X′_{3}(t) = −X_{3}(t) − δX_{1}(t)X_{3}(t).

Övning 5. Bestäm de stationära punkterna och avgör deras stabilitet för följande system:

\begin{aligned}(a) X′_{1} & = X_{1}(1 − X_{2}) \\ X′_{2} & = X_{2}(1 − X_{1}) \end{aligned}  \begin{aligned}(b) X′_{1} & = −2(X_{1} − 10) + X_{2}e^{X_{1}} \\ X′_{2} & = −2X_{2} − X_{2}e^{X_{1}} \end{aligned}

\begin{aligned}(c) X′_{1} & = X_{1} + X_{1}X_{2}^{2} + X_{1}X_{3}^{2} \\ X′_{2} & = −X_{1} + X_{2} − X_{2}X_{3} + X_{1}X_{2}X_{3} \\ X′_{3} & = X_{2} + X_{3} − X_{1}^{2} \end{aligned}

Övning 6. Skriv följande begynnelsevärdesproblem som ett system x′ = Ax, x(0) = b genom att införa nya variabler x_{1} = u, x_{2} = u′. Lös systemet med diagonaliseringsmetoden.

u′′ + 4u′ + 13u = 0, (dämpad harmonisk svängning)

u(0) = 0, u′(0) = −1

Övning 7. Lös ekvationerna (21) och (23) med α = 0 och β = 0. Övning 8. Explicit/implicit metod, styva differentialekvationer. Skriv två Matlab-program som löser begynnelsevärdesproblemet x′ = Ax för t > 0, x(0) = x_{0} med explicit respektive implicit Eulers metod och konstant steg. Provkör med matriserna

A = [\begin{matrix}−1000 & −999 \\ −999 & −999 \end{matrix}] , A = [\begin{matrix}−1000 & −999 \\ −999 & −998 \end{matrix}] .

Prova med olika steglängder från h = 0.2 och mindre. Välj sluttid så att hela förloppet kommer med. Hur stora steg använder Matlabs ode23 och ode15s? Ser du någon skillnad mellan ode23 och ode15s när det gäller att lösa styva problem?

Tips: Matlab-programmet cond(A) beräknar konditionstalet och diff(t) kan användas för att beräkna stegen när man kört ode23 eller ode15s.

3. Svar och lösningar till vissa övningar

2. (a) asymptotiskt stabilt, (b) instabilt, (c) asymptotiskt stabilt

3. De stationära lösningarna ges av ekvationssystemet

X̄_{2}(1 − X̄_{1}^{2}) = 0
2 − X̄_{1}X̄_{2} = 0

med två lösningar (X̄_{1}, X̄_{2}) = (1, 2) och (X̄_{1}, X̄_{2}) = (−1, −2). Jacobimatrisen är

[\begin{matrix}−2X̄_{1}X̄_{2} & 1 − X̄_{1}^{2} \\ −X̄_{2} & −X̄_{1} \end{matrix}] .

Linjärisering kring (X̄_{1}, X̄_{2}) = (1, 2) leder till det linjära systemet

[\begin{matrix}x′_{1}(t) \\ x′_{2}(t) \end{matrix}] = [\begin{matrix}−4 & 0 \\ −2 & −1 \end{matrix}] [\begin{matrix}x_{1}(t) \\ x_{2}(t) \end{matrix}] ,

vars matris har negativa egenvärden −4 och −1, och alltså är asymptotiskt stabilt. Linjärisering kring (X̄_{1}, X̄_{2}) = (−1, −2) leder till systemet

[\begin{matrix}x′_{1}(t) \\ x′_{2}(t) \end{matrix}] = [\begin{matrix}−4 & 0 \\ 2 & 1 \end{matrix}] [\begin{matrix}x_{1}(t) \\ x_{2}(t) \end{matrix}] .

Egenvärdena är −4 och 1 och systemet är instabilt.

5. (a) (0, 0) instabil, (1, 1) instabil; (b) (10, 0) asymptotiskt stabil; (c) (0, 0, 0) och (0, −1, 1) båda instabila 6.

x(t) = c_{1}e^{λ_{1}t}g_{1} + c_{2}e^{λ_{2}t}g_{2} = c_{1}e^{(−2+3i)t} [\begin{matrix}1 \\ −2 + 3i \end{matrix}] + c_{2}e^{(−2−3i)t} [\begin{matrix}1 \\ −2 − 3i \end{matrix}] .

Begynnelsevillkoret x(0) = [\begin{matrix}0 \\ −1 \end{matrix}] ger c_{1} = −c_{2} = i/6, så att

x(t) = [\begin{matrix}−\frac{1}{3}e^{−2t} sin 3t \\ \frac{2}{3}e^{−2t} sin 3t − e^{−2t} cos 3t \end{matrix}] .

7. θ(t) = Ae^{iωt} + Be^{−iωt} = C cos ωt + D sin ωt, respektive θ(t) = Ae^{ωt} + Be^{−ωt} = C cosh ωt + D sinh ωt