EExpert · deel E1 · Wiskundige gereedschapskist

Handboek · hoofdstuk 4Differentiaalvergelijkingen en simulatie

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

Bijna elk natuurkundig of biologisch model zegt niet wat een grootheid is, maar hoe snel ze verandert: de kas warmt op met de zon en verliest warmte naar buiten, algen groeien tot hun licht op is, een condensator laadt tot de bron. Zo'n model is een differentiaalvergelijking. Je leert ze opstellen als balans, lezen in een richtingsveld, oplossen door veranderlijken te scheiden, en, als dat niet lukt, simuleren met de methode van Euler, met oog voor stapgrootte, fout en stabiliteit. Het hoofdstuk eindigt met de eerste simulatie van de meetkas, naast haar echte uurwaarden.

8× uitleg6× uitgewerkt voorbeeld1× ontleed1× naslag1× verhaal1× het geheel3× code16 opdrachten in het werkboek± 12 lestijden

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.

Na dit hoofdstuk

Je kan een veranderingsproces vertalen in een differentiaalvergelijking dy/dt = f(y, t) met een beginwaarde, als balans van wat erin komt en eruit gaat, en de eenheden en evenwichten ervan interpreteren.
Je kan een richtingsveld lezen en schetsen en er oplossingskrommen in tekenen.
Je kan een oplossing controleren door ze in te vullen, en exponentiële groei en verval, de afkoelingswet van Newton, de RC-kring en de logistische vergelijking oplossen door scheiden van veranderlijken.
Je kan tijdconstante, halverings- en verdubbelingstijd, groeisnelheid en draagkracht berekenen en uit metingen schatten.
Je kan de methode van Euler afleiden uit de raaklijn en toepassen, met de hand en in Python, ook op een stelsel.
Je kan de fout van Euler tegenover de stapgrootte verklaren (orde 1) en de stabiliteitsgrens h < 2τ afleiden en gebruiken.
Je kan een simulatie vergelijken met metingen en uit de residuen afleiden of een parameter slecht gekozen is of een mechanisme ontbreekt.
Je kan een eerste warmtemodel van de meetkas simuleren en er een onderbouwde ontwerpbeslissing mee nemen.

Verhaal · 4.1

Rekenen in stapjes: van Euler tot het weerbericht

1/21

In 1768 verscheen in Sint-Petersburg het eerste deel van de Institutiones calculi integralis, het leerboek integraalrekening van Leonhard Euler. Euler wist dat de meeste vergelijkingen die zeggen hoe iets verandert, geen oplossing in een formule hebben. Hij beschreef een eenvoudige uitweg: zet een klein stapje langs de raaklijn, reken op de nieuwe plaats de nieuwe helling uit, zet het volgende stapje, en zo verder. Hoe kleiner de stapjes, hoe dichter de gebroken lijn bij de echte oplossing ligt. In de jaren 1820 bewees Augustin-Louis Cauchy in zijn lessen aan de École polytechnique in Parijs dat dat klopt als de stap naar nul gaat. Die gebroken lijn heet vandaag de methode van Euler.

Rekenen in stapjes was toen al ernst. In 1757 en 1758 rekenden de wiskundige Alexis Clairaut, de astronoom Joseph Jérôme Lalande en de rekenaarster Nicole-Reine Lepaute maandenlang met de hand uit hoe Jupiter en Saturnus de komeet van Halley zouden vertragen. In november 1758 kondigden ze aan dat de komeet rond half april 1759 het dichtst bij de zon zou komen, met een marge van een maand. Het werd 13 maart: net binnen de marge, en een triomf voor de mechanica van Newton.

Tijdens de Eerste Wereldoorlog probeerde de Britse natuurkundige Lewis Fry Richardson, als vrijwilliger bij een ambulancedienst in Frankrijk, het weer van zes uur verder uit te rekenen met de vergelijkingen van de atmosfeer, opgedeeld in cellen en tijdstappen. Na weken rekenwerk kwam er voor één punt een drukstijging van 145 hPa in zes uur uit: onmogelijk. Toch publiceerde hij zijn methode in 1922, met de droom van een zaal met 64 000 rekenaars die het weer sneller zouden uitrekenen dan het zelf veranderde. Veel later bleek uit een reconstructie dat niet zijn methode de schuldige was, maar de begintoestand: kleine onevenwichten in de gemeten wind en druk. In 1928 toonden Richard Courant, Kurt Friedrichs en Hans Lewy aan dat een tijdstap ook te groot kan zijn: dan explodeert de berekening, hoe goed de vergelijkingen ook zijn. In 1950 lukte de eerste weersvoorspelling met een computer, op de ENIAC, door Jule Charney, Ragnar Fjørtoft en John von Neumann.

Beginwaarde, stapgrootte, stabiliteit: de drie struikelblokken uit dit verhaal zijn precies de onderwerpen van dit hoofdstuk.

Waarom dit hoofdstuk ertoe doet

Het brein van het station moet vooruitkijken: hoe warm wordt de kas vanmiddag, hoe lang houdt de batterij het, wanneer wordt het reservoir groen? Dat kan alleen met een model dat zegt hoe snel iets verandert, en een computer die dat stap voor stap uitrekent. De digitale tweeling van hoofdstuk 49 is een groot stelsel zulke vergelijkingen, opgelost met de methoden van dit hoofdstuk en hoofdstuk 36.

Clairaut, Lalande en Lepaute rekenden de baan van de komeet van Halley stap voor stap met de hand uit.

Uitleg · 4.2

Wat je al weet: veranderingen die van de toestand afhangen

2/21

wat je al kan	wat nu komt
groei en verval met een groeifactor per tijdstap, verdubbelings- en halveringstijd (hoofdstuk 5 van I)	dezelfde processen met een verandering per ogenblik: $d y / d t = k y$
RC-kring: per stapje een vaste fractie van wat ontbreekt; $τ = R C$ (hoofdstuk 24 van I)	de RC-vergelijking uit de wetten van Kirchhoff, exact en stap voor stap
de afgeleide als helling, de raaklijn als benadering (hoofdstuk 1)	een vergelijking die de helling geeft; stappen langs de raaklijn (Euler)
integreren als ophopen, de linker Riemann-som (hoofdstuk 3)	een simulatie telt de veranderingen op, stap na stap
warmtecapaciteit $C$ en warmteverlies $U A \cdot Δ T$ (hoofdstuk 18 en 19 van I)	de warmtebalans als differentiaalvergelijking: het eerste model van de meetkas

In al die voorbeelden hangt de verandering af van de toestand zelf: hoe meer algen, hoe meer er bijkomen; hoe warmer de kas, hoe meer warmte ze verliest. Zo'n verband tussen een grootheid en haar afgeleide is een differentiaalvergelijking. Hier gaat het over de eenvoudigste soort, van de eerste orde: alleen de eerste afgeleide komt voor.

beginwaardeprobleem

\frac{d y}{d t} = f (y, t), y (t_{0}) = y_{0}

f

geeft in elk punt

(t, y)

de helling van de onbekende functie

y (t)

. Een differentiaalvergelijking heeft een hele familie oplossingen; de beginwaarde kiest er één uit.

Een model opstellen is bijna altijd een balans: verandering = wat erin komt − wat eruit gaat. Voor de warmte in de meetkas: de zon brengt warmte binnen, een eventuele verwarming ook, en door de wanden lekt er $U A (T - T_{b})$ weg (hoofdstuk 19 van I). Wat overblijft, warmt de kas op met $C d T / d t$ (hoofdstuk 18 van I):

C \frac{d T}{d t} = a G (t) + P - U A (T - T_{b} (t)) .

Controleer altijd de eenheden: links J/K maal K/s = W, rechts drie vermogens. Twee woorden die terugkomen: een vergelijking is lineair als $y$ er alleen in de eerste macht in staat (zoals hierboven), en autonoom als $f$ niet expliciet van $t$ afhangt. De kas is niet autonoom: de zon en de buitentemperatuur veranderen met de klok. Een evenwicht is een waarde $y^{*}$ waar $f (y^{*}) = 0$ : daar staat alles stil.

model	vergelijking	evenwicht	in Habitat
exponentiële groei	$\dot{N} = r N$	$N = 0$ (onstabiel)	algen in de eerste dagen
exponentieel verval	$\overset{y}{˙} = - y / τ$	$y = 0$ (stabiel)	condensator die ontlaadt
afkoelen (Newton)	$\dot{T} = - (T - T_{b}) / τ$	$T = T_{b}$	station na een stroomuitval
RC laden	$R C \overset{u}{˙} = U_{0} - u$	$u = U_{0}$	ontdenderfilter, filter op de batterijmeting
logistische groei	$\dot{N} = r N (1 - N / K)$	$0$ en $K$	algen in het reservoir
warmtebalans	$C \dot{T} = a G (t) + P - U A (T - T_{b} (t))$	verschuift met $t$	meetkas, station

De warmtebalans van de meetkas: wat de zon en een verwarming binnenbrengen, min wat door de wanden weglekt, warmt de kas op.

Uitleg · 4.3

Het richtingsveld: de vergelijking tekent haar eigen oplossingen

3/21

Een differentiaalvergelijking geeft in elk punt $(t, y)$ een helling. Teken in een rooster van punten een kort lijnstukje met die helling: dat is het richtingsveld. Een oplossing is een kromme die overal aan de lijnstukjes raakt, zoals een blad dat met de stroming meedrijft.

Het voorbeeld is het station na een stroomuitval (hoofdstuk 39): buiten is het −8 °C, de tijdconstante is $τ = C / U A = 540 kJ/K /50 W/K = 3 h$ :

\frac{d T}{d t} = - \frac{T + 8}{3} (T in^{\circ} C, t in h) .

De helling hangt alleen van $T$ af (autonoom): op elke horizontale lijn staan de lijnstukjes evenwijdig.
Op $T = - 8$ °C zijn ze horizontaal: het evenwicht. Erboven dalen ze, eronder stijgen ze, en hoe verder van −8 °C, hoe steiler. Alle oplossingen lopen naar −8 °C: een stabiel evenwicht.
De oplossing door $(0; 20)$ is het station dat op 20 °C stond; door $(0; 6)$ een station dat al kouder was; door $(0; - 12)$ een doos die kouder was dan buiten (wiskundig mag dat) en opwarmt naar −8 °C.
Twee oplossingen kruisen elkaar nooit. Als $f$ netjes is (afleidbaar), gaat door elk punt precies één oplossing: de beginwaarde bepaalt alles wat volgt. Dat is de stelling van Picard en Lindelöf, en de reden dat een simulatie een betrouwbare voorspelling kan zijn.

Een differentiaalvergelijking heeft één oplossing

Ze heeft er oneindig veel: één door elk beginpunt. Wie een model opstelt maar de beginwaarde vergeet (hoe warm was het station toen de stroom uitviel?), kan niets voorspellen.

Het richtingsveld van dT/dt = −(T + 8)/3: de lijnstukjes wijzen de weg, de oplossingen volgen ze en lopen allemaal naar −8 °C.

Uitgewerkt voorbeeld · 4.4

Klopt deze oplossing?

4/21

Toon aan dat $T (t) = - 8 + 28 e^{- t /3}$ de oplossing is van $\frac{d T}{d t} = - \frac{T + 8}{3}$ met $T (0) = 20^{\circ} C$ . Geef ook de algemene oplossing en bereken wanneer het station 10 °C bereikt.

Gegeven

$\dot{T} = - (T + 8) /3$ , $T (0) = 20^{\circ} C$
voorstel: $T (t) = - 8 + 28 e^{- t /3}$ ( $t$ in h)

Gevraagd

controle van vergelijking en beginwaarde
de algemene oplossing
het tijdstip waarop T = 10 °C

Oplossing

1
Linkerlid: afleiden met de kettingregel (hoofdstuk 1): $\dot{T} (t) = 28 \cdot (- \frac{1}{3}) e^{- t /3} = - \frac{28}{3} e^{- t /3}$ .
2
Rechterlid: $- \frac{T + 8}{3} = - \frac{28 e ^{- t /3}}{3}$ . Beide zijn gelijk voor elke $t$ : de functie voldoet aan de vergelijking. En $T (0) = - 8 + 28 = 20$ : ook de beginwaarde klopt.
3
Algemeen: $T (t) = - 8 + A e^{- t /3}$ voldoet voor elke constante $A$ (dezelfde rekening). De beginwaarde geeft $A = T (0) + 8$ . Dat is de familie krommen in het richtingsveld.
4
10 °C: $- 8 + 28 e^{- t /3} = 10$ geeft $e^{- t /3} = \frac{18}{28}$ , dus $t = 3 ln \frac{28}{18} = 1, 325 h$ , ongeveer 80 minuten.

Antwoord

T (t) = - 8 + 28 e^{- t /3}

voldoet aan vergelijking en beginwaarde; algemeen

T = - 8 + A e^{- t /3}

. Na ongeveer 80 minuten is het nog 10 °C in het station.

Klopt dit? In het richtingsveld zakt de kromme door

(0; 20)

door 10 °C net na

t = 1, 3 h

. Het verschil met buiten halveert elke

3 ln 2 = 2, 08 h

: van 28 naar 14 °C (station op 6 °C) na 2,1 h. Invullen is de snelste controle van elke oplossing, ook van een simulatie: reken

\dot{T}

numeriek uit en vergelijk met

f

Uitleg · 4.5

Scheiden van veranderlijken: groei en verval

5/21

Waar komt zo'n oplossing vandaan? Neem de eenvoudigste vergelijking, $d y / d t = k y$ : de verandering is evenredig met wat er is. Breng alle $y$ naar links en alle $t$ naar rechts, alsof $d y$ en $d t$ kleine getallen zijn (Leibniz, hoofdstuk 1), en integreer beide kanten (hoofdstuk 3):

\frac{d y}{y} = k d t \Rightarrow \int \frac{d y}{y} = \int k d t \Rightarrow ln ∣ y ∣ = k t + c \Rightarrow y = y_{0} e^{k t} .

Streng genomen deel je door $y$ (dus $y \neq = 0$ ; $y = 0$ is apart een oplossing, het evenwicht) en gebruik je de kettingregel: $\frac{d}{d t} ln ∣ y ∣ = \overset{y}{˙} / y = k$ . De constante $c$ volgt uit de beginwaarde: $e^{c} = ∣ y_{0} ∣$ .

exponentiële groei (k > 0) en verval (k < 0)

\frac{d y}{d t} = k y \Rightarrow y = y_{0} e^{k t} τ = \frac{1}{∣ k ∣} t_{1/2} = τ ln 2 t_{2} = \frac{ln 2}{k}

De tijdconstante

τ

is de tijd waarin een verval tot

1/ e = 37 %

zakt; de halveringstijd is

0, 693 τ

. Verband met de groeifactor per tijdstap

Δ t

van hoofdstuk 5 van I:

g = e^{k Δ t}

, dus

k = ln g /Δ t

. Algen die elke dag verdubbelen, hebben

k = ln 2 = 0, 693

per dag.

Met één substitutie los je meteen alle modellen met een evenwicht op. Staat er $d y / d t = - (y - y_{\infty}) / τ$ , noem dan $u = y - y_{\infty}$ : omdat $y_{\infty}$ constant is, is $\overset{u}{˙} = \overset{y}{˙} = - u / τ$ , en dus $u = u_{0} e^{- t / τ}$ :

naar een evenwicht, met tijdconstante τ

\frac{d y}{d t} = \frac{y _{\infty} - y}{τ} \Rightarrow y (t) = y_{\infty} + (y_{0} - y_{\infty}) e^{- t / τ}

Het afkoelende station (

y_{\infty} = - 8

°C,

τ = 3 h

), de thee, de condensator die laadt (

y_{\infty} = U_{0}

τ = R C

): dezelfde formule. Na

τ

is 63 % van de weg afgelegd, na

3 τ

95 %, na

5 τ

meer dan 99 %.

Scheiden lukt voor elke vergelijking van de vorm $d y / d t = g (t) h (y)$ , als je de twee integralen kan uitrekenen. Dat is minder vaak dan je zou willen; daarom de tweede helft van dit hoofdstuk.

Na t½ = 0,69 τ is nog de helft over, na τ 37 %, na 3τ 5 %; de raaklijn in t = 0 bereikt nul na precies τ (hoofdstuk 1).

Uitgewerkt voorbeeld · 4.6

De afkoelingswet van Newton, getoetst aan een beker thee

6/21

In hoofdstuk 5 van Intermediate koelde een beker thee af in een kamer van 21,0 °C:

$t$ (min)	0	4	8	12	16	20	24	28	32	36	40
$T$ (°C)	85,0	75,0	66,1	59,2	53,0	47,9	43,9	40,1	37,0	34,6	32,5

Stel een differentiaalvergelijking op vanuit een warmtebalans, los ze op, bepaal $τ$ uit alle metingen en schat hoeveel watt de beker per graad verschil verliest (thee en beker samen $C \approx 1, 3 kJ/K$ , een schatting). Hoe goed doet Euler het met een stap van 4 minuten, de stap van de meting?

Boven: meting, model en Euler met de meetstap. Onder: ln(T − 21 °C) ligt op een rechte met helling −1/τ.

Gegeven

metingen om de 4 min, $T_{k} = 21, 0^{\circ} C$
$C \approx 1, 3 kJ/K$

Gevraagd

de vergelijking en haar oplossing
$τ$ en $h A$
Euler met h = 4 min

Oplossing

1
Balans: de enige warmtestroom gaat naar de kamer, evenredig met het verschil: $C \dot{T} = - h A (T - T_{k})$ . Dat is de afkoelingswet van Newton (1701). Dus $\dot{T} = - (T - T_{k}) / τ$ met $τ = C / (h A)$ .
2
Oplossing (vorige blad): $T (t) = T_{k} + (T_{0} - T_{k}) e^{- t / τ}$ .
3
Lineariseren: $ln (T - T_{k}) = ln (T_{0} - T_{k}) - t / τ$ is een rechte in $t$ . De kleinste-kwadratenrechte door de elf punten (hoofdstuk 3 van I, onder in de figuur) heeft helling $- 0, 0430$ per minuut, dus $τ = 23, 2 min$ , en snijpunt $T_{0} - T_{k} = 63, 9$ °C (gemeten 64,0).
4
$h A = C / τ = 1300 J/K / (23, 2 \cdot 60 s) = 0, 93 W/K$ : bij 50 °C verschil verliest de beker ongeveer 47 W.
5
Euler met $h = 4 min$ : $T_{n + 1} = T_{n} - \frac{4}{23 , 2} (T_{n} - 21)$ . Per stap blijft $1 - 4/ τ = 0, 828$ van het verschil over, in plaats van $e^{- 4/ τ} = 0, 842$ . Na 40 minuten geeft dat 30,7 °C tegenover 32,4 °C exact (gemeten 32,5 °C).

Antwoord

τ \approx 23 min

(het verschil met de kamer halveert elke 16 minuten),

h A \approx 0, 9 W/K

. Euler met de meetstap zit na 40 minuten 1,7 °C te laag: de stap is te groot tegenover

τ

Klopt dit? De residuen van het model zijn hoogstens 0,22 °C en wisselen van teken zonder patroon: het model van Newton past. Bij een open beker is het begin meestal iets sneller dan het model (verdamping); hier zie je dat niet.

Uitgewerkt voorbeeld · 4.7

De RC-kring: een vergelijking uit Kirchhoff, getoetst aan een meting

7/21

In hoofdstuk 24 van Intermediate laadde een Arduino een condensator (opdruk 100 µF, echt 94 µF) via 10 kΩ op tot 5,0 V en mat om de 0,1 s met de 10-bits ADC. Stel de differentiaalvergelijking op met de wetten van Kirchhoff, los ze op, en simuleer met Euler met de meetstap van 0,1 s. Is die simulatie goed genoeg om met de meting te vergelijken?

De meetpunten liggen op de exacte laadcurve; Euler met h = 0,1 s ligt er zichtbaar boven, met h = 0,3 s nog meer.

Gegeven

$R = 10 kΩ$ , $C = 94 μ F$ , $U_{0} = 5, 0 V$ , $u (0) = 0$
ADC-stap $5/1024 = 4, 9 mV$

Gevraagd

de vergelijking en $u (t)$
Euler met h = 0,1 s en de fout

Oplossing

1
Maaswet (hoofdstuk 21 van I): $U_{0} = R i + u$ . Stroom in een condensator (hoofdstuk 1): $i = C d u / d t$ . Samen: $R C \frac{d u}{d t} = U_{0} - u$ .
2
Dezelfde vorm als bij Newton, met $u_{\infty} = U_{0}$ en $τ = R C = 1 0^{4} \cdot 94 \cdot 1 0^{- 6} = 0, 94 s$ : $u (t) = U_{0} (1 - e^{- t / τ})$ . De meetpunten liggen op die kromme, op de afronding van de ADC na (minder dan één stap).
3
Euler: $u_{n + 1} = u_{n} + \frac{h}{τ} (U_{0} - u_{n})$ . Met $h = 0, 1 s$ blijft per stap $1 - h / τ = 0, 894$ van wat ontbreekt over, exact is dat $e^{- h / τ} = 0, 899$ . Euler laadt dus te snel.
4
Na 1,0 s: Euler $3, 376 V$ , exact $3, 274 V$ : 102 mV te hoog, 21 ADC-stappen. De grootste afwijking is 102 mV (bij $t \approx 0, 9 s$ ). Met $h = 0, 01 s$ wordt het 10 mV: tien keer kleiner (orde 1, verderop).

Antwoord

R C \overset{u}{˙} = U_{0} - u

, dus

u = U_{0} (1 - e^{- t / R C})

met

τ = 0, 94 s

. Euler met de meetstap zit tot 102 mV naast de exacte kromme, 21 keer de meetresolutie: te grof om met deze meting te vergelijken.

Klopt dit? Vuistregel: de fout van een simulatie moet veel kleiner zijn dan de meetfout, anders vergelijk je je rekenfout met de werkelijkheid. De stap van de meting (0,1 s) heeft niets te maken met de stap die de simulatie nodig heeft; die volgt uit

τ

en de gewenste nauwkeurigheid.

Uitleg · 4.8

Logistische groei: groeien tot een plafond

8/21

Geen groei is eeuwig exponentieel: algen nemen elkaar het licht af, een voedingsstof raakt op (hoofdstuk 5 van I). In 1838 stelde de Brusselse wiskundige Pierre-François Verhulst een eenvoudige correctie voor, om de groei van een bevolking te beschrijven: laat de relatieve groeisnelheid lineair dalen tot nul bij een grens $K$ , de draagkracht. In 1845 noemde hij de kromme logistique.

logistische vergelijking en haar oplossing

\frac{d N}{d t} = r N (1 - \frac{N}{K}) N (t) = \frac{K}{1 + \frac{K - N _{0}}{N _{0}} e ^{- r t}}

r

: groeisnelheid bij kleine

N

(per dag);

K

: draagkracht (zelfde eenheid als

N

). Voor

N ≪ K

\dot{N} \approx r N

: exponentieel. Evenwichten:

N = 0

(onstabiel) en

N = K

(stabiel).

Afleiding. Scheiden: $\frac{d N}{N ( 1 - N / K )} = r d t$ . Splits de breuk in partiële breuken: $\frac{1}{N ( 1 - N / K )} = \frac{K}{N ( K - N )} = \frac{1}{N} + \frac{1}{K - N}$ (tel maar op). Integreren geeft $ln N - ln (K - N) = r t + c$ , dus $\frac{N}{K - N} = \frac{N _{0}}{K - N _{0}} e^{r t}$ . Los op naar $N$ : de formule in het kader.

Het snelste punt. De groeisnelheid $f (N) = r N - r N^{2} / K$ is een bergparabool in $N$ met top waar $f^{'} (N) = r - 2 r N / K = 0$ : bij $N = K /2$ , met $\dot{N}_{max} = r K /4$ . Daar heeft $N (t)$ een buigpunt: de S-kromme gaat van hol naar bol. De relatieve groeisnelheid $\dot{N} / N = r (1 - N / K)$ is een dalende rechte in $N$ : dat is de sleutel om $r$ en $K$ uit metingen te halen (volgend blad).

In het richtingsveld: onder $K$ stijgen de lijnstukjes, boven $K$ dalen ze (te veel algen voor het licht: ze sterven af tot $K$ ). Van heel dicht bij 0 duurt het lang voor er iets te zien is; dan gaat het snel, en bij $K$ vlakt het af.

Alle oplossingen met N(0) > 0 lopen naar K; de onderste start met de algen van het testreservoir en is S-vormig.

Uitgewerkt voorbeeld · 4.9

Habitat: de draagkracht van het testreservoir

9/21

In het doorzichtige testreservoir (hoofdstuk 5 van I, opdracht 5.13) groeiden de algen tot dag 6 exponentieel en daarna trager. Bepaal $r$ en $K$ uit de relatieve groeisnelheid, vergelijk het logistische model met de metingen, en bereken wanneer het alarm ( $5 \cdot 1 0^{5}$ cellen per mL) afgaat volgens het exponentiële en volgens het logistische model. Hoe goed doet Euler het met een stap van één dag?

De relatieve groeisnelheid tegen N: waar de rechte de horizontale as snijdt, ligt de draagkracht K.

Gegeven

$N$ per dag: 0,019; 0,038; 0,082; 0,143; 0,307; 0,540; 1,160; 1,960; 3,300; 4,410; 5,850 ( $\times 1 0^{5}$ per mL, dag 0 tot 10)
alarm bij $5 \cdot 1 0^{5}$ per mL

Gevraagd

$r$ en $K$
model tegenover meting
tijdstip van het alarm
Euler met h = 1 dag

Oplossing

1
Relatieve groeisnelheid: $\dot{N} / N = d (ln N) / d t$ , met een centrale differentie over twee dagen (hoofdstuk 1): $(ln N_{i + 1} - ln N_{i - 1}) /2$ .
dag 1 2 3 4 5 6 7 8 9
$N$ (10⁵ /mL) 0,038 0,082 0,143 0,307 0,540 1,16 1,96 3,30 4,41
groei (/dag) 0,73 0,66 0,66 0,66 0,66 0,64 0,52 0,41 0,29
2
Tot $N \approx 1 0^{5}$ blijft ze rond 0,66 per dag, daarna daalt ze. Een rechte $\dot{N} / N = r - (r / K) N$ door de negen punten (kleinste kwadraten) geeft $r = 0, 704$ per dag en snijdt de $N$ -as in $K = 7, 68 \cdot 1 0^{5}$ per mL.
3
Model: $N (t) = \frac{7 , 68 \cdot 1 0 ^{5}}{1 + 403 e ^{- 0, 704 t}}$ . Het ligt hoogstens 10 % naast elke meting, over een factor 300 in $N$ .
4
Alarm: $N (t) = 5 \cdot 1 0^{5}$ geeft $t = \frac{1}{r} ln (\frac{K - N _{0}}{N _{0}} \cdot \frac{5 \cdot 1 0 ^{5}}{K - 5 \cdot 1 0 ^{5}}) = 9, 4$ dagen. Het exponentiële model met dezelfde $r$ zegt $ln (5 \cdot 1 0^{5} /1900) / r = 7, 9$ dagen. Gemeten: tussen dag 9 en 10.
5
Euler met $h = 1$ dag: in het begin groeit $N$ per stap met de factor $1 + r = 1, 70$ in plaats van $e^{r} = 2, 02$ . Op dag 10 geeft Euler $2, 97 \cdot 1 0^{5}$ in plaats van $5, 68 \cdot 1 0^{5}$ ; met $h = 0, 1$ dag $5, 44 \cdot 1 0^{5}$ .

Antwoord

r \approx 0, 70

per dag (verdubbeling per dag),

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

per mL. Het logistische model voorspelt het alarm na 9,4 dagen, zoals gemeten; het exponentiële 1,5 dagen te vroeg.

Klopt dit? De draagkracht ligt maar 1,5 keer boven de alarmgrens: in dit reservoir begrenst het licht de algen al vroeg. Een lichtdicht reservoir (het besluit van hoofdstuk 5 van I) zet

r

bijna op nul, en dan komt geen enkele grens in zicht. Het buigpunt (

N = K /2

) valt volgens het model op dag 8,5.

Exponentieel past tot dag 6; het logistische model volgt ook de afvlakking. Euler met h = 1 dag blijft achter.

Uitleg · 4.10

De methode van Euler

10/21

De meeste vergelijkingen uit de praktijk, zoals die van de meetkas met haar zon en buitentemperatuur die met de klok veranderen, kan je niet scheiden. Dan rekent een computer de oplossing stap voor stap uit. Het idee komt uit hoofdstuk 1: dicht bij een punt valt een kromme bijna samen met haar raaklijn, $y (t + h) \approx y (t) + h \overset{y}{˙} (t)$ . En de helling $\overset{y}{˙}$ ken je in elk punt: dat is precies wat de differentiaalvergelijking zegt.

methode van Euler

t_{n + 1} = t_{n} + h y_{n + 1} = y_{n} + h f (y_{n}, t_{n})

Vertrek in

(t_{0}, y_{0})

, zet een stap

h

langs de raaklijn, reken in het nieuwe punt de nieuwe helling uit, en herhaal. Het resultaat is een gebroken lijn die de oplossing benadert.

In de figuur koelt het station af met stappen van 1,5 h en van 0,5 h. Elk lijnstuk heeft de helling van het richtingsveld in zijn beginpunt. Omdat de kromme bol (naar boven gekromd) daalt, schiet elke stap onder de echte oplossing door: Euler koelt het station te snel af. Met een kleinere stap ligt de gebroken lijn dichter bij de kromme.

Twee manieren om Euler te begrijpen die je al kent:

Als integraal. $y (t_{n}) = y_{0} + \int_{t_{0}}^{t_{n}} f d t$ , en Euler benadert die integraal met de linker Riemann-som (hoofdstuk 3): in elke strook de hoogte van het beginpunt.
Als vaste fractie per stap. Voor $\overset{y}{˙} = - y / τ$ geeft Euler $y_{n + 1} = (1 - h / τ) y_{n}$ : per stapje gaat er een vaste fractie $h / τ$ af, zoals bij de condensator van hoofdstuk 24 van I. Na $n = t / h$ stappen is $y_{n} = (1 - \frac{t}{n τ})^{n} y_{0}$ , en voor $n \to \infty$ is dat precies $e^{- t / τ} y_{0}$ .

Elk lijnstuk van Euler volgt de helling in zijn beginpunt; met h = 1,5 h schiet de simulatie ver onder de echte afkoeling.

Uitgewerkt voorbeeld · 4.11

Euler met de hand

11/21

Simuleer het afkoelende station, $\dot{T} = - (T + 8) /3$ met $T (0) = 20$ °C, met de methode van Euler en een stap van een half uur tot $t = 2 h$ . Vergelijk met de exacte oplossing en verklaar het teken van de fout.

Gegeven

$f (T, t) = - (T + 8) /3$ (°C/h)
$T_{0} = 20^{\circ} C$ , $h = 0, 5 h$

Gevraagd

$T_{1}$ tot $T_{4}$ en de fout

Oplossing

1
Stap 1: $f (T_{0}) = - (20 + 8) /3 = - 9, 333$ °C/h, dus $T_{1} = 20 + 0, 5 \cdot (- 9, 333) = 15, 333$ °C.

Stap 2:

f (T_{1}) = - (15, 333 + 8) /3 = - 7, 778

, dus

T_{2} = 11, 444

°C. Zo verder:

$n$	$t_{n}$ (h)	$T_{n}$ (°C)	$f (T_{n}) = - (T_{n} + 8) /3$ (°C/h)	exact (°C)	fout (°C)
0	0,0	20,000	−9,333	20,000	0,000
1	0,5	15,333	−7,778	15,701	−0,368
2	1,0	11,444	−6,481	12,063	−0,618
3	1,5	8,204	−5,401	8,983	−0,779
4	2,0	5,503	−4,501	6,376	−0,873

3
Handig: elke stap vermenigvuldigt het verschil met buiten met $1 - h / τ = 1 - 0, 5/3 = \frac{5}{6}$ . Na vier stappen: $28 \cdot (5/6)^{4} = 13, 503$ , dus $T_{4} = 5, 503$ °C. Exact: $28 e^{- 2/3} = 14, 376$ .
4
De fout is negatief en groeit: elke stap gebruikt de helling aan het begin van de stap, en die is steiler dan gemiddeld over de stap (de afkoeling vertraagt).

Antwoord

T_{4} = 5, 50

°C tegenover

6, 38

°C exact: 0,87 °C te koud na vier stappen.

Klopt dit? De verhouding

5/6 = 0, 833

tegenover

e^{- 1/6} = 0, 846

per stap verklaart alles. Met een half zo grote stap wordt de fout ongeveer half zo groot: dat bekijkt het volgende blad.

Code · 4.12

Euler in Python

12/21

Een functie euler(f, y0, t0, t1, h) werkt voor elke vergelijking: je geeft $f$ mee als Python-functie. Ernaast de verbeterde Euler of methode van Heun: zet eerst een proefstap met Euler, reken ook daar de helling uit, en gebruik het gemiddelde van beide hellingen (de trapeziumregel van hoofdstuk 3 in plaats van de linker Riemann-som). Het programma rekent het station uit tot $t = 3 h$ en deelt de fout door $h$ en door $h^{2}$ .

Pythoneuler.py38 regelsDownload

# euler.py - de methode van Euler, en de verbeterde Euler (Heun) als voorsmaakje van hoofdstuk 36import numpy as np  def euler(f, y0, t0, t1, h):    """Lost dy/dt = f(y, t) op van t0 tot t1 met stap h; geeft de tijden en de waarden terug."""    n = round((t1 - t0) / h)    t = t0 + h * np.arange(n + 1)    y = np.empty(n + 1)    y[0] = y0    for k in range(n):        y[k + 1] = y[k] + h * f(y[k], t[k])                 # een stap langs de raaklijn    return t, y  def heun(f, y0, t0, t1, h):    """Gemiddelde van de helling aan het begin en (geschat) aan het eind van elke stap."""    n = round((t1 - t0) / h)    t = t0 + h * np.arange(n + 1)    y = np.empty(n + 1)    y[0] = y0    for k in range(n):        k1 = f(y[k], t[k])        k2 = f(y[k] + h * k1, t[k] + h)        y[k + 1] = y[k] + h / 2 * (k1 + k2)    return t, y  def station(T, t):                     # het station na een stroomuitval: dT/dt in graden per uur    return -(T + 8) / 3  exact = -8 + 28 * np.exp(-1)          # T(3 h) exactprint("    h     Euler T(3)    fout    fout/h   Heun: fout   fout/h²")for h in (1.0, 0.5, 0.25, 0.1, 0.01):    fe = euler(station, 20, 0, 3, h)[1][-1] - exact    fh = heun(station, 20, 0, 3, h)[1][-1] - exact    print(f"{h:6.2f} {exact + fe:12.4f} {fe:8.4f} {fe / h:8.3f} {fh:12.2e} {fh / h**2:8.3f}")

Uitvoer:

    h     Euler T(3)    fout    fout/h   Heun: fout   fout/h²
  1.00       0.2963  -2.0043   -2.004     2.47e-01    0.247
  0.50       1.3771  -0.9235   -1.847     5.42e-02    0.217
  0.25       1.8559  -0.4447   -1.779     1.27e-02    0.203
  0.10       2.1265  -0.1741   -1.741     1.96e-03    0.196
  0.01       2.2834  -0.0172   -1.719     1.91e-05    0.191

Bij Euler is fout/h vrijwel constant: de fout is evenredig met $h$ . Die constante kan je zelfs voorspellen: per stap is de fout ongeveer $\frac{1}{2} h^{2} \ddot{T}$ , en opgeteld over $t / h$ stappen geeft dat voor dit model $- \frac{h}{2} \cdot \frac{t}{τ ^{2}} \cdot 28 e^{- t / τ} = - 1, 717 h$ (opdracht 4.13). Bij Heun is fout/h² constant: tien keer kleinere stap, honderd keer kleinere fout. Heun rekent twee hellingen per stap, maar dat loont: met $h = 0, 1$ (60 hellingen) is de fout 0,0020 °C, kleiner dan die van Euler met $h = 0, 01$ (300 hellingen, 0,017 °C).

numpy en scipy

In de praktijk schrijf je Euler zelden zelf. scipy.integrate.solve_ivp kiest de methode (standaard een Runge-Kutta van orde 4-5) én de stap zelf, zodat de fout onder een tolerantie blijft. Wie begrijpt wat Euler doet, weet wat zo'n functie moet bewaken: de fout per stap en de stabiliteit (hoofdstuk 36).

Uitleg · 4.13

Stapgrootte en fout: Euler is van orde 1

13/21

Waarom halveert de fout als je de stap halveert? De reeks van Taylor (hoofdstuk 1) geeft de echte oplossing over één stap:

y (t + h) = y (t) + h \overset{y}{˙} (t) + \frac{1}{2} h^{2} \overset{y}{¨} (t) + \dots

Euler houdt alleen de eerste twee termen. De lokale fout per stap is dus ongeveer $\frac{1}{2} h^{2} \overset{y}{¨}$ . Over een vast interval $T$ zet je $T / h$ stappen, en die fouten tellen op (en planten zich voort):

Euler: orde 1

globale fout \approx \frac{T}{h} \cdot \frac{1}{2} h^{2} ∣ \overset{y}{¨} ∣ \sim h

Halveer

h

: de fout halveert, het rekenwerk verdubbelt. Tien keer nauwkeuriger kost tien keer zoveel stappen. Heun heeft een lokale fout

\sim h^{3}

en een globale fout

\sim h^{2}

(orde 2); de klassieke Runge-Kutta van hoofdstuk 36 is van orde 4.

Op een log-loggrafiek is een fout $\sim h^{p}$ een rechte met helling $p$ (hoofdstuk 1 deed hetzelfde voor differenties). Euler en Heun op het afkoelende station tonen hellingen 1 en 2. Anders dan bij numeriek differentiëren speelt de afronding hier bijna geen rol: pas bij miljarden stappen tellen de afrondingsfouten van $1 0^{- 16}$ merkbaar op.

Een kleine stap maakt een simulatie juist

Een kleine stap maakt de numerieke fout klein: het verschil tussen de simulatie en de exacte oplossing van het model. Is het model zelf fout (een ontbrekende warmtestroom, een verkeerde $τ$ ), dan reken je heel precies het verkeerde uit. Beide fouten schat je apart: de numerieke door de stap te halveren, de modelfout door te vergelijken met metingen.

Euler: tien keer kleinere stap, tien keer kleinere fout (helling 1). Heun: honderd keer kleinere fout (helling 2).

Uitleg · 4.14

Instabiliteit: als de stap te groot is

14/21

Neem opnieuw $\overset{y}{˙} = - y / τ$ . Euler geeft $y_{n + 1} = (1 - h / τ) y_{n}$ , dus $y_{n} = (1 - h / τ)^{n} y_{0}$ . Wat er gebeurt, hangt volledig af van de factor $q = 1 - h / τ$ :

$0 < h < τ$ : $0 < q < 1$ , een dalende rij zonder tekenwissel (te snel, maar redelijk);
$h = τ$ : $q = 0$ , na één stap is alles weg;
$τ < h < 2 τ$ : $- 1 < q < 0$ , de rij springt om beurten boven en onder nul, maar dempt uit: een afkoelend station dat afwisselend kouder en warmer wordt dan buiten;
$h > 2 τ$ : $∣ q ∣ > 1$ , de rij explodeert, terwijl de echte oplossing naar nul gaat.

stabiliteit van Euler bij verval

∣1 - h / τ ∣ < 1 ⟺ 0 < h < 2 τ

Voor een systeem met meer tijdconstanten beslist de kleinste. Nauwkeurig wordt het pas bij

h ≪ τ

;

h < 2 τ

is alleen de grens waaronder het niet ontploft.

Een Habitat-voorbeeld: het station heeft $τ = 3 h$ , maar zijn temperatuursensor heeft een tijdconstante van 60 s (hoofdstuk 39). Wie station en sensor samen simuleert met Euler, moet een stap onder 120 s nemen, hoewel het station zelf uren nodig heeft om te veranderen. Zo'n stelsel met sterk verschillende tijdconstanten heet stijf. De simulatie van hoofdstuk 39 rekent daarom met stappen van 10 s. Hoofdstuk 36 toont impliciete methoden die voor stijve stelsels veel grotere stappen toelaten (opdracht 4.14).

Een simulatie die ontploft, heeft een fout model

Niet noodzakelijk: het model kan perfect zijn en de methode met die stap onstabiel. Halveer de stap: verdwijnt het probleem, dan was het numeriek. Hetzelfde geldt op een microcontroller: het exponentiële filter $y \leftarrow y + α (x - y)$ van hoofdstuk 21 is een Euler-stap met $α = Δ t / τ$ , en is alleen stabiel voor $0 < α < 2$ .

Met h = 0,5 τ daalt Euler netjes, met 1,5 τ springt hij van teken maar dempt uit, met 2,5 τ explodeert hij.

Code · 4.15

Tweede orde en stelsels: dezelfde Euler met vectoren

15/21

De tweede wet van Newton is een vergelijking van de tweede orde: $m \overset{x}{¨} = F (x, \overset{x}{˙})$ . Voer de snelheid in als tweede onbekende, dan heb je twee vergelijkingen van de eerste orde:

\overset{x}{˙} = v, \overset{v}{˙} = F (x, v) / m .

Schrijf de toestand als vector $y = (x, v)$ ; dan is het stelsel weer $\dot{y} = f (y, t)$ en is Euler letterlijk dezelfde regel: $y_{n + 1} = y_{n} + h f (y_{n}, t_{n})$ . Met numpy werkt dezelfde code voor vectoren.

Voorbeeld. Een pakket van 2,0 kg (een doos van 30 cm) valt uit een bevoorradingsdrone van 40 m hoogte (hoofdstuk 9). De luchtweerstand is $F = \frac{1}{2} ρ c_{d} A v^{2} = k v^{2}$ met $ρ = 1, 2 kg/ m^{3}$ , $c_{d} \approx 1, 05$ en $A = 0, 09 m^{2}$ : $k = 0, 0567 kg/m$ . Naar beneden positief:

\overset{x}{˙} = v, \overset{v}{˙} = g - \frac{k}{m} v ∣ v ∣.

De snelheid stopt met stijgen als de weerstand het gewicht opheft: de eindsnelheid $v_{e} = m g / k = 18, 6 m/s$ . (De vergelijking voor $v$ alleen is scheidbaar: $v = v_{e} tanh (g t / v_{e})$ , zie de naslag.)

Pythonpakket.py15 regelsDownload

import numpy as np m, g, k = 2.0, 9.81, 0.0567           # massa (kg), zwaarteveld (N/kg), k = ½ · ρ · cd · A (kg/m)  def f(y, t):                          # y = [gevallen afstand x (m), snelheid v (m/s)], naar beneden positief    x, v = y    return np.array([v, g - k / m * v * abs(v)])  y, t, h = np.array([0.0, 0.0]), 0.0, 0.1while y[0] < 40.0:                    # tot het pakket 40 m gevallen is    y = y + h * f(y, t)               # dezelfde Euler-stap, nu met een vector    t = t + hprint(f"na {t:.1f} s: {y[0]:.1f} m gevallen, v = {y[1]:.1f} m/s")

Uitvoer:

na 3.5 s: 41.1 m gevallen, v = 17.8 m/s

Exact valt het pakket 3,41 s en raakt het de grond met 17,6 m/s; zonder lucht zou het na 2,86 s met 28,0 m/s neerkomen. Euler met $h = 0, 1 s$ zit er een tiende van een seconde naast (de simulatie stopt aan de eerste stap voorbij 40 m). Zo simuleer je in hoofdstuk 8 banen, in hoofdstuk 17 twee warmteknopen en in hoofdstuk 32 prooien en predatoren: telkens een vector en dezelfde stap.

Zonder lucht stijgt v rechtlijnig; met luchtweerstand buigt ze af naar de eindsnelheid. Euler met h = 0,5 s volgt grof.

Naslag · 4.16

Differentiaalvergelijkingen: formules en werkwijze

16/21

Werkwijze: (1) balans opstellen, eenheden controleren; (2) evenwichten zoeken ( $f = 0$ ) en het richtingsveld schetsen; (3) exact oplossen als het kan (scheiden), anders simuleren; (4) numerieke fout controleren door de stap te halveren; (5) vergelijken met metingen.

beginwaardeprobleem

\frac{d y}{d t} = f (y, t), y (t_{0}) = y_{0}

scheiden van veranderlijken

\frac{d y}{d t} = g (t) h (y) \Rightarrow \int \frac{d y}{h ( y )} = \int g (t) d t

lineair met evenwicht

y (t) = y_{\infty} + (y_{0} - y_{\infty}) e^{- t / τ}

symbool	betekenis	eenheid
$τ$	tijdconstante: C/(UA), RC, V/Q	s
$y_{\infty}$	evenwichtswaarde	eenheid van y

halverings- en verdubbelingstijd

t_{1/2} = τ ln 2 t_{2} = \frac{ln 2}{r} g = e^{k Δ t}

logistisch

N (t) = \frac{K}{1 + \frac{K - N _{0}}{N _{0}} e ^{- r t}} \dot{N}_{max} = \frac{r K}{4} bij N = \frac{K}{2}

symbool	betekenis	eenheid
$r$	groeisnelheid	1/s, 1/dag
$K$	draagkracht	eenheid van N

val met kwadratische weerstand

v (t) = v_{e} tanh \frac{g t}{v _{e}}, v_{e} = m g / k

Euler (orde 1)

y_{n + 1} = y_{n} + h f (y_{n}, t_{n})

symbool	betekenis	eenheid
$h$	tijdstap	s

Heun (orde 2)

y_{n + 1} = y_{n} + \frac{h}{2} (f (y_{n}, t_{n}) + f (y_{n} + h f_{n}, t_{n + 1}))

stabiliteit van Euler

0 < h < 2 τ_{min}

systeem	tijdconstante	waarde in Habitat
station (één knoop)	$τ = C / U A$	540 kJ/K / 50 W/K = 3 h
meetkas (snelle knoop)	$τ = C / U A$	30 kJ/K / 16,7 W/K ≈ 30 min
klimaatmodule 's nachts	$τ = C / U A$	1,1 MJ/K / 20,2 W/K ≈ 15 h
RC-meting (hoofdstuk 24 van I)	$τ = R C$	10 kΩ · 94 µF = 0,94 s
filter op de batterijmeting	$τ = (R_{1} ∥ R_{2}) C$	24,8 kΩ · 100 nF = 2,5 ms
temperatuursensor (hoofdstuk 39)	$τ_{s}$	60 s
CO₂ in de woonmodule	$τ = V / Q$	62 m³ / 50 m³/h ≈ 1,2 h

Code · 4.17

Habitat: de kastemperatuur simuleren

17/21

Nu het echte werk: het brein van het station moet kunnen voorspellen hoe warm de meetkas wordt. Het model is de warmtebalans van het tweede blad, met één temperatuur voor de hele kas (één knoop):

C \frac{d T}{d t} = a G (t) + P - U A (T - T_{b} (t)) .

$U A = 16, 7 W/K$ : het warmteverlies van de meetkas, nagerekend in hoofdstuk 19 van I.
$T_{b} (t)$ : de gemeten buitentemperatuur van de zonnige aprildag (hoofdstuk 4 van F), lineair geïnterpoleerd tussen de uurwaarden.
$G (t)$ : de instraling op een horizontaal vlak bij heldere hemel, met het eenvoudige zonnemodel van hoofdstuk 3 (20 april, Brussel). De klok van het station staat het hele jaar op wintertijd (hoofdstuk 30 van I): de zon staat om 12:42 in het zuiden.
$a$ (m²): de effectieve zonnewinst, alsof $a$ m² zonlicht volledig in kaswarmte omgezet wordt. Bovengrens: het dak (0,72 m²) maal de doorlaat van de platen (80 %) is 0,58 m².
$C$ : de warmtecapaciteit. Alles samen (tabel) is ongeveer 248 kJ/K, maar het water onder de tafel en de diepere grond volgen de dag niet. De snelle massa (lucht, platen, frame, bakjes, planten en de bovenste centimeter grond) is ongeveer 29 kJ/K.

deel van de meetkas	$C = m c$ (kJ/K)
water in het reservoir (24 kg)	100,5
vochtige potgrond (21 kg water, 22 kg droog)	109,9
tafel, hout (15,1 kg)	24,2
12 bakjes, PP (3,6 kg)	6,8
platen, PC (2,9 kg)	3,5
frame, aluminium (3,1 kg)	2,8
lucht (0,69 kg)	0,7
samen	248

Het programma simuleert eerst drie keuzes voor $C$ en $a$ en zoekt daarna in een raster welke combinatie de uurwaarden het best volgt, met de rms-fout (de wortel uit het gemiddelde kwadraat van de verschillen) als maat.

Pythonkas.py53 regelsDownload

# kas.py - eerste simulatie van de meetkas op de zonnige aprildag (klok in wintertijd, UTC+1)import mathimport numpy as np # uurwaarden van 0:00 tot 24:00 (hoofdstuk 4 van Fundamental)KAS = [9.5, 9.0, 8.6, 8.2, 7.9, 7.7, 7.7, 8.2, 11.0, 15.5, 20.5, 25.0, 28.5, 31.0, 32.0, 30.5, 27.5, 23.5,       19.5, 16.0, 13.5, 12.0, 11.0, 10.3, 9.8]BUITEN = [8.5, 8.0, 7.6, 7.2, 6.9, 6.7, 6.6, 7.2, 8.6, 10.4, 12.3, 14.0, 15.4, 16.5, 17.3, 17.8, 17.9, 17.3,          16.0, 14.3, 12.6, 11.3, 10.3, 9.5, 8.8]UA = 16.7                                # warmteverlies van de meetkas (W/K)  def zon(t, dag=110, lat=50.8, lon=4.35):    """Instraling op een horizontaal vlak bij heldere hemel (W/m²); t = kloktijd in wintertijd (h)."""    d = math.radians(23.44 * math.sin(math.radians(360 / 365 * (284 + dag))))            # declinatie    B = math.radians(360 / 365 * (dag - 81))    tv = 9.87 * math.sin(2 * B) - 7.53 * math.cos(B) - 1.5 * math.sin(B)                 # tijdvereffening (min)    w = math.radians(15 * (t - 1 + lon / 15 + tv / 60 - 12))                            # uurhoek    s = (math.sin(math.radians(lat)) * math.sin(d)         + math.cos(math.radians(lat)) * math.cos(d) * math.cos(w))                     # sinus van de zonnehoogte    if s <= 0:        return 0.0    direct = 1361 * 0.7 ** ((1 / max(s, 0.02)) ** 0.678)                                # loodrecht op de stralen    return direct * s + 0.12 * direct * s + 20 * math.sqrt(s)                           # direct + diffuus  def simuleer(C, a, dt=60.0, T0=KAS[0], P_verw=0.0):    """Euler voor C dT/dt = a G(t) + P - UA (T - T_buiten(t)), een dag lang."""    n = round(24 * 3600 / dt)    t = np.arange(n + 1) * dt / 3600                          # tijd in uren    T = np.empty(n + 1)    T[0] = T0    for k in range(n):        T_buiten = np.interp(t[k], range(25), BUITEN)        P = a * zon(t[k]) + P_verw - UA * (T[k] - T_buiten)    # netto warmtestroom (W)        T[k + 1] = T[k] + dt * P / C                           # Euler-stap    return t, T  def rms(T, dt=60.0):    """Wortel uit het gemiddelde kwadraat van (model - meting) op de hele uren."""    return math.sqrt(np.mean((T[::round(3600 / dt)] - np.array(KAS)) ** 2))  for naam, C, a in (("alle massa, heel het dak ", 248e3, 0.58), ("snelle massa, heel het dak", 30e3, 0.58),                   ("snelle massa, a = 0,30   ", 30e3, 0.30)):    t, T = simuleer(C, a)    print(f"{naam}: max {T.max():5.1f} °C om {t[T.argmax()]:5.2f} h, min {T.min():4.1f} °C, rms {rms(T):4.2f} °C") print("\nrms (°C)   a = " + "  ".join(f"{a:.2f}" for a in (0.20, 0.25, 0.30, 0.35, 0.40)))for C in (20e3, 30e3, 40e3, 60e3, 100e3, 150e3, 250e3):    fouten = [rms(simuleer(C, a, dt=120.0)[1], dt=120.0) for a in (0.20, 0.25, 0.30, 0.35, 0.40)]    print(f"C = {C / 1000:3.0f} kJ/K  " + "  ".join(f"{f:4.2f}" for f in fouten))

Uitvoer:

alle massa, heel het dak : max  34.8 °C om 16.13 h, min  7.6 °C, rms 6.56 °C
snelle massa, heel het dak: max  43.4 °C om 13.67 h, min  6.7 °C, rms 6.42 °C
snelle massa, a = 0,30   : max  30.5 °C om 14.03 h, min  6.7 °C, rms 0.82 °C

rms (°C)   a = 0.20  0.25  0.30  0.35  0.40
C =  20 kJ/K  2.39  1.36  0.81  1.45  2.50
C =  30 kJ/K  2.50  1.47  0.82  1.36  2.38
C =  40 kJ/K  2.67  1.66  1.01  1.38  2.33
C =  60 kJ/K  3.08  2.19  1.61  1.72  2.43
C = 100 kJ/K  4.03  3.31  2.85  2.77  3.09
C = 150 kJ/K  5.06  4.49  4.10  3.95  4.07
C = 250 kJ/K  6.46  6.02  5.70  5.51  5.48

Met alle massa is de kas veel te traag: de top valt om 16:08 in plaats van 14:00 en de fout is 6,5 °C. Met de snelle massa en het hele dak wordt ze 43 °C: veel te warm. Met $C = 30 kJ/K$ en $a = 0, 30 m^{2}$ zakt de fout onder 1 °C. Het raster toont een lang dal: voor elke $C$ tussen 20 en 40 kJ/K is $a \approx 0, 30 m^{2}$ het best, en de fout verandert er nauwelijks. Maar zodra $C$ groter wordt, helpt geen enkele $a$ nog. We kiezen $C = 30 kJ/K$ (de schatting van de snelle massa) en $a = 0, 30 m^{2}$ : de tijdconstante van de kas is dan $τ = C / U A = 30 min$ .

Waarom een stap van een minuut?

Met $τ \approx 30 min$ is Euler stabiel voor stappen onder een uur. Met 45 minuten ( $h = 1, 5 τ$ ) schommelt de simulatie al; met 75 minuten ( $2, 5 τ$ ) explodeert ze, tot 1010 °C en −1495 °C (figuur). Een minuut is ruim veilig en kost de Pi niets: 1440 stappen per dag.

Hetzelfde model, drie stappen: met 1 min glad, met 45 min schokkerig, met 75 min onstabiel.

Ontleed · 4.18

Habitat: model en meting naast elkaar

18/21

Boven: de uurwaarden van de aprildag, de simulatie met $C = 30 kJ/K$ en $a = 0, 30 m^{2}$ , de simulatie met alle massa en het hele dak, en de buitentemperatuur. Onder: de instraling waarmee het model rekent.

Het eenvoudige model volgt de dag op ongeveer een graad; waar het afwijkt, vertelt het wat er nog ontbreekt.

Wijs een nummer aan: de figuur zoomt in. Klik om vast te zetten.

1
De nacht
Zonder zon zakt het model tot op de buitentemperatuur (gemiddeld 0,8 °C te koud tussen 1 en 6 uur). De echte kas blijft ongeveer 1 °C warmer: de grond en het water geven 's nachts warmte af, en die trage massa zit niet in het model.
2
Opwarmen
Het model warmt het snelst op om 9:11: 4,7 °C/h. Hoofdstuk 1 haalde uit de uurwaarden 5,0 °C/h om 9:27. De ochtend klopt goed.
3
De top
Model 30,5 °C om 14:02, meting 32,0 °C om 14:00: het model blijft 1,5 °C onder de top. Mogelijk meet de sensor in de zon iets te hoog (stralingsfout), of krijgt de kas rond de middag meer zon dan een horizontaal vlak.
4
De avond
Om 18 uur is het model 1,8 °C te warm: de echte kas koelt sneller af. Kandidaten: de verdamping van ongeveer 1,8 L water per dag (gemiddeld 89 W overdag) en de trage massa die 's middags warmte opneemt.
5
Alle massa: te traag
Met C = 248 kJ/K (τ = 4,1 h) en het hele dak piekt het model pas om 16:08 op 34,7 °C: een kas die de hele massa in één keer moet opwarmen, loopt uren achter.
6
De zon
Om 12:42 staat de zon in het zuiden: 785 W/m². Zonsopgang 5:46, zonsondergang 19:37 (wintertijd); samen 6,2 kWh per m² op een horizontaal vlak.

Uitleg · 4.19

Simulatie tegenover meting: wat vertellen de verschillen?

19/21

Een model vergelijk je niet op het oog. Drie gereedschappen:

1. Eén getal: de rms-fout. $\frac{1}{n} \sum (T_{model} - T_{meting})^{2}$ vat de afwijking samen. De niveaulijnen tonen haar over alle combinaties van $C$ en $a$ (uit dezelfde simulatie). De kleinste fout, 0,80 °C, ligt in een lang, smal dal: een grotere $C$ met een iets grotere $a$ doet het bijna even goed. Deze meting kan $C$ en $a$ dus niet goed scheiden. Wie $C$ echt wil kennen, meet een avond zonder zon: dan geldt $\dot{T} = - (T - T_{b}) / τ$ en geeft de afkoeling rechtstreeks $τ = C / U A$ (opdracht 4.8).

2. De residuen: model min meting, uur per uur. Ze vertellen meer dan hun rms. Hier hebben ze een patroon: 's nachts te koud, rond de middag te koud, 's avonds te warm. Toevallige meetfouten wisselen willekeurig van teken; een patroon betekent dat er iets in het model ontbreekt, en geen enkele keuze van $C$ en $a$ haalt het weg (de rms zakt nergens onder 0,8 °C).

3. Fysische verklaringen. Wat ontbreekt er?

Trage massa. Grond en water nemen overdag warmte op en geven ze 's nachts af: warmere nachten, een koelere namiddag. Dat vraagt een tweede knoop met zijn eigen temperatuur (hoofdstuk 17).
Verdamping. 1,8 L water per dag kost $1, 8 kg \cdot 2454 kJ/kg = 4, 4 MJ$ , overdag gemiddeld 89 W: een flinke extra warmtestroom naar buiten, vooral als het warm is.
De zon op de wanden. Laag in de ochtend en de avond valt de zon door de zijwanden, rond de middag door het dak; een horizontaal vlak is maar een benadering.
De sensor. Een temperatuursensor in de zon meet te hoog; de houder van hoofdstuk 32 van F heeft daarom een dakje.

Fitten en testen op verschillende dagen

We hebben $C$ en $a$ gekozen op dezelfde dag waarop we het model beoordelen. Eerlijker is: parameters schatten op één dag en het model testen op een andere (een bewolkte dag, een warme dag). Dat is de kern van hoofdstuk 36 (curve fitting) en hoofdstuk 41 (train en test). De statisticus George Box vatte het samen: alle modellen zijn fout, maar sommige zijn nuttig.

De rms-fout over C en a: een lang dal met de kleinste fout bij kleine C en a ≈ 0,30 m²; alle massa met het hele dak ligt ver buiten het dal.

Residuen met een patroon: 's nachts en rond de middag te koud, 's avonds te warm. Dat is geen ruis.

Uitgewerkt voorbeeld · 4.20

Habitat: halen we eis K7 met watt of met water?

20/21

Eis K7 van de meetkas (wens, hoofdstuk 30 van F) vraagt dat de kas 's nachts minstens 1 °C warmer blijft dan buiten. Het model met $C = 30 kJ/K$ en $a = 0, 30 m^{2}$ haalt dat niet. Vergelijk twee oplossingen: (a) een verwarming die van 21:00 tot 6:00 werkt, (b) extra water in donkere jerrycans in de zon, 20 L of 40 L. Simuleer drie keer dezelfde dag na elkaar en beoordeel de nacht van 22:00 tot 6:00 op de derde dag (dan is de begintoestand vergeten).

Hoe meer water in de zon, hoe vlakker de dag: minder warm overdag, warmer 's nachts (model met één knoop).

Gegeven

$U A = 16, 7 W/K$ , $C = 30 kJ/K$ , $a = 0, 30 m^{2}$
water: 4,19 kJ/K per liter
energie meetkas: 44 tot 63 Wh per dag, batterij ±250 Wh (hoofdstuk 21 en 33 van F)

Gevraagd

verwarming: vermogen en energie
water: effect volgens het model
een onderbouwde keuze

Oplossing

1
Verwarming: 1 °C extra vraagt in evenwicht $P = U A \cdot 1 K = 16, 7 W$ , negen uur lang: $150 Wh$ per nacht. Dat is twee à drie keer het hele dagbudget van de meetkas en meer dan de helft van de batterij. Elektrisch verwarmen valt af.

Water: 20 L voegt

20 \cdot 4, 19 = 84 kJ/K

toe, 40 L

167 kJ/K

. Het model met die grotere

C

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

3
Het water neemt overdag warmte op (de top zakt tot 27 °C met 40 L) en geeft ze 's nachts af: met 40 L blijft de kas de hele nacht minstens 1,4 °C warmer dan buiten. Bonus: minder oververhitting, dus minder raambewegingen.
4
Maar dit is een bovengrens. Het model met één knoop veronderstelt dat water en lucht altijd dezelfde temperatuur hebben; in werkelijkheid geeft een jerrycan zijn warmte maar traag af (twee knopen, hoofdstuk 17). Het water moet overdag in de zon staan en 's nachts aan de lucht: tussen de plantenbakjes is geen plaats, dus bv. zwarte flessen langs de achterwand.

Antwoord

Watt: 17 W en ongeveer 150 Wh per nacht, onhaalbaar met een zonnepaneel van 50 Wp. Water: volgens het model haalt ongeveer 40 L in de zon K7 ruim, en de kas wordt overdag minder heet. Advies: een proef met 20 tot 40 L in donkere flessen, gemeten met de logger, en een model met twee knopen.

Klopt dit? Hoofdstuk 16 van F schatte met een ruwe rekening dat 24 L die 's nachts 5 °C afkoelt, de kas ongeveer 0,8 °C warmer houdt; het model geeft met 20 L gemiddeld 1,5 °C over de hele nacht: dezelfde orde. Opdracht 4.11 werkt de keuze uit.

Donkere flessen water langs de achterwand: overdag vangen ze zon, 's nachts geven ze warmte af.

Het geheel · 4.21

Waar differentiaalvergelijkingen vandaan komen en waar ze heen gaan

21/21

Een differentiaalvergelijking verbindt de afgeleide (hoofdstuk 1) met het integreren (hoofdstuk 3): ze zegt hoe snel iets verandert, en de oplossing telt die veranderingen op. Ze is de taal van elk dynamisch model, van warmte (E) en elektronica (L) tot groei in de biologie (B) en regelen (R). Simuleren is programmeren (P) in dienst van die modellen.

I5 Groei en verval: exponentiële modellen

Groeifactor, verdubbelings- en halveringstijd, de thee en de algen: hier als dy/dt = ky.

I18 Warmte berekenen

Warmtecapaciteit C = m · c en de nachtbuffer van de module.

I19 Warmtetransport en isolatie

UA van de meetkas (16,7 W/K) en de module (20,2 W/K).

I24 Condensatoren en de RC-tijd

De RC-meting en de vaste fractie per stapje: de voorloper van Euler.

E1 Afgeleiden: verandering meten

De raaklijn als benadering is de Euler-stap; centrale differenties schatten r uit de algen.

E3 Integralen: ophopen

Scheiden van veranderlijken is integreren; Euler is een linker Riemann-som, Heun een trapeziumregel.

E17 Warmteoverdracht kwantitatief

Het warmtemodel met twee knopen: de trage massa die hier ontbrak.

E32 Ecologische modellen

Logistische groei en Lotka-Volterra: stelsels van groeivergelijkingen.

E36 Numerieke methoden en wetenschappelijk rekenen

Runge-Kutta, stapgrootteregeling, impliciete methoden voor stijve stelsels, parameters fitten.

E38 Dynamische systemen modelleren

Eerste- en tweedeordesystemen, tijdconstante en stapresponsie.

E39 PID-regeling

De PID-simulatie van het station is een Euler-simulatie met τ = 3 h.

E49 Integratieproject: het brein van het station

De digitale tweeling: een stelsel differentiaalvergelijkingen dat meeloopt met de metingen.

HabitatHet brein van het station heeft nu een eerste model van de meetkas:

C \dot{T} = a G (t) + P - U A (T - T_{b})

met

U A = 16, 7 W/K

C \approx 30 kJ/K

(snelle massa) en

a \approx 0, 30 m^{2}

, gesimuleerd met Euler in stappen van een minuut; het volgt de aprildag op ongeveer een graad. Het toont dat eis K7 met water haalbaar is en niet met elektriciteit. Hoofdstuk 17 voegt de trage knoop toe, hoofdstuk 36 fit de parameters op echte logs. Opdracht 4.11 werkt de buffer uit.

Samenvatting

In het kort

Een differentiaalvergelijking $d y / d t = f (y, t)$ met een beginwaarde beschrijft een verandering; je stelt ze op als balans (in − uit) en controleert de eenheden.
Het richtingsveld toont in elk punt de helling; oplossingen raken eraan en kruisen elkaar niet. Evenwichten liggen waar $f = 0$ .
Scheiden van veranderlijken geeft $y = y_{0} e^{k t}$ en $y = y_{\infty} + (y_{0} - y_{\infty}) e^{- t / τ}$ (Newton, RC, station); $t_{1/2} = τ ln 2$ .
Logistische groei $\dot{N} = r N (1 - N / K)$ vlakt af naar de draagkracht $K$ ; de relatieve groei $\dot{N} / N$ is lineair in $N$ , zo vind je $r$ en $K$ uit metingen.
Euler: $y_{n + 1} = y_{n} + h f (y_{n}, t_{n})$ , stappen langs de raaklijn; werkt ook voor stelsels (vectoren).
De globale fout van Euler is evenredig met $h$ (orde 1), die van Heun met $h^{2}$ . Euler op verval is alleen stabiel voor $h < 2 τ$ ; de kleinste tijdconstante beslist.
Een simulatie vergelijk je met metingen via de rms-fout én de residuen: een patroon in de residuen wijst op een ontbrekend mechanisme.
De meetkas volgt $C \dot{T} = a G + P - U A (T - T_{b})$ met $C \approx 30 kJ/K$ , $a \approx 0, 30 m^{2}$ en $τ \approx 30 min$ ; eis K7 haal je met water, niet met watt.

Begrippen

Wat je nu kent

differentiaalvergelijking: Een vergelijking tussen een onbekende functie en haar afgeleide(n), bv. $d y / d t = f (y, t)$ : ze zegt hoe snel $y$ verandert in elke toestand.
beginwaardeprobleem: Een differentiaalvergelijking samen met een beginwaarde $y (t_{0}) = y_{0}$ ; als $f$ netjes is, ligt de oplossing daarmee vast.
evenwicht: Een constante oplossing $y^{*}$ met $f (y^{*}) = 0$ . Stabiel als oplossingen in de buurt ernaartoe gaan, onstabiel als ze ervan weglopen.
richtingsveld: Een tekening met in elk punt $(t, y)$ een kort lijnstuk met helling $f (y, t)$ ; oplossingen raken overal aan die lijnstukjes.
scheiden van veranderlijken: Oplossingsmethode voor $d y / d t = g (t) h (y)$ : alle $y$ naar één kant, alle $t$ naar de andere, en beide kanten integreren.
tijdconstante: De tijd $τ$ waarin een exponentieel proces $d y / d t = - (y - y_{\infty}) / τ$ nog 37 % ( $1/ e$ ) van zijn afstand tot het evenwicht over heeft.
logistische groei: Groei volgens $d N / d t = r N (1 - N / K)$ : eerst exponentieel, daarna afremmend naar de draagkracht $K$ (Verhulst, 1838).
methode van Euler: Numerieke oplossing met stappen langs de raaklijn: $y_{n + 1} = y_{n} + h f (y_{n}, t_{n})$ .
orde van een methode: De macht van $h$ waarmee de globale fout daalt: Euler orde 1 (fout $\sim h$ ), Heun orde 2 (fout $\sim h^{2}$ ).
numerieke stabiliteit: Een methode is stabiel als fouten niet aangroeien. Euler op $\overset{y}{˙} = - y / τ$ is stabiel voor $h < 2 τ$ .
stijf stelsel: Een stelsel met sterk verschillende tijdconstanten: de snelste bepaalt de toegelaten stap, ook als ze niet interessant is.
residu: Het verschil tussen model en meting in één punt. Een patroon in de residuen wijst op een ontbrekend mechanisme.

Aan de slag in het werkboek16 opdrachten, van oefenen tot uitdagen

dag	1	2	3	4	5	6	7	8	9
$N$ (10⁵ /mL)	0,038	0,082	0,143	0,307	0,540	1,16	1,96	3,30	4,41
groei (/dag)	0,73	0,66	0,66	0,66	0,66	0,64	0,52	0,41	0,29