TMA855 MATEMATIK KX 1999. KOMPLETTERANDE FÖRELÄSNINGS-ANTECKNINGAR OCH ÖVNINGAR.

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 6 1.8. Inhomogent system 7 1.9. Styva differentialekvationer 8 2. Värmeledning, diffusion, reaktion 9 2.1. Härledning av värmeledningsekvationen 9 2.2. Diffusionsekvationen 12 2.3. Reaktions-diffusionsekvationer 12 3. Övningar 12 4. Svar och lösningar till vissa övningar 17 5. Tankreaktorn 19 5.1. Lösning 19

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

mlθ′′(t) = −mg sin θ(t) − blθ′(t) + M(t)/l, (1)

med beteckningar enlig Figur 1, och där b [Ns/m] är en dämpningskoefficient (“luftmotståndet”) och M(t) [Nm] är ett vridmoment (“en elmotor på axeln”).

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

\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} (2)

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}],

Date: September 1999.

PSfrag replacements

M(t) m l mg θ

Figur 1. Pendeln.

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

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

där X(t) ∈ R^{n} är den obekanta lösningen, F ∈ 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). I detta sammanhang kallas X tillståndsvariabel och U styrvariabel. 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 (3) har ett stationärt (tidsoberoende) tillstånd X̄, Ū. Eftersom X̄′ = 0 ger (3) att

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

I tillämpningar kan X̄ representera en arbetspunkt (eller jämviktsläge), som man vill uppnå genom att välja 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

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

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: (1) med styrvariabeln Ū given, bestäm vilka jämviktstillstånd X̄ som är möjliga; (2) 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 (3) som är nära X̄. Med

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

kan vi betrakta ΔX(t) som en störning i 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)).

Taylorutveckling (MA3 2.9) av F_{i} kring X̄, Ū ger

\begin{aligned}F_{i}(X̄ + ΔX(t), Ū + ΔU(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), \end{aligned} (7)

där R_{i}(t) är resttermen. 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

A = \{\frac{∂F_{i}}{∂X_{j}} (X̄, Ū)\}_{i=1, j=1}^{n, n} , B = \{\frac{∂F_{i}}{∂U_{j}} (X̄, Ū)\}_{i=1, j=1}^{n, m} , (8)

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

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

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

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

för den approximativa störningen x(t) ≈ ΔX(t) orsakad av datastörningarna x_{0} = ΔX_{0} och u(t) = ΔU(t).

Det linjära systemet (10) 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 (3). Observera att triviala lösningen x(t) = 0, u(t) = 0 till (10) motsvarar stationära lösningen X(t) = X̄, U(t) = Ū till (3).

Exempel. För pendeln blir Jacobimatriserna

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

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

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

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

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

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

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

vilket motsvarar

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

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

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 (12) 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);

Kommandot plot(x(:,1),x(:,2)) plottar istället 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 systemet (10). 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

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

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

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

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

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}. (18)

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 (18) framställer lösningen x(t) som en linjärkombination av egenmoder c_{k}e^{λ_{k}t}g_{k}. Karaktären hos dessa bestäms av egenvärdena λ_{k}. Om vi skriver λ_{k} = α_{k} + jω_{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} + jω_{k} är inte reellt (ω_{k} ≠ 0). Då är

e^{λ_{k}t} = e^{α_{k}t+jω_{k}t} = e^{α_{k}t}e^{jω_{k}t} = e^{α_{k}t}(cos ω_{k}t + j 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.

Bland annat kan vi nu dra slutsatser om lösningens beteende då t → ∞. Vi säger att systemet (16) är

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

• asymptotiskt stabilt om x(t) → 0 då t → ∞ för alla begynnelsevärden b,

• instabilt annars.

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 beror på tecknet på realdelen Re λ_{k} = α_{k}. Vi formulerar detta som en sats.

Sats 1. Om A är diagonaliserbar så är systemet x′ = Ax asymptotiskt stabilt om och endast om Re λ_{k} < 0 för alla egenvärden λ_{k}, 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.

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 (16) i

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

vilken är identisk med (17), 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

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

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

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

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.

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

\begin{aligned}& 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}, \\ (20) \end{aligned} ; \begin{aligned}& 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}. \\ (21) \end{aligned}

Antag nu att alla λ_{k} ≤ 0, dvs att A är negativt semidefinit (TLA 8.4). I så fall ordnar vi egenvärdena λ_{n} ≤ … ≤ λ_{k} ≤ … ≤ λ_{1} ≤ 0 så att e^{λ_{k}t} ≤ e^{λ_{1}t} ≤ 1 och

∥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} ≤ ∥b∥^{2}, (22)

dvs

∥x(t)∥ ≤ ∥b∥, t ≥ 0.

Vi drar slutsatsen att en störning x(t) inte kan bli större än begynnelsestörningen x(0) = b. Ett symmetriskt system är alltså stabilt om alla λ_{k} ≤ 0. Om A är negativt definit så har vi λ_{k} ≤ λ_{1} < 0 och

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

dvs störningen avtar exponentiellt då t → ∞; asymptotisk stabilitet.

Alltså: Om A är symmetrisk så är systemet x′ = Ax stabilt om och endast om λ_{k} ≤ 0 för alla egenvärden λ_{k}, asymptotiskt stabilt om och endast om λ_{k} < 0 för alla egenvärden λ_{k}, och instabilt om och endast om λ_{k} > 0 för något egenvärde λ_{k}.

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 (TLA 8.7). 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 (TLA 8.7).

En enkel räkning 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.

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 (10) av

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

där matrisen e^{tA} är lösningsoperatorn till det homogena systemet (16), se övning 4. 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(t) är liten.

Exempel. Vi återvänder nu till pendeln. Jacobimatrisen för det övre jämviktsläget är, se (14),

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 (12), 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

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

är asymptotiskt stabilt. 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

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

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 (24) är (MA2 9.8)

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

dvs vi har rekursionsekvationen

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

Vi får

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

Matrisen B = E + 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

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

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

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

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

Matrisen B = (E −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.

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 12.

2. Värmeledning, diffusion, reaktion

2.1. Härledning av värmeledningsekvationen. Vi studerar värmeledning i en kropp V ⊂ R^{3} med begränsningsytan S. Vi inför följande beteckningar:

T = T (x, y, z, t) temperaturen [K] i punkten (x, y, z) [m] vid tiden t[s],

j = j(x, y, z, t) värmeflödesvektorn [J/(m^{2}s)],

c_{p} = c_{p}(x, y, z) värmekapaciteten [J/(kgK)],

ρ = ρ(x, y, z) masstätheten [kg/m^{3}],

λ = λ(x, y, z) värmeledningsförmågan [J/(mKs)],

q = q(x, y, z, t) värmekällor/sänkor [J/(m^{3}s)].

Vi betraktar värmebalansen i en godtycklig delvolym V_{0} med begränsningsytan S_{0}. Totala värmeenergin är ∫_{V_{0}} ρc_{p}T dV [J] och nettoflödet ut genom S_{0} blir ∫_{S_{0}} j · n dS [J/s]. Här är n den utåtriktade enhetsnormalen till ytan S_{0}.

Energiprincipen säger att ändringen per tidsenhet av den totala värmeenergin i V_{0} är lika med värmeflödet in genom S_{0} plus värmeproduktionen per tidsenhet i V_{0}. Eftersom den totala värmeenergin är ∫_{V_{0}} ρc_{p}T dV och värmeproduktionen per tidsenhet är ∫_{V_{0}} q dV, får vi

\frac{d}{dt} ∫_{V_{0}} ρc_{p}T dV = − ∫_{S_{0}} j · n dS + ∫_{V_{0}} q dV.

Vi tillämpar nu divergenssatsen (MA3 4.14) för vektorfältet j och derivering under integraltecknet:

∫_{S_{0}} j · n dS = ∫_{V_{0}} ∇ · j dV och \frac{d}{dt} ∫_{V_{0}} ρc_{p}T dV = ∫_{V_{0}} ρc_{p} \frac{∂T}{∂t} dV,

vilket leder till

∫_{V_{0}} (ρc_{p} \frac{∂T}{∂t} + ∇ · j − q) dV = 0.

Eftersom V_{0} ⊂ V är godtycklig, så måste integranden vara identiskt lika med 0, dvs

ρc_{p} \frac{∂T}{∂t} + ∇ · j = q.

Fouriers lag säger att värmeflödet är proportionellt mot temperaturgradienten och att värme strömmar från varmt till kallt. Därför blir flödesvektorn

j = −λ∇T = −λ grad T

och vi får

ρc_{p} \frac{∂T}{∂t} − ∇ · (λ∇T) = q. (26)

Detta är värmeledningsekvationen.

Ekvation (26) måste kombineras med ett begynnelsevillkor

T (x, y, z, 0) = T_{0}(x, y, z) (27)

samt ett randvillkor, som uttrycker att värmeflödet ut genom randen är proportionellt mot temperaturdifferensen, j · n = k(T − T_{a}), där T_{a} är omgivningens temperatur och k [J/(m^{2} sK)] en värmeöverföringskoefficient. Eftersom

j · n = −λ∇T · n = −λ\frac{∂T}{∂n} på S,

där ∂T /∂n = ∇T · n betecknar normalderivatan av T, dvs riktningsderivatan i normalriktningen n, blir randvillkoret

(28) λ\frac{∂T}{∂n} + k(T − T_{a}) = 0 på S, (Robins randvillkor)

Gränsfallet k = 0 svarar mot att begränsningsytan är isolerad:

\frac{∂T}{∂n} = 0 på S. (Neumanns randvillkor)

Om vi å andra sidan dividerar med k i (28) och sedan låter k → ∞, får vi att temperaturen är föreskriven:

T = T_{a} på S. (Dirichlets randvillkor)

Gränsfallet k = ∞ betyder alltså att begränsningsytan är oisolerad, dvs värme kan obehindrat strömma genom ytan och kroppens yttemperatur blir då lika med omgivningens temperatur.

Det är ofta lämpligt att skriva om ekvationerna på dimensionslös form. Vi inför då dimensionslösa variabler:

t̃= t/τ, (x̃, ỹ, z̃) = (x/L, y/L, z/L)

u(x̃, ỹ, z̃, t̃) = \frac{T (x̃L, ỹL, z̃L, t̃τ) − T_{f}}{T_{f}} ,

där T_{f}, L, τ är referenskonstanter. Då fås

\frac{∂u}{∂t̃} = τ \frac{∂}{∂t}(\frac{T − T_{f}}{T_{f}} ), \frac{∂u}{∂x̃} = L\frac{∂}{∂x}(\frac{T − T_{f}}{T_{f}} ),

och så vidare. Om vi för enkelhets skull antar att koefficienterna ρ, c_{p}, λ är konstanta, och om vi dividerar (26) med λT_{f} får vi

\frac{ρc_{p}}{λ} \frac{∂}{∂t}(\frac{T − T_{f}}{T_{f}} ) − ∇ · ∇(\frac{T − T_{f}}{T_{f}} ) = \frac{q}{λT_{f}} ,

dvs

\frac{ρc_{p}}{λτ} \frac{∂u}{∂t̃} − L^{−2}∇ · ∇u = \frac{q}{λT_{f}} .

Med τ = L^{2}ρc_{p}/λ, f(x, y, z, t) = L^{2}q(x, y, z, t)/(λT_{f}) och med Laplaceoperatorn

Δu = ∇ · ∇u = \frac{∂^{2}u}{∂x̃^{2}} + \frac{∂^{2}u}{∂ỹ^{2}} + \frac{∂^{2}u}{∂z̃^{2}},

har vi värmeledningsekvationen med konstanta koefficienter:

\frac{∂u}{∂t̃} − Δu = f. (29)

Randvillkoret (28) och begynnelsevillkoret (27) blir på liknande vis:

(30) \frac{∂u}{∂n} + h(u − u_{a}) = 0 på S, (Robins randvillkor)

respektive

u(x, y, z, 0) = u_{0}(x, y, z). (31)

Här är h = Lk/λ och u_{a} = (T_{a} − T_{f})/T_{f} och motsvarande för u_{0}.

Den partiella differentialekvationen (29) tillsammans med begynnelsevillkoret (31) och randvillkoret (30) säges utgöra ett begynnelse-randvärdesproblem. I fortsättningen skriver vi för enkelhets skull t i stället för t̃ och så vidare.

Om f inte beror på t uppnås då t → ∞ ett stationärt tillstånd u = u(x, y, z). Eftersom nu ∂u/∂t = 0 får vi Poissons ekvation

−Δu = f.

Om dessutom f = 0 har vi Laplaces ekvation

−Δu = 0.

Dessa måste kombineras med lämpliga randvillkor enligt ovan. Vi har då ett så kallat randvärdesproblem.

Om vi betraktar värmeledning i en lång cylinder orienterad parallellt med x-axeln och vars mantelyta är isolerad, och där begynnelsetemperaturen u_{0} och källtermen f är oberoende av y och z, är det rimligt att anta att temperaturen också är oberoende av y och z, dvs u = u(x, t). Värmeledningsekvationen (29) blir då

u_{t} − u_{xx} = f.

Denna studeras i ett intervall t ex 0 < x < 1. Eftersom ∂u/∂n = −∂u/∂x vid x = 0 och ∂u/∂n = ∂u/∂x vid x = 1 blir randvillkoret (30) nu

−u_{x} + h_{0}(u − u_{a0}) = 0, x = 0; u_{x} + h_{1}(u − u_{a1}) = 0, x = 1,

där h_{0}, u_{a0} respektive h_{1}, u_{a1} är värmeöverföringskoefficienten och omgivningens temperatur vid ytan x = 0 respektive x = 1. Gränsfallen h_{0} = 0, h_{0} = ∞ etc tolkas analogt med ovan.

2.2. Diffusionsekvationen. Vid en diffusionsprocess är massflödet proportionellt mot gradienten av koncentrationen c = c(x, y, z, t) [mol/m^{3}] av det diffunderande ämnet (Ficks lag):

(32) j = −D∇c = −D grad c. [mol/(m^{2}s)]

Här är D [m^{2}/s] diffusionskoefficienten. Samma härledning som ovan leder till diffusionsekvationen:

\frac{∂c}{∂t} − ∇ · (D∇c) = f, [mol/(m^{3}s)] (33)

eller, om D är konstant,

\frac{∂c}{∂t} − DΔc = f. (34)

Matematiskt sett är detta samma ekvation som värmeledningsekvationen. Randvillkoret blir

(35) D \frac{∂c}{∂n} + d(c − c_{a}) = 0 på S, (Robins randvillkor)

där koefficienten d [m/s] beror på ytans permeabilitet och c_{a} är concentrationen utanför S.

2.3. Reaktions-diffusionsekvationer. Betrakta en enkel exoterm kemisk reaktion där reaktionshastigheten ges av Arrhenius lag k_{0}c exp (−E/RT) [mol/(m^{3}s)]. Här är E [J/mol] aktiveringsenergin, R [ 8.31 J/(mol K)] allmänna gaskonstanten, T [K] temperaturen och c [mol/m^{3}] koncentrationen av det reagerande ämnet. Reaktionsentalpin är −ΔH [J/mol]. Vi får två kopplade partiella differentialekvationer

\begin{aligned}\frac{∂T}{∂t} − κΔT & = βk_{0}ce^{−E/RT} , β = \frac{−ΔH}{ρc_{p}} , κ = \frac{λ}{ρc_{p}} , \\ \frac{∂c}{∂t} − DΔc & = −k_{0}ce^{−E/RT} , \end{aligned} (36)

med randvillkoren

\frac{∂T}{∂n} + h(T − T_{0}) = 0, \frac{∂c}{∂n} + d(c − c_{0}) = 0 på S. (37)

Observera att differentialekvationerna är ickelinjära och kopplade till varandra.

3. Övningar

Svar och lösningar till vissa övningar finns i nästa avsnitt.

Övning 1. Bestäm lokala maxima och minima till funktionen f(x, y) = x^{3} + y^{2} − xy.

Övning 2. Kvadratisk form. (TLA 8.4, Matlab) Bestäm största och minsta värdet av den kvadratiska formen

f(x, y, z) = 9x^{2} + 5y^{2} + 6z^{2} + 8xy + xz + 3yz

på sfären x^{2} + y^{2} + z^{2} = 1 och ange i vilka punkter maximum och minimum inträffar.

Övning 3. Minimering. (Matlab) Sök minimum av funktionen

f(x) = x_{1}^{4} + x_{2}^{4} − 4x_{1}x_{2} + x_{1} + 8, x = [x_{1}, x_{2}]^{t}

(Tips: använd kommandot fmins; skriv » help fmins för dokumentation. Detta kommando fordrar att man skriver en funktionsfil (kallad t ex funk.m) som definierar funktionen f, se TLA sid App:27.)

Övning 4. Bevisa (23) genom att multiplicera ekvation (10) med den integrerande faktorn e^{−tA} (se (19)) och integrera.

Övning 5. 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 6. 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 7. 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 8. Bestäm de stationära punkterna och avgör deras stabilitet för följande system:

\begin{aligned}(a) X′ & = X(1 − Y) \\ Y′ & = Y (1 − X) \end{aligned}  \begin{aligned}(b) X′ & = −2(X − 10) + Y e^{X} \\ Y′ & = −2Y − Y e^{X} \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 9. 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 10. (Repetition av Laplacetransformen, (MA2 kap 15).) Lös differentialekvationen i övning 9 med Laplacetransformen.

Övning 11. Lös ekvationerna (13) och (15) med α = 0 och β = 0.

Övning 12. 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.

Övning 13. Laplaceoperatorn i sfärisk symmetri. Inför sfäriska koordinater (r, θ, ϕ) (MA3, sid 3.26):

x = r sin θ cos ϕ
y = r sin θ sin ϕ
z = r cos θ

Antag att funktionen u ej beror på θ och ϕ, dvs u = u(r). Visa att

Δu = \frac{1}{r^{2}} \frac{∂}{∂r} (r^{2} \frac{∂u}{∂r})

Övning 14. Laplaceoperatorn i cylindrisk symmetri. Inför cylindriska koordinater (ρ, φ, z) (MA3, sid 2.41):

x = ρ cos φ
y = ρ sin φ
z = z

Antag att funktionen u ej beror på φ, dvs u = u(ρ, z). Visa att

Δu = \frac{1}{ρ} \frac{∂}{∂ρ}(ρ\frac{∂u}{∂ρ}) + \frac{∂^{2}u}{∂z^{2}}

Övning 15. Visa att följande funktioner satisfierar Laplaces ekvation Δu = 0.

(a) u(r, θ, ϕ) = \frac{1}{r} (sfäriska koordinater)

(b) u(ρ, φ, z) = ln ρ + z (cylindriska koordinater)

(c) u(x, y, z) = 2x^{2} − y^{2} − z^{2}

Övning 16. Lös följande randvärdesproblem. Tolka lösningen fysikaliskt i termer av värmeledning eller diffusion samt rita dess graf.

(a) − u′′_{xx} = 1, 0 < x < L, u(0) = u(L) = 0

(b) − u′′_{xx} = 1, 0 < x < L, u(0) = 0, u′(L) = 0

(c) − u′′_{xx} = 0, 0 < x < L, −u′(0) + \frac{2}{L}(u(0) − 1) = 0, u′(L) + \frac{3}{L}(u(L) − 2) = 0

(d) − u′′_{xx} = 1, 0 < x < L, u′(0) = u′(L) = 0

Övning 17. Visa att funktionen u(x, y, z) = sin x sin y sin z safisfierar randvärdesproblemet

−Δu = 3u i K, u = 0 på S

där K är kuben {0 < x < π, 0 < y < π, 0 < z < π} med randen S. Alltså är u en egenfunktion till Laplaceoperatorn i K med egenvärde 3.

Övning 18. Betrakta stationär värmeledning genom en vägg:

−(d(x)u′(x))′ = 0, −L < x < L; u(−L) = g_{0}, u(L) = g_{1}.

Yttertemperaturerna g_{0}, g_{1} samt värmeledningsförmågan d(x) i isoleringsmaterialet är givna. Uppgiften är att bestämma en effektiv värmeöverföringskoefficient k sådan att värmeflödet genom väggen kan skrivas j = k(g_{0} − g_{1}).

(a) Antag att d(x) = ϵ(1 + (βx)^{2}) och lös problemet för hand.

(b) Lös randvärdesproblemet med Adfem. För att bestämma värmeflödet j(x) = −d(x)u′(x) gör vi så här: Nodvärdena av den beräknade temperaturen finns i vektorn u, nätstorleken finns i vektorn h. Om vi bildar differensen mellan konsekutiva u-värden med programmet diff (» help diff) och dividerar med h-värdena, så får vi temperatur-derivatan u′. Den är en styckvis konstant funktion och fås alltså i vektorn diff(u)./h. Värdena av värmeledningskoefficienten i mittpunkterna på elementen finns i vektorn d. Värmeflödet ges alltså av j=-d.*diff(u)./h. Testa först med ett specialfall där exakta lösningen är känd, dvs välj d(x) så att u(x) kan beräknas symboliskt som t ex i (a). Prova sedan med en annan värmeledningskoefficient, t ex d(x) = 1.01 + cos 7x. (Försök inte bestämma den exakta lösningen i detta fall.) Detta problem är svårare att lösa även för Adfem. Varför?

Övning 19. Reaktion-diffusion i porös katalysator-pellet. Ett ämne i gasform diffunderar in i en porös katalysator-pellet och omsätts vid isoterma och stationära betingelser. Vi får reaktionsdiffusions ekvationen, se (36),

\begin{matrix}− DΔc = −kc i V , \\ c = c_{f} på S, \end{matrix} (38)

med koncentrationen c [mol/m^{3}] i katalysatorn, den konstanta koncentrationen c_{f} i omgivningen, diffusionskonstanten D [m^{2}/s], reaktionshastigheten k [s^{−1}], storleken L [m], och där Δ = \frac{∂^{2}}{∂x^{2}} + \frac{∂^{2}}{∂y^{2}} + \frac{∂^{2}}{∂z^{2}} .

PSfrag replacements X(t)

X_{0}
U(t)

c

c_{f}
V

S L

Figur 5. Pellet.

(a) Inför dimensionslösa variabler (x^{∗}, y^{∗}, z^{∗}) = (x/L, y/L, z/L), c^{∗} = c/c_{f} så att (38) övergår i

\begin{matrix}− Δ^{∗}c^{∗} = −ϕ^{2}c^{∗} i V^{∗}, \\ c^{∗} = 1 på S^{∗}, \end{matrix} (39)

där Δ^{∗} = \frac{∂^{2}}{∂x^{∗2}} + \frac{∂^{2}}{∂y^{∗2}} + \frac{∂^{2}}{∂z^{∗2}} och ϕ kallas Thiele-modulen. Uttryck ϕ i de givna storheterna.

PSfrag replacements X(t) X_{0} U(t) c c_{f} V S L c^{∗} 1 V^{∗} S^{∗} 1

Figur 6. Pellet i dimensionlösa variabler.

(b) För olika värden på ϕ vill vi undersöka dels hur långt in i katalysatorn reaktionen når, dels beräkna kvoten mellan den uppnådda reaktionshastigheten och den ideala reaktionshastighet man skulle få utan hindrande diffusion, dvs om c = c_{f}. Denna kvot kallas effektivitetsfaktorn och ges alltså av

η = \frac{∫_{V} kc dV}{kc_{f}V} = \frac{∫_{V} c dV}{c_{f}V} = \frac{1}{V^{∗}} ∫_{V^{∗}} c^{∗} dV^{∗}. (40)

Visa detta! Visa också att

η = \frac{1}{ϕ^{2}} \frac{1}{V^{∗}} ∫_{S^{∗}} \frac{∂c^{∗}}{∂n^{∗}} dS^{∗}. (41)

Tips: (39) och Gauss sats.

(c) Speciell geometri: (i) platta (eller stav med tät mantelyta)

\begin{matrix}− \frac{d^{2}c^{∗}}{dr^{2}} = −ϕ^{2}c^{∗}, −1 < r < 1, \\ c^{∗}(−1) = c^{∗}(1) = 1, \end{matrix} (42)

eller på grund av symmetri (c^{∗} är en jämn funktion, dvs c^{∗}(−r) = c^{∗}(r), och för varje jämn funktion gäller att derivatan är udda och därför är derivatan = 0 i origo)

\begin{matrix}− \frac{d^{2}c^{∗}}{dr^{2}} = −ϕ^{2}c^{∗}, 0 < r < 1, \\ \frac{dc^{∗}}{dr} (0) = 0, c^{∗}(1) = 1. \end{matrix} (43)

16 STIG LARSSON

PSfrag replacements X(t) X_{0} U(t) c c_{f} V S L c^{∗} 1 V^{∗} S^{∗} 1 −1 1 0 r

Figur 7. Platta eller stav.

(ii) cirkulär cylinder

\begin{matrix}− \frac{1}{r} \frac{d}{dr} (r \frac{dc^{∗}}{dr} ) = −ϕ^{2}c^{∗}, 0 < r < 1, \\ \frac{dc^{∗}}{dr} (0) = 0, c^{∗}(1) = 1. \end{matrix} (44)

(iii) sfär

\begin{matrix}− \frac{1}{r^{2}} \frac{d}{dr} (r^{2} \frac{dc^{∗}}{dr} ) = −ϕ^{2}c^{∗}, 0 < r < 1, \\ \frac{dc^{∗}}{dr} (0) = 0, c^{∗}(1) = 1. \end{matrix} (45)

(Den dimensionlösa variabeln r borde egentligen tecknas r^{∗}, men det blir alltför svårläst.)

PSfrag replacements X(t) X_{0} U(t) c c_{f} V S L c^{∗} 1 V^{∗} S^{∗} 1 −1 1 0 r 1 0 r

Figur 8. Cylinder och sfär.

(d) Lös (43) exakt och bestäm och plotta η som funktion av ϕ.

(e) Lös (45) exakt och bestäm och plotta η som funktion av ϕ. (Se lösning i facit.)

(f) (Adfem, Matlab) Lös (43), (44), (45) numeriskt med programmet Adfem i Matlab och bestäm η för några olika värden på ϕ (små och stora). Jämför med exakta lösningar i (d) och (e). Vilken geometri är fördelaktigast? Tips: Funktionsvärdena c^{∗}(r_{i}) till lösningen finns i vektorn u, beräkningsnoderna r_{i} i vektorn x och nätstorleken h_{i} i vektorn h. Kommandot integrate(u.*x,h) beräknar integralen ∫_{0}^{1} c^{∗}(r)r dr.

(g) Varianter för den ambitiöse teknologen. Gör sammalunda för någon av följande modifikationer. (i) Robins randvillkor: byt ut randvillkoret i (38) mot villkoret (35). Visa att i dimensionslösa variabler får vi

\frac{dc^{∗}}{dr} (0) = 0, \frac{dc^{∗}}{dr} (1) + ν(c^{∗}(1) − 1) = 0.

Uttryck ν i de ursprungliga variablerna. Plotta η som funktion av ν för fixt ϕ.

(ii) Tät kärna:

\begin{matrix}− \frac{1}{r^{2}} \frac{d}{dr} (r^{2} \frac{dc^{∗}}{dr} ) = −ϕ^{2}c^{∗}, ω < r < 1, \\ \frac{dc^{∗}}{dr} (ω) = 0, c^{∗}(1) = 1. \end{matrix} (46)

(iii) Förgiftning:

\begin{matrix}− \frac{1}{r^{2}} \frac{d}{dr} (r^{2} \frac{dc^{∗}}{dr} ) = −ϕ^{2}c^{∗}, 0 < r < ω, \\ − \frac{1}{r^{2}} \frac{d}{dr} (r^{2} \frac{dc^{∗}}{dr} ) = 0, ω < r < 1, \\ \frac{dc^{∗}}{dr} (0) = 0, c^{∗}(1) = 1. \end{matrix} (47)

4. Svar och lösningar till vissa övningar

1. De stationära punkterna till funktionen f(x, y) = x^{3} + y^{2} − xy ges av ekvationssystemet

f′_{x}(x, y) = 3x^{2} − y = 0 (48)

f′_{y}(x, y) = 2y − x = 0 (49)

Den andra ekvationen ger x = 2y vilket insättes i den första ekvationen: 12y^{2} − y = 0, vilken har två lösningar y = 0 och y = 1/12. Vi har alltså två stationära punkter nämligen (0, 0) och

(1/6, 1/12). Hessematrisen är f′′(x, y) = [\begin{matrix}f′′_{xx} & f′′_{xy} \\ f′′_{xy} & f′′_{yy} \end{matrix}] = [\begin{matrix}6x & −1 \\ −1 & 2 \end{matrix}]. I de två punkterna får vi

f′′(0, 0) = [\begin{matrix}0 & −1 \\ −1 & 2 \end{matrix}] med egenvärdena 1 ± \sqrt{2}, så att matrisen är indefinit, och f′′(1/6, 1/12) =

[\begin{matrix}1 & −1 \\ −1 & 2 \end{matrix}] med egenvärdena (3 ± \sqrt{5})/2, så att matrisen är positivt definit. Vi drar slutsatsen att

(0, 0) är en sadelpunkt och att (1/6, 1/12) är en minimipunkt.

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

6. 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.

8. (a) (0, 0) instabil, (1, 1) instabil; (b) (10, 0) asymptotiskt stabil; (c) (0, 0, 0) och (0, −1, 1) båda instabila 10. −\frac{1}{3}e^{−2t} sin 3t

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

16. (a) \frac{1}{2}L^{2}(x/L)(1−(x/L)). (b) L^{2}(x/L)(1−\frac{1}{2}(x/L)). (c) \frac{6}{11}(x/L)+\frac{14}{11}. (d) Det finns ingen lösning. Problemet är inte “rätt ställt” (not “well posed”). Detta värmeledningsproblem är ofysikaliskt och saknar stationär lösning. Man pumpar in värme men släpper inte ut någon.

18. (a) Genom att integrera differentialekvationen får vi

d(x)u′(x) = A, u′(x) = \frac{A}{d(x)} = \frac{A}{ϵ} \frac{1}{1 + (βx)^{2}} , u(x) = \frac{A}{ϵβ} arctan(βx) + B,

g_{0} = u(−L) = \frac{A}{ϵβ} arctan(−βL) + B = −\frac{A}{ϵβ} arctan(βL) + B,

g_{1} = u(L) = \frac{A}{ϵβ} arctan(βL) + B,

A = ϵβ \frac{g_{1} − g_{0}}{2 arctan(βL)}, B = \frac{g_{1} + g_{0}}{2} ,

u(x) = \frac{g_{1} − g_{0}}{2} \frac{arctan(βx)}{arctan(βL)} + \frac{g_{1} + g_{0}}{2} ,

j(x) = −d(x)u′(x) = −A = −ϵβ \frac{g_{1} − g_{0}}{2 arctan(βL)}, k = \frac{ϵβ}{2 arctan(βL)}.

19. (a) ϕ = L\sqrt{k/D}. (d) c^{∗}(r) = cosh(ϕr)/ cosh(ϕ), η = tanh(ϕ)/ϕ, c^{∗}(0) = 1/ cosh(ϕ). Obs att c^{∗}(r) ≡ 1 om ϕ = 0.

(e) Med u = rc^{∗}(r) övergår (45) i

u′′(r) = ϕ^{2}u(r),

med allmänna lösningen

u(r) = A cosh(ϕr) + B sinh(ϕr),

så att

c^{∗}(r) = \frac{1}{r}(A cosh(ϕr) + B sinh(ϕr)),
\frac{dc^{∗}(r)}{dr} = −\frac{1}{r^{2}} (A cosh(ϕr) + B sinh(ϕr)) + \frac{ϕ}{r}(A sinh(ϕr) + B cosh(ϕr))

= Aϕ^{2} \frac{ϕr sinh(ϕr) − cosh(ϕr)}{(ϕr)^{2}} + Bϕ^{2} \frac{ϕr cosh(ϕr) − sinh(ϕr)}{(ϕr)^{2}}.

Taylorutvecklingarna cosh x = 1 + \frac{1}{2}x^{2} + O(x^{4}), sinh x = x + \frac{1}{6}x^{3} + O(x^{5}) visar att uttrycket efter A (i derivatan) går mot ∞, och uttrycket efter B går mot 0, då r går mot 0. Därför måste A = 0 och randvillkoret c^{∗}(1) = 1 ger sedan B = 1/ sinh(ϕ). Lösningen blir

c^{∗}(r) = \frac{sinh(ϕr)}{r sinh(ϕ)}.

Effektivitetsfaktorn blir enligt (40)

η = \frac{1}{V^{∗}} ∫_{V^{∗}} c^{∗} dV^{∗} = \frac{1}{4π/3} ∫_{0}^{1} c^{∗}(r)4πr^{2} dr = 3 ∫_{0}^{1} c^{∗}(r)r^{2} dr

= \frac{3}{sinh(ϕ)} ∫_{0}^{1} r sinh(ϕr) dr = \frac{3}{ϕ^{2} sinh(ϕ)} ∫_{0}^{ϕ} s sinh(s) ds

= \{ partiell integration \} = \frac{3}{ϕ^{2}} (\frac{ϕ}{tanh(ϕ)} − 1).

Eller enligt (41)

η = \frac{1}{ϕ^{2}} \frac{1}{V^{∗}} ∫_{S^{∗}} \frac{∂c^{∗}}{∂n^{∗}} dS^{∗} = \frac{1}{ϕ^{2}} \frac{1}{4π/3}c^{∗}′(1)4π

= \frac{3}{ϕ^{2}} \frac{ϕ cosh(ϕ) − sinh(ϕ)}{sinh(ϕ)} = \frac{3}{ϕ^{2}} (\frac{ϕ}{tanh(ϕ)} − 1).

En intressant observation: Eftersom sinh(x)/x → 1 då x → 0, och

c^{∗}(r) = \frac{ϕ}{sinh(ϕ)} \frac{sinh(ϕr)}{ϕr},

ser vi att c^{∗}(0) = ϕ/ sinh(ϕ), och att c^{∗}(r) → 1 då ϕ → 0 med r fixt, vilket stämmer överens med att c^{∗}(r) ≡ 1 om ϕ = 0.

5. Tankreaktorn

I en tankreaktor på 1.2 m^{3} skall ett ämne (c_{f} = 5 kmol/m^{3}) omsättas till 50% enligt en första ordningens reaktion med k = 2.8 · 10^{−5} s^{−1} vid 70 ^{∘}C. Mätningar visar att värdet på k fördubblas för var tionde ^{∘}C. Reaktionsvärmet är ΔH = −83.7 kJ/mol. Värmekapaciteten och densiteten hos materialströmmen varierar ej med temperaturen och sammansättningen, och är c_{p} = 4.19 kJ/(kg K) respektive ρ = 1000 kg/m^{3}. Flödeshastigheten är q = 2.79 · 10^{−4} m^{3}/s. Inflödets temperatur är T_{f} = 70^{∘}C. För värmeväxlaren i tanken kan man räkna med värmegenomgångskoefficienten κ = 0.58 kJ/(m^{2} s K).

PSfrag replacements X(t)

X_{0}
U(t)
c
c_{f}
q
T
T_{f}
T_{K}
A

Figur 9. Tankreaktor.

Beräkna gränsvärdet på erforderlig värmeväxlararea A för stabil drift. Välj därefter lämplig area och beräkna temperaturen T_{K} i kylvattnet vid jämvikt vid 50% omsättningsgrad. (Omsättningsgraden= (c_{f} − c)/c_{f}.) (Från “Kemisk reaktionsteknik. Övningsuppgifter”, Kemisk reaktionsteknik, CTH 1993, uppgift 7.1.)

5.1. Lösning. De använda m-filerna finns på kemi-pc och på kursens web-sida, se

  http://www.math.chalmers.se/~stig

På kemi-pc blir filerna tillgängliga efter att man exekverat kommandot “tillmat” i Matlab.

(a) Vi ställer först upp balansekvationer för massa och energi. Omsättningshastigheten av massa är V ck [mol/s], där reaktionshastigheten ges av

k = k_{0} exp(−E/(RT)) [s^{−1}]

enligt Arrhenius lag. Värmeproduktionshastigheten blir då Q = (−ΔH)V ck [J/s]. Balansekvationerna för massa V c [mol] och värmeenergi ρc_{p}V T [J] blir

\frac{d}{dt}(V c) = q(c_{f} − c) − V ck_{0} exp(−E/(RT)),
\frac{d}{dt}(ρc_{p}V T) = ρc_{p}q(T_{f} − T) + (−ΔH)V ck_{0} exp(−E/(RT)) − κA(T − T_{K}).

(b) För att kunna räkna (numeriskt) måste man skala variablerna så att de får rimliga storleksordningar (det visar sig t ex att exp(−E/(RT)) ≈ e^{−30}). Vi inför därför dimensionslösa variabler enligt

s = t/τ, τ = V/q [s],
X_{1}(s) = c(sτ)/c_{f}, X_{2}(s) = T (sτ)/T_{f}, U(s) = T_{K}(sτ)/T_{f},

där τ = V/q [s] kallas uppehållstiden.

Den första ekvationen divideras med qc_{f}, den andra med ρc_{p}qT_{f}. Med τ = V/q får vi

τ \frac{d}{dt}(\frac{c}{c_{f}} ) = 1 − \frac{c}{c_{f}} − \frac{c}{c_{f}} τk_{0} exp (−\frac{E}{RT_{f}} \frac{T_{f}}{T} ),
τ \frac{d}{dt}(\frac{T}{T_{f}} ) = 1 − \frac{T}{T_{f}} + \frac{(−ΔH)c_{f}}{ρc_{p}T_{f}} \frac{c}{c_{f}} τk_{0} exp (−\frac{E}{RT_{f}} \frac{T_{f}}{T} ) − \frac{κAτ}{ρc_{p}V} (\frac{T}{T_{f}} − \frac{T_{K}}{T_{f}} ).

Med de dimensionslösa variablerna

s = t/τ, X_{1}(s) = \frac{c(sτ)}{c_{f}}, X_{2}(s) = \frac{T (sτ)}{T_{f}}, U(s) = \frac{T_{K}(sτ)}{T_{f}},
α = \frac{(−ΔH)c_{f}}{ρc_{p}T_{f}} , β = \frac{κAτ}{ρc_{p}V} , γ = \frac{E}{RT_{f}} , δ = τk_{0}e^{−γ}, f(x) = δe^{γ−γ/x},

får vi det ickelinjära differentialekvationssystemet

\frac{dX_{1}}{ds} = 1 − X_{1} − X_{1}f(X_{2})  ( = F_{1}(X, U) ),
\frac{dX_{2}}{ds} = 1 − X_{2} + αX_{1}f(X_{2}) − β(X_{2} − U)  ( = F_{2}(X, U) ).

Talet γ = E/(RT_{f}) kallas Arrhenius tal. Observera att e^{γ−γ/X_{2}} ≈ 1 då X_{2} ≈ 1, vilket är en hanterlig storleksordning. (Alternativt kan man använda variablerna Y_{1} = 1 − X_{1} = (c_{f} − c)/c_{f} (omsättningsgraden) och Y_{2} = γ(X_{2}−1) = γ(T −T_{f})/T_{f}. Då blir e^{γ−γ/X_{2}} = e^{Y_{2}/(1+Y_{2}/γ)} ≈ e^{Y_{2}} ≈ 1 då Y_{2} ≈ 0.)

(c) Vi bestämmer parametrarna δ, γ genom att anpassa funktionen

τk = δe^{γ−γ/X_{2}}

till givna mätdata för k och T. Dessa finns i filen data.m. När man exekverar denna läses alla data in. Se exemplet “fitting” i Matlabs demo. Vi använder Matlab-programmet fmins och filen fel.m. Mätdata plottas tillsammans med den anpassade funktionen för att se hur bra det stämmer.

Med givna data får man δ = 0.0761, γ = 30.9896. Detta betyder att reaktionshastighetskonstanten k_{0} = e^{γ}δ/τ ≈ 5.0883 · 10^{8} [s^{−1}] och aktiveringsenergin E = γRT_{f} ≈ 8.8376 · 10^{4} [J/mol] (gaskonstanten R ≈ 8.3143 [J/molK]). Beräkningarna görs av programmet step1.m.

Anmärkning: Det visar sig vara omöjligt att anpassa Arrhenius lag i sin ursprungliga form k = k_{0} exp(−E/(RT)) till mätdata för k och T. Den numeriska algoritmen i fmins konvergerar inte. Det är orsaken till att vi skriver om Arrhenius lag på formen kτ = δe^{γ−γ/X}.

(d) Kravet att omsättningen ska vara 50% innebär att arbetspunkten ska vara c̄ = 0.5c_{f}, dvs X̄_{1} = 0.5. Vi bestämmer X̄_{2} och Ū så att X̄_{1} = 0.5. Dessa uttrycks i α, β, γ, δ, dvs numeriska värden sätts inte in ännu.

Stationära lösningar ges av

1 − X̄_{1} − X̄_{1}f(X̄_{2}) = 0,
1 − X̄_{2} + αX̄_{1}f(X̄_{2}) − β(X̄_{2} − Ū) = 0.

Om X̄_{1} är fixerad ger den första ekvationen

f(X̄_{2}) = (1 − X̄_{1})/X̄_{1},

dvs

X̄_{2} = \frac{γ}{γ + ln(δX̄_{1}/(1 − X̄_{1}))}.

Den andra ekvationen ger sedan

Ū = X̄_{2} − \frac{1 − X̄_{2} + αX̄_{1}f(X̄_{2})}{β} = X̄_{2} − \frac{1 − X̄_{2} + α(1 − X̄_{1})}{β}.

(e) Nu linjäriserar vi kring X̄, Ū. Det linjäriserade systemet blir

x′(s) = Ax(s) + Bu(s), s > 0; x(0) = x_{0},

med Jacobimatriserna

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}−1 − f(X̄_{2}) & −X̄_{1}f′(X̄_{2}) \\ αf(X̄_{2}) & −1 + αX̄_{1}f′(X̄_{2}) − β \end{matrix}] ,

där f′(x) = \frac{γ}{x^{2}} f(x), och

B = [\begin{matrix}\frac{∂F_{1}}{∂U} (X̄, Ū) \\ \frac{∂F_{2}}{∂U} (X̄, Ū) \end{matrix} ] = [\begin{matrix}0 \\ β \end{matrix}] .

(f) Nu bestämmer vi en undre gräns β_{0} för β sådan att det linjäriserade systemet blir stabilt för β > β_{0}. Vi gör detta genom att systematiskt beräkna egenvärdena till Jacobimatrisen numeriskt för olika värden på β.

Vi sätter X̄_{1} = 0.5 och beräknar att f(X̄_{2}) = 1, X̄_{2} = 1.0906, α = 0.2912. Vi beräknar sedan A och Ū för olika värden på β. Numeriska experiment visar att A har två egenvärden λ_{1} < 0, λ_{2} > 0 om 0 < β < β_{0} ≈ 0.898; ett konjugerat par av komplexa egenvärden med Re λ_{1,2} < 0 om β_{0} < β < β_{1} ≈ 8.688; och två negativa egenvärden om β > β_{1}. Arbetspunkten X̄ = (0.5, 1.0906) är alltså stabil om β > β_{0} ≈ 0.898 och instabil annars. Beräkningarna görs av programmet step2.m. Där görs en parameterstudie: β varieras i ett intervall och real- och imaginärdelarna av egenvärdena plottas som funktion av β.

(g) Här löser vi numeriskt både det linjäriserade och det icke-linjära systemet med små störningsdata x_{0} och u(s) = u_{0} cos(s) kring X̄, Ū. Vi väljer β både under och över stabilitetsgränsen.

Beräkningarna görs av programmen step3.m och step4.m. Vi använder data x(0) = x_{0}, u(s) = u_{0} cos(s) för det linjäriserade systemet och X(0) = X̄ + x_{0}, U(s) = Ū + u(s) för det ickelinjära. Vi tar först u_{0} = 0 och t ex x_{0} = (0, 0), x_{0} = 0.001 (1, ±1), och sedan t ex u_{0} = 0.001. Man ser tydligt att arbetspunkten X̄_{1} = 0.5, X̄_{2} = 1.0906 är stabil om β > β_{0} ≈ 0.898 och instabil annars. I det senare fallet tippar systemet över till en annan (icke önskvärd) stabil stationär punkt. Det verkar finnas två sådana.

(h) Slutligen ger vi en rekommendation om hur reaktorn ska köras uttryckt i de ursprungliga variablerna c, T, T_{K}, A. Vi beräknar även hur mycket värme som måste kylas bort.

Den kritiska arean blir

A_{0} = \frac{β_{0}ρc_{p}V}{κτ} = 1.8099 m^{2}.

För stabil drift väljer vi t ex β = 3, vilket ger arbetspunkten

c = X̄_{1}c_{f} = 0.5c_{f} = 2500 mol/m^{3},

T_{f} = 343.15 K = 70.00 ^{∘}C,

T = X̄_{2}T_{f} = 374.24 K = 101.09 ^{∘}C,

T_{K} = ŪT_{f} = 367.95 K = 94.80 ^{∘}C, T − T_{K} = 6.28 K = 6.28 ^{∘}C,

A = \frac{βρc_{p}V}{κτ} = 6.047 m^{2},

kyleffekten = κA(T − T_{K}) = 22.04 kJ/s.

Om man vill ha ett mycket snabbt insvängningsförlopp vid start av reaktorn, bör man välja β > β_{1}, t ex β = 10.

Beräkningarna görs av programmet step5.m.

Anmärkning: I “Kemisk reaktionsteknik. Övningsuppgifter”, CTH 1993, uppgift 7.1, uppskattas kritiska arean till 1.5 m^{2}, vilket svarar mot β = 0.7442. För att göra en jämförelse beräknar vi Q = (−ΔH)kcV = ρc_{p}qT_{f}αX_{1}f(X_{2}) [J/s] och i arbetspunkten

\frac{∂Q}{∂T} = \frac{∂Q}{∂X_{2}} \frac{∂X_{2}}{∂T} = ρc_{p}qαX̄_{1}f′(X̄_{2}) = ρc_{p}qαX̄_{1} \frac{γ}{X̄_{2}^{2}} f(X̄_{2}) = 4.43 · 10^{3} J/(s K).

I kompendiet uppskattas ∂Q/∂T till 2.03·10^{3} J/(s K). Detta leder till en underskattning av kritiska arean. Felet verkar vara att man beräknar

\frac{∂Q}{∂T}|_{c_{f} −c−ckτ=konst} i stället för \frac{∂Q}{∂T}|_{c=konst}.