EExpert · deel E1 · Wiskundige gereedschapskist

Oplossingen · hoofdstuk 4Differentiaalvergelijkingen en simulatie

Wie weet hoe snel iets verandert, kan uitrekenen hoe het verder gaat: exact of stap voor stap

Uitgewerkte oplossingen bij de 16 opdrachten van het werkboek. Voor de leerkracht, en voor wie zichzelf wil verbeteren nadat hij het eerst zelf probeerde.

5× oefenen6× toepassen5× uitdagen± 480 min

WWiskundig modelleren PProgrammeren en computationeel denken EEnergie BLeven en systemen

Wat de sensor in de meetkas meet, kan het brein van het station nu ook voorspellen: met een model van hoe snel de kas opwarmt en afkoelt.

Voor de leerkrachtDeze pagina bevat de antwoorden. Druk het werkboek af zonder deze pagina.

Oefenen

één stap, direct toepassen

4.1

Een oplossing controleren

Oefenen6Rekenen15 min

Controleer telkens of de functie een oplossing is van de vergelijking én voldoet aan de beginwaarde. Toon je berekening.

a
$y (t) = 3 e^{- 2 t}$ voor $\overset{y}{˙} = - 2 y$ , $y (0) = 3$ .
$\overset{y}{˙} = 3 \cdot (- 2) e^{- 2 t} = - 2 \cdot 3 e^{- 2 t} = - 2 y$ en $y (0) = 3$ : ja.
b
$u (t) = 12 (1 - e^{- 50 t})$ voor $0, 02 \overset{u}{˙} = 12 - u$ , $u (0) = 0$ . Wat is $τ$ ?
$\overset{u}{˙} = 12 \cdot 50 e^{- 50 t} = 600 e^{- 50 t}$ , dus $0, 02 \overset{u}{˙} = 12 e^{- 50 t}$ ; en $12 - u = 12 e^{- 50 t}$ . Gelijk, en $u (0) = 0$ : ja. $τ = 0, 02 s = 20 ms$ (en $1/ τ = 50$ ).
c
$T (t) = 20 + 5 e^{- t}$ voor $\dot{T} = - (T - 20)$ , $T (0) = 30$ .
$\dot{T} = - 5 e^{- t} = - (T - 20)$ : de vergelijking klopt, maar $T (0) = 25 \neq = 30$ . Geen oplossing van dit beginwaardeprobleem; de juiste is $20 + 10 e^{- t}$ .
d
$N (t) = \frac{10}{1 + 9 e ^{- t}}$ voor $\dot{N} = N (1 - N /10)$ , $N (0) = 1$ .
$\dot{N} = \frac{10 \cdot 9 e ^{- t}}{( 1 + 9 e ^{- t} ) ^{2}}$ (kettingregel). Rechterlid: $N (1 - N /10) = \frac{10}{1 + 9 e ^{- t}} \cdot \frac{9 e ^{- t}}{1 + 9 e ^{- t}}$ : hetzelfde. $N (0) = 10/10 = 1$ : ja (logistisch met $r = 1$ , $K = 10$ ).

4.2

Van woorden naar een vergelijking

Oefenen6Begrijpen20 min

Schrijf voor elke situatie een differentiaalvergelijking. Geef de eenheid van elke constante, het evenwicht, en zeg of de vergelijking lineair is.

a
De batterij van de meetkas verliest in rust 3 % van haar lading per maand (zelfontlading).
$d Q / d t = - k Q$ met $k = - ln 0, 97 = 0, 0305$ per maand (ongeveer 3 % per maand); evenwicht $Q = 0$ (leeg); lineair.
b
Een regenton loopt leeg door een kraantje onderaan; het debiet is evenredig met de wortel uit de waterhoogte h (Torricelli).
$A d h / d t = - c h$ , dus $d h / d t = - k h$ met $k$ in $m^{1/2} /s$ ; evenwicht $h = 0$ ; niet lineair ( $h$ ). Scheiden geeft $h$ lineair dalend in de tijd: de ton is na een eindige tijd leeg.
c
Een fles uit de koelkast (4 °C) warmt op in een kamer van 21 °C met een tijdconstante van 40 minuten.
$d T / d t = (21 - T) /40$ ( $t$ in min); evenwicht 21 °C (stabiel); lineair. Oplossing $T = 21 - 17 e^{- t /40}$ .
d
In een afgesloten kamer met één bewoner stijgt de CO₂-concentratie met 760 ppm per uur (hoofdstuk 14 van I). Wat verandert er als je verlucht met een debiet Q?
Zonder verluchting: $d c / d t = 760 ppm/h$ , een constante: geen evenwicht, de concentratie stijgt lineair. Met verluchting: $V d c / d t = G - Q (c - c_{b})$ , lineair, met evenwicht $c_{b} + G / Q$ en $τ = V / Q$ (opdracht 4.7).

4.3

Een richtingsveld lezen en schetsen

Oefenen6Tekenen15 min

De figuur toont het richtingsveld van

\frac{d y}{d t} = t - y

Het richtingsveld van dy/dt = t − y.

a
Waar zijn de lijnstukjes horizontaal? Teken die lijn in de figuur.
Waar t − y = 0: op de rechte y = t (een isokline met helling 0). Boven die rechte dalen de oplossingen, eronder stijgen ze.
b
Schets de oplossingen door (0; 0), (0; 2) en (0; 3).
Zie de figuur hieronder. De oplossing door (0; 2) daalt eerst en heeft een minimum waar ze de rechte y = t kruist, in (1,10; 1,10); die door (0; 3) ook, iets verder; die door (0; 0) vertrekt horizontaal en stijgt.
c
Eén oplossing is een rechte. Welke? Controleer door in te vullen.
$y = t - 1$ : $\overset{y}{˙} = 1$ en $t - y = t - (t - 1) = 1$ . Klopt, voor elke $t$ (door (0; −1)). De algemene oplossing is $y = t - 1 + C e^{- t}$ .
d
Wat gebeurt er met alle oplossingen voor grote t?
$C e^{- t} \to 0$ : alle oplossingen naderen de rechte $y = t - 1$ , ongeacht de beginwaarde.

4.4

Scheiden van veranderlijken

Oefenen6Rekenen20 min

Los op en beantwoord de vraag. Geef eenheden.

a
$\overset{y}{˙} = - 0, 2 y$ ( $t$ in s), $y (0) = 50$ . Bereken de halveringstijd en $y (10)$ .
$y = 50 e^{- 0, 2 t}$ ; $t_{1/2} = ln 2/0, 2 = 3, 47 s$ ; $y (10) = 50 e^{- 2} = 6, 77$ .
b
$\dot{N} = 0, 35 N$ ( $t$ in h), $N (0) = 1000$ bacteriën. Bereken de verdubbelingstijd en $N$ na 24 h.
$N = 1000 e^{0, 35 t}$ ; $t_{2} = ln 2/0, 35 = 1, 98 h$ ; $N (24) = 1000 e^{8, 4} = 4, 45 \cdot 1 0^{6}$ .
c
$\dot{T} = - (T - 20) /15$ ( $t$ in min), $T (0) = 80$ °C. Wanneer is $T = 30$ °C?
$T = 20 + 60 e^{- t /15}$ ; $30 = 20 + 60 e^{- t /15}$ geeft $e^{- t /15} = 1/6$ , $t = 15 ln 6 = 26, 9 min$ .
d
$\overset{y}{˙} = y^{2}$ met $y (0) = 1$ . Wat is er vreemd aan de oplossing?
$d y / y^{2} = d t$ geeft $- 1/ y = t + c$ , met $y (0) = 1$ : $c = - 1$ , dus $y = \frac{1}{1 - t}$ . De oplossing ontploft na een eindige tijd: $y \to \infty$ voor $t \to 1$ . Een niet-lineaire vergelijking kan dat; een lineaire niet.

4.5

Euler met de hand

Oefenen6Rekenen20 min

Gegeven

\overset{y}{˙} = y - t

met

y (0) = 2

. De exacte oplossing is

y = t + 1 + e^{t}

a
Zet drie stappen Euler met h = 0,5 tot t = 1,5. Noteer telkens de helling.
Tabel: t = 0,0: y = 2,000 (exact 2,000, fout 0,000); t = 0,5: y = 3,000 (exact 3,149, fout −0,149); t = 1,0: y = 4,250 (exact 4,718, fout −0,468); t = 1,5: y = 5,875 (exact 6,982, fout −1,107). Hellingen: 2; 2,5; 3,25.
b
Bereken y(0,5) opnieuw met twee stappen van 0,25. Hoeveel keer kleiner is de fout?
$y (0, 25) = 2 + 0, 25 \cdot 2 = 2, 5$ ; helling $2, 5 - 0, 25 = 2, 25$ ; $y (0, 5) = 3, 0625$ . Fout $- 0, 086$ tegenover $- 0, 149$ : $1, 72$ keer kleiner, bijna 2 (orde 1).
c
Waarom ligt Euler hier onder de echte oplossing?
De oplossing is hol ( $\overset{y}{¨} = \overset{y}{˙} - 1 = e^{t} > 0$ ): de helling neemt toe tijdens elke stap, en Euler gebruikt de helling van het begin. Dus blijft Euler achter.

Toepassen

meerdere stappen, in een context

4.6

De klimaatmodule in een lentenacht

Toepassen7Rekenen25 min

Om 22 uur gaan de lampen van de klimaatmodule uit; het is dan 26 °C binnen en 4 °C buiten. Met de nachtbuffer (hoofdstuk 18 van I) is C ≈ 1,1 MJ/K; het warmteverlies is 20,2 W/K, met nachtscherm 16,3 W/K (hoofdstuk 19 van I). De sla mag niet onder 16 °C.

a
Stel de differentiaalvergelijking op en bereken τ.
$C \dot{T} = - U A (T - 4)$ , dus $\dot{T} = - (T - 4) / τ$ met $τ = C / U A = 1, 1 \cdot 1 0^{6} /20, 2 = 54455 s = 15, 1 h$ .
b
Bereken de temperatuur om 6 uur, zonder en met nachtscherm.
$T (8 h) = 4 + 22 e^{- 8/15, 13} = 17, 0$ °C. Met scherm: $τ = 18, 7 h$ , $T = 18, 4$ °C. Beide boven 16 °C.
c
Simuleer zonder nachtscherm met Euler met stappen van 2 h. Hoe groot is de fout om 6 uur? Is die stap stabiel?
Per stap blijft $1 - 2/15, 13 = 0, 8678$ van het verschil over: $4 + 22 \cdot 0, 867 8^{4} = 16, 5$ °C, $0, 5$ °C te koud. Stabiel: $h = 2 h$ is veel kleiner dan $2 τ \approx 30 h$ .
d
Hoe lang mag de nacht (zonder scherm) duren voor de sla onder 16 °C zakt?
$16 = 4 + 22 e^{- t / τ}$ geeft $t = τ ln (22/12) = 9, 2 h$ : tot ongeveer 7:10.

4.7

CO₂ in de woonmodule

Toepassen7Rekenen30 min

De woonmodule is 6,0 × 4,0 × 2,6 m. Twee bewoners ademen elk ongeveer 1,0 kg CO₂ per dag uit (hoofdstuk 14 van I). CO₂ heeft bij 20 °C een dichtheid van 1,83 kg/m³. Buiten is er 420 ppm. De ventilatie ververst Q m³ lucht per uur. Waarschuwing bij 1000 ppm.

a
Bereken de CO₂-productie G in m³/h en stel een vergelijking op voor de concentratie c (als volumefractie).
$G = 2 \cdot 1, 0 kg /24 h /1, 83 kg/ m^{3} = 0, 0455 m^{3} /h$ . Balans voor het CO₂-volume in de module: $V d c / d t = G - Q (c - c_{b})$ met $V = 62, 4 m^{3}$ .
b
Zonder ventilatie: hoeveel ppm per uur stijgt c, en na hoeveel minuten wordt 1000 ppm bereikt vanaf 420 ppm?
$d c / d t = G / V = 7, 30 \cdot 1 0^{- 4}$ per uur $= 730 ppm/h$ . Van 420 naar 1000 ppm: $580/730 = 0, 79 h$ , ongeveer 48 minuten.
c
Met Q = 50 m³/h: bereken het evenwicht, τ en de tijd tot 1000 ppm vanaf 420 ppm.
Evenwicht $c_{b} + G / Q = 420 + 911 = 1331 ppm$ ; $τ = V / Q = 1, 25 h$ . $1000 = 1331 - 911 e^{- t / τ}$ geeft $t = τ ln \frac{911}{331} = 1, 26 h$ .
d
Welk debiet is minstens nodig om onder 1000 ppm te blijven? Hoeveel keer per uur is dat de lucht van de module?
Evenwicht onder 1000 ppm: $G / Q \leq 580 \cdot 1 0^{- 6}$ , dus $Q \geq 79 m^{3} /h$ : $1, 26$ keer per uur het volume van de module. Met warmteterugwinning (basisfiche) kost dat weinig warmte.

4.8

τ uit een afkoelcurve

Toepassen7Analyseren25 min

Na een stroomuitval logt het station om het half uur de binnentemperatuur. Buiten meet de sensor −8,0 °C. De snelle warmtecapaciteit van het station is 540 kJ/K.

$t$ (h)	0,0	0,5	1,0	1,5	2,0	2,5	3,0	3,5	4,0
$T$ (°C)	20,0	15,6	12,1	9,1	6,4	4,2	2,3	0,9	−0,8

a
Bereken $ln (T + 8)$ voor elke meting en toon dat die punten ongeveer op een rechte liggen.
ln(T + 8): 3,332; 3,161; 3,001; 2,839; 2,667; 2,501; 2,332; 2,186; 1,974. Ze dalen gelijkmatig, ongeveer 0,17 per half uur: een rechte.
b
Bepaal τ uit de helling (met een rekenblad of Python), en ook uit enkel de eerste en de laatste meting.
Kleinste kwadraten: helling $- 0, 3344$ per uur, dus $τ = 2, 99 h$ . Uit twee punten: $τ = 4/ ln \frac{28 , 0}{7 , 2} = 2, 95 h$ .
c
Hoe groot is het warmteverlies UA volgens deze meting?
$U A = C / τ = 540 000/ (2, 99 \cdot 3600) = 50 W/K$ , zoals het dossier (50 W/K).
d
Stel dat het buiten eigenlijk −7,0 °C was. Welke τ vind je dan met ln(T + 7)? Wat leer je daaruit?
Met $- 7$ °C: $τ = 2, 77 h$ , en de punten liggen niet meer op een rechte: $T + 7$ gaat niet naar nul maar naar −1, dus de punten buigen naar beneden af. Een fout van 1 °C in het evenwicht verandert $τ$ met −7 %: meet de buitentemperatuur even zorgvuldig als de binnentemperatuur.

4.9

Euler in Python op de algen

Toepassen7Programmeren25 min

De algen in het testreservoir volgen

\dot{N} = r N (1 - N / K)

met

r = 0, 70

per dag,

K = 7, 7 \cdot 1 0^{5}

per mL en

N (0) = 1900

per mL.

Schrijf een functie euler(f, y0, t0, t1, h) en simuleer tot dag 10 met h = 1; 0,5; 0,25; 0,1 en 0,01 dag.

Bijvoorbeeld:

import math

r, K, N0 = 0.70, 7.7e5, 1900.0                       # per dag, cellen per mL


def f(N, t):
    return r * N * (1 - N / K)


def euler(f, y0, t0, t1, h):
    n = round((t1 - t0) / h)
    t, y = t0, y0
    for k in range(n):
        y = y + h * f(y, t)
        t = t + h
    return y


exact = K / (1 + (K - N0) / N0 * math.exp(-r * 10))
print(f"exact op dag 10: {exact:.4g}")
for h in (1.0, 0.5, 0.25, 0.1, 0.01):
    N = euler(f, N0, 0, 10, h)
    print(f"h = {h:5.2f} dag: N(10) = {N:.4g}   relatieve fout {N / exact - 1:+.3f}")

Vergelijk met de exacte oplossing op dag 10. Bevestigt de fout dat Euler van orde 1 is?

Echte uitvoer:

exact op dag 10: 5.626e+05
h =  1.00 dag: N(10) = 2.919e+05   relatieve fout -0.481
h =  0.50 dag: N(10) = 4.266e+05   relatieve fout -0.242
h =  0.25 dag: N(10) = 4.986e+05   relatieve fout -0.114
h =  0.10 dag: N(10) = 5.384e+05   relatieve fout -0.043
h =  0.01 dag: N(10) = 5.603e+05   relatieve fout -0.004

Voor kleine h halveert de relatieve fout ongeveer bij elke halvering van h (orde 1). Bij h = 1 dag is de fout veel groter dan lineair zou voorspellen: daar is de stap niet meer klein.

c
Waarom is een stap van 1 dag hier zo slecht, terwijl de algen maar één keer per dag verdubbelen?
Per stap vermenigvuldigt Euler in het begin met $1 + r h = 1, 70$ in plaats van $e^{r h} = 2, 01$ : 15 % te weinig per dag, en dat stapelt op (samengestelde interest). Een stap is pas klein als $r h ≪ 1$ .

4.10

Het filter op de batterijmeting

Toepassen7Rekenen25 min

De batterijspanning van de module komt via een deler (100 kΩ boven, 33 kΩ onder) op A3, met 100 nF over de onderste weerstand (hoofdstuk 22 van I). Als de lampen aangaan, zakt de batterij van 13,30 V naar 12,80 V.

a
Stel met de knooppuntwet een vergelijking op voor de spanning u over de condensator. Toon dat τ = (R₁ ∥ R₂) · C, en bereken τ.
Knooppunt boven C: stroom door $R_{1}$ erin = stroom door $R_{2}$ + stroom in C: $\frac{U _{bat} - u}{R _{1}} = \frac{u}{R _{2}} + C \overset{u}{˙}$ . Herschikken: $C \overset{u}{˙} = \frac{U _{th} - u}{R _{th}}$ met $U_{th} = U_{bat} \frac{R _{2}}{R _{1} + R _{2}}$ en $R_{th} = R_{1} ∥ R_{2} = 24, 8 kΩ$ . $τ = 2, 48 ms$ .
b
Hoe lang duurt het tot u minder dan één ADC-stap (4,88 mV) van de nieuwe waarde af ligt?
Van $3, 300$ naar $3, 176 V$ : een sprong van $124, 1 mV$ . $124, 1 e^{- t / τ} = 4, 88$ geeft $t = τ ln \frac{124 , 1}{4 , 88} = 8, 0 ms$ .
c
Simuleer de sprong met Euler met stappen van 1 ms en van 5 ms (5 stappen volstaan). Wat gebeurt er?
Met 1 ms: factor $1 - 1/2, 48 = 0, 597$ per stap, de simulatie daalt netjes: $3, 300; 3, 250; 3, 220; 3, 202; 3, 192; 3, 185 V$ . Met 5 ms: $h / τ = 2, 015 > 2$ , factor $- 1, 015$ : de simulatie springt om beurten boven en onder $3, 176 V$ en de sprongen worden groter ( $3, 300; 3, 050; 3, 304; 3, 046; 3, 308; 3, 042 V$ ). Onstabiel, al is de stap maar 5 ms.

4.11

Habitat: een warmtebuffer voor de meetkas

Toepassen7Habitat60 min

Het team wil eis K7 halen: 's nachts minstens 1 °C warmer dan buiten. Gebruik kas.py uit het handboek (C = 30 kJ/K, a = 0,30 m², UA = 16,7 W/K) en simuleer drie keer dezelfde aprildag na elkaar. Water: 4,19 kJ/K per liter. Het energiebudget van de meetkas is 44 tot 63 Wh per dag.

a
Pas kas.py aan zodat het drie dagen simuleert en P_verw en extra water (in liter) als parameter heeft.
Bijvoorbeeld een parameter dagen=3 (de tijd modulo 24 voor zon en buitentemperatuur), P_verw als functie van de klok en C = 30e3 + liter * 4186; evalueer de laatste 1441 punten.

Vul de tabel in voor: niets, 17 W verwarming van 21 tot 6 uur, 20 L water, 40 L water.

Model:

variant (model)	hoogste T overdag	22-6 uur: kleinste T − T_b	22-6 uur: gemiddeld T − T_b
zonder buffer	30,5 °C	0,1 °C	0,3 °C
17 W verwarming 's nachts	30,5 °C	1,1 °C	1,3 °C
+ 20 L water	28,8 °C	0,5 °C	1,5 °C
+ 40 L water	26,6 °C	1,4 °C	3,4 °C

Energie: 0 Wh, 150 Wh, 0 Wh, 0 Wh per nacht.

c
Welke oplossing kies je? Weeg K7, het energiebudget, de hoogste temperatuur en de praktische inpassing af.
Verwarmen haalt K7 maar kost ongeveer 150 Wh per nacht: drie keer het dagbudget, onhaalbaar met 50 Wp. 20 L haalt K7 volgens het model niet (kleinste verschil 0,5 °C); 40 L wel (1,4 °C) en verlaagt bovendien de top van 30,5 naar 26,6 °C. Keuze: ongeveer 40 L in donkere flessen langs de achterwand, in de zon (let op het gewicht: 40 kg extra, wel gunstig tegen kantelen).
d
Welke twee beperkingen van het model kunnen je conclusie veranderen, en hoe test je ze?
Eén knoop: water en lucht hebben in het model dezelfde temperatuur, in werkelijkheid geeft een fles zijn warmte traag af (tweede knoop, hoofdstuk 17), dus het echte effect is kleiner. En de trage massa van grond en reservoir ontbreekt al in het basismodel (dat 's nachts te koud is). Test: flessen plaatsen en een week loggen, met en zonder, en het model met twee knopen fitten.

Zo wordt dit beoordeeld

Simulatie	Drie dagen, laatste dag geëvalueerd; P_verw alleen 's nachts; water als extra C.
Resultaten	Tabel volledig en juist, met eenheden.
Afweging	K7, energie (17 W · 9 h ≈ 150 Wh tegenover ±50 Wh per dag), temperatuur overdag en plaats/gewicht afgewogen.
Kritiek	Eén knoop als bovengrens; koppeling water-lucht en trage massa; een proef met de logger voorgesteld.

Uitdagen

transfer, open problemen, leerlijnen combineren

4.12

De logistische oplossing afleiden

Uitdagen8Analyseren30 min

Het handboek gaf de oplossing van

\dot{N} = r N (1 - N / K)

met een schets van de afleiding. Maak ze volledig.

a
Toon dat $\frac{K}{N ( K - N )} = \frac{1}{N} + \frac{1}{K - N}$ en leid daarmee $N (t)$ af voor $0 < N_{0} < K$ .
Samen: $\frac{( K - N ) + N}{N ( K - N )} = \frac{K}{N ( K - N )}$ . Scheiden: $\frac{K d N}{N ( K - N )} = r d t$ , integreren: $ln N - ln (K - N) = r t + c$ (voor $0 < N < K$ ). Dus $\frac{N}{K - N} = A e^{r t}$ met $A = \frac{N _{0}}{K - N _{0}}$ . Oplossen: $N = A e^{r t} (K - N)$ , $N (1 + A e^{r t}) = A K e^{r t}$ , $N = \frac{K}{1 + e ^{- r t} / A} = \frac{K}{1 + \frac{K - N _{0}}{N _{0}} e ^{- r t}}$ .
b
Toon dat $N (t)$ een buigpunt heeft waar $N = K /2$ , op $t^{*} = \frac{1}{r} ln \frac{K - N _{0}}{N _{0}}$ . Bereken $t^{*}$ en $\dot{N} (t^{*})$ voor de algen ( $r = 0, 704$ per dag, $K = 7, 68 \cdot 1 0^{5}$ , $N_{0} = 1900$ ).
$\ddot{N} = \frac{d}{d t} f (N) = f^{'} (N) \dot{N} = r (1 - \frac{2 N}{K}) \dot{N}$ (kettingregel). Nul (en van teken wisselend) waar $N = K /2$ : een buigpunt. $N = K /2$ invullen: $1 + \frac{K - N _{0}}{N _{0}} e^{- r t} = 2$ , dus $t^{*} = \frac{1}{r} ln \frac{K - N _{0}}{N _{0}} = 8, 52$ dagen. $\dot{N} (t^{*}) = r K /4 = 1, 35 \cdot 1 0^{5}$ per mL per dag.
c
Wat verandert er in de afleiding als N₀ > K? Is de formule dan nog geldig?
Voor $N > K$ is $K - N < 0$ : gebruik $ln ∣ K - N ∣$ . Dan is $A = \frac{N _{0}}{K - N _{0}} < 0$ en dezelfde formule geldt: de noemer $1 + \frac{K - N _{0}}{N _{0}} e^{- r t}$ ligt tussen 0 en 1, dus $N > K$ en dalend naar $K$ (de bovenste kromme in het richtingsveld).

4.13

De fout van Euler, exact

Uitdagen9Analyseren30 min

Voor

\overset{y}{˙} = - y / τ

kan je de Euler-benadering exact opschrijven:

y_{n} = (1 - h / τ)^{n} y_{0}

t = nh

a
Toon met $ln (1 - x) = - x - \frac{1}{2} x^{2} - \dots$ dat $\frac{y _{n}}{y ( t )} \approx 1 - \frac{t h}{2 τ ^{2}}$ voor kleine $h / τ$ .
$ln y_{n} = ln y_{0} + n ln (1 - h / τ) = ln y_{0} - \frac{nh}{τ} - \frac{n h ^{2}}{2 τ ^{2}} - \dots = ln y (t) - \frac{t h}{2 τ ^{2}} - \dots$ Dus $\frac{y _{n}}{y ( t )} = e^{- t h / (2 τ^{2})} \approx 1 - \frac{t h}{2 τ ^{2}}$ : de relatieve fout is evenredig met $h$ (orde 1) en groeit met $t$ .
b
Leid daaruit de fout van het station af (het verschil met buiten is $28 e^{- t / τ}$ ) op $t = 3 h$ , en vergelijk met de kolom fout/h van euler.py.
Absolute fout: $- \frac{t h}{2 τ ^{2}} \cdot 28 e^{- t / τ}$ . Op $t = τ = 3$ : $- \frac{3 h}{18} \cdot 28 e^{- 1} = - 1, 717 h$ . euler.py toont fout/h $= - 1, 719$ bij $h = 0, 01$ : klopt.
c
Voor welke t is de absolute fout het grootst (bij vaste h)? Interpreteer.
$E (t) \propto t e^{- t / τ}$ is maximaal waar $\frac{d}{d t} (t e^{- t / τ}) = e^{- t / τ} (1 - t / τ) = 0$ : op $t = τ$ . Eerst groeien de fouten aan (ze stapelen), later dooft alles uit naar het evenwicht, ook de fout.
d
Wat zegt $(1 - h / τ)^{n}$ als $h > 2 τ$ ?
$∣1 - h / τ ∣ > 1$ : $∣ y_{n} ∣$ groeit geometrisch en wisselt van teken. De reeksontwikkeling geldt dan niet meer: de fout is niet klein maar onbegrensd (instabiliteit).

4.14

Onderzoek: een stijf stelsel

Uitdagen9Onderzoeken45 min

Het station (

τ = 3 h

K = 0, 020

°C/W) wordt verwarmd met constant 1000 W; zijn sensor volgt met

τ_{s} = 60 s

(hoofdstuk 39). Met

y

de stijging boven buiten en

y_{m}

de meting:

τ \overset{y}{˙} = - y + K u

τ_{s} \overset{y}{˙}_{m} = y - y_{m}

Simuleer twee uur met Euler voor dt = 10, 60, 110 en 130 s en druk y en y_m af.

Bijvoorbeeld (beide toestanden in één toewijzing, zodat de oude waarden gebruikt worden):

TAU, K, TAU_S = 10800.0, 0.02, 60.0     # station (s), versterking (°C/W), sensor (s)
u = 1000.0                                # constante verwarming (W)

for dt in (10.0, 60.0, 110.0, 130.0):
    y, ym, t = 0.0, 0.0, 0.0                  # stijging boven buiten: echt en gemeten (°C)
    while t < 7200.0 - 1e-9:                  # twee uur
        y, ym = y + dt * (-y + K * u) / TAU, ym + dt * (y - ym) / TAU_S
        t += dt
    print(f"dt = {dt:5.0f} s: y = {y:6.3f} °C, gemeten {ym:10.3f} °C")

Echte uitvoer:

dt =    10 s: y =  9.735 °C, gemeten      9.677 °C
dt =    60 s: y =  9.751 °C, gemeten      9.693 °C
dt =   110 s: y =  9.824 °C, gemeten      9.767 °C
dt =   130 s: y =  9.849 °C, gemeten    636.691 °C

b
Verklaar het resultaat met de stabiliteitsgrens. Welke tijdconstante beslist, en waarom?
Het station zelf is na 2 uur 9,7 °C gestegen (exact 20 · (1 − e^(−2/3)) = 9,73 °C): voor y zijn alle stappen ruim klein genoeg tegenover τ = 3 h. De sensor volgt netjes tot dt = 110 s; bij 130 s is h/τ_s = 2,17 > 2 en ontploft y_m. De snelste tijdconstante (60 s) beslist, ook al verandert de meting na een minuut niet meer wezenlijk: dat is een stijf stelsel.
c
De impliciete Euler gebruikt de helling aan het eind van de stap: $y_{n + 1} = y_{n} + h f (y_{n + 1})$ . Toon dat die voor $\overset{y}{˙} = - y / τ$ geeft $y_{n + 1} = y_{n} / (1 + h / τ)$ , en dat ze voor elke $h > 0$ stabiel is.
$y_{n + 1} = y_{n} - \frac{h}{τ} y_{n + 1}$ geeft $y_{n + 1} (1 + h / τ) = y_{n}$ . De factor $\frac{1}{1 + h / τ}$ ligt tussen 0 en 1 voor elke $h > 0$ : altijd dalend, nooit ontploffend. Voor een lineair stelsel moet je per stap wel een stelsel oplossen (hoofdstuk 5, en hoofdstuk 36).

Zo wordt dit beoordeeld

Code	Correct stelsel, beide toestanden tegelijk bijgewerkt (met de oude waarden).
Analyse	Grens 2τ_s = 120 s herkend; de snelste tijdconstante beslist, al is ze oninteressant.
Impliciet	Afleiding correct; \|1/(1 + h/τ)\| < 1 voor elke h > 0.

4.15

Open modelleervraag: algen oogsten

Uitdagen8Onderzoeken50 min

In hoofdstuk 31 kweekt het station algen voor zuurstof. Stel dat ze logistisch groeien zoals in het testreservoir (r = 0,70 per dag, K = 7,7 · 10⁵ per mL) in een vat van 60 L, en dat je ze oogst.

a
Continu oogsten met een fractie $E$ per dag: $\dot{N} = r N (1 - N / K) - E N$ . Zoek de evenwichten en de voorwaarde op $E$ opdat de kweek niet uitsterft.
$\dot{N} = N (r - E - r N / K) = 0$ : $N = 0$ of $N^{*} = K (1 - E / r)$ . Positief alleen als $E < r$ ; dan is $N^{*}$ stabiel en $N = 0$ onstabiel. Bij $E \geq r$ sterft de kweek uit.
b
De oogst per dag is E · N*. Welke E geeft de grootste oogst (hoofdstuk 2)? Hoeveel cellen per dag is dat voor het hele vat?
$Y (E) = E K (1 - E / r)$ : $Y^{'} (E) = K (1 - 2 E / r) = 0$ bij $E = r /2 = 0, 35$ per dag, $N^{*} = K /2$ , $Y = r K /4 = 1, 35 \cdot 1 0^{5}$ per mL per dag. Voor 60 L (60 000 mL): $8, 1 \cdot 1 0^{9}$ cellen per dag.
c
Eén keer per dag oogsten is praktischer: elke dag een fractie f weghalen, daartussen logistisch groeien. Welke f is best (simuleer)? Vergelijk met continu.
Simulatie (exacte logistische groei over één dag, dan oogsten, 60 dagen tot een vaste cyclus): de beste fractie is ongeveer $f = 0, 30$ , met $1, 34 \cdot 1 0^{5}$ cellen per mL per dag, 99 % van het continue optimum. Dagelijks oogsten kost dus maar weinig.
d
Waarom zegt de K van het testreservoir weinig over een echte algenreactor? Welke grootheden zou je eerst meten?
In het testreservoir begrensde het licht de algen (ze beschaduwen elkaar in een diep, weinig belicht vat). Een reactor is ontworpen voor licht (dunne lagen, lampen), CO₂ en voeding: K en r kunnen tien tot honderd keer groter zijn. Meet in de reactor zelf de groei van een kleine beginwaarde (r) en de eindwaarde zonder oogst (K), bij het gekozen licht.

Zo wordt dit beoordeeld

Model	Evenwichten 0 en K(1 − E/r); voorwaarde E < r; stabiliteit besproken.
Optimum	E = r/2, N* = K/2, oogst rK/4, per vat omgerekend.
Simulatie	Dagelijkse oogst correct gemodelleerd; optimum gevonden en vergeleken.
Kritiek	K hangt af van licht (en voeding, CO₂); r en K meten in de reactor zelf.

4.16

Programmeeronderzoek: hoe klein moet de stap zijn?

Uitdagen9Programmeren45 min

Je weet niet altijd wat de exacte oplossing is. Dan schat je de numerieke fout door de stap te halveren. Gebruik kas.py (C = 30 kJ/K, a = 0,30 m²) en bekijk de temperatuur om 14:00.

a
Bereken T(14:00) met stappen van 60, 30, 15, 7,5, 3,75, ... minuten (zeven stappen). Maak een tabel.
Een lus over dt = 3600 / 2**i met simuleer tot 14 h. Resultaat:
stap (min) 60 30 15 7,5 3,75 1,875 0,9375
$T$ om 14:00 (°C) 32,11 30,69 30,62 30,58 30,56 30,55 30,55
b
Toon dat het verschil tussen twee opeenvolgende resultaten telkens ongeveer halveert. Wat zegt dat?
Verschillen: −1,418, −0,070, −0,037, −0,018, −0,009, −0,005 °C: telkens ongeveer de helft. De fout is evenredig met de stap: orde 1, zoals voorspeld. (Bij 60 min is h net iets groter dan 2τ: de simulatie is al zwak onstabiel en 1,4 °C te warm; daar klopt de verhouding niet.)
c
Bij orde 1 is $T (h) \approx T^{*} + c h$ . Leid af dat $T^{*} \approx 2 T (h /2) - T (h)$ (extrapolatie van Richardson) en schat zo de exacte waarde.
$T (h) = T^{*} + c h$ en $T (h /2) = T^{*} + c h /2$ ; aftrekken: $c h /2 = T (h) - T (h /2)$ , dus $T^{*} = T (h /2) - c h /2 = 2 T (h /2) - T (h)$ . Met de twee kleinste stappen: $T^{*} \approx 30, 543$ °C (met een stap van 5 s: $30, 543$ °C).
d
Welke stap volstaat voor een fout onder 0,05 °C? Waarom rekent de Pi toch met 1 minuut?
De fout bij stap h is ongeveer T(h) − T*: bij 15 min 0,07 °C, bij 7,5 min 0,04 °C. Een stap van 7,5 minuten volstaat dus. De Pi rekent met 1 minuut omdat het niets kost (1440 stappen per dag) en een ruime marge geeft tegenover de stabiliteitsgrens (2τ ≈ 60 min), ook als C kleiner blijkt.

Zo wordt dit beoordeeld

Code	Stappen correct gehalveerd, telkens tot 14:00 gesimuleerd.
Orde	Halverende verschillen herkend als orde 1.
Richardson	Afleiding en toepassing correct.
Keuze	Stap met onderbouwde marge; stabiliteit (τ ≈ 30 min) mee in rekening gebracht.