Modellierung dynamischer Prozesse mit Differentialgleichungen (DGL) und deren Parameterschätzung mit dem Excel® Solver

Qualitätsüberprüfung

Um die Qualität des numerischen Verfahrens für die Integration von DGLs zu prüfen, wird das numerische Ergebnis den analytischen Lösungen von sattsam bekannten Funktionen, wie der logistischen Wachstumskurve oder auch der Richards Funktion, der numerischen Lösung der RKF Methode für einen gegebenen Parametersatz gegenübergestellt (Abb. 1). Die zu schätzenden Parameter sind neben dem Anfangswert der DGL zum Zeitpunkt t=0, die Wachstumsrate r, sowie der Kapazitätsterm K. Bei der Richardsfunktion kommt ein Exponent γ, ein Formparameter, hinzu. Die maximale Abweichung der numerischen gelösten DGLs liegt im Bereich von 1E-8.

Logistic & Richards function — Abb. 1: Vergleich der numerisch gelösten DGLs (Linien) zur exakten Lösung (Punkte)

Logistic growth equation — Gl. 1,2: Logistisches Wachstum und Richards Gleichung

The related VBA code:
dy(1) = r * y(1) * (1 - y(1) / K)
dy(2) = r * y(2) * (1 - (y(2) / K) ^ gamma)

Der zweite Vergleich ist Güte der geschätzten Parameter jeweils für die analytische als auch die numerische Lösung beider DGLs.

Tabelle 1: Vergleich der Parameter der Solver angepassten DGL-Modelle, zuerst für die analytische, gefolgt für die numerische Lösung

Parameter	Exakte Parameter	Anpassung an analytisches Modell	Anpassung an numerisch gelöstes Modell
y₀	2.2	2.20004	2.20030
r	0.1	0.10002	0.09999
K	40	39.99996	40.00018
γ	1.4	1.40001	1.40004

Die Ergebnisse zeigen minimal Abweichungen vom wahren Wert, wobei das Ergebnis für die numerische Lösung ein klein wenig schlechter erscheint im Vergleich zur analytischen Lösung für beide Modelle. Da zwei numerische Verfahren hintereinandergeschaltet sind, verbleibt es die Entscheidung des Anwenders, ob die erreichte Präzision für seine Bedürfnisse ausreichend ist.

Anmerkung: Es ist leicht möglich völlig falsche Anpassungsergebnisse zu erhalten, wenn man mit extrem ungeeigneten Startwerten den Schätzprozess startet. Man kann das System zu einem Zero Division Problem treiben, wenn z.B. die Wachstumsrate r größer ist als der Kapazitätsterm K, oder einer der Parameter Null approximiert. Da beide Gleichung ein chaotisches Verhalten oberhalb gewisser Parameterbereiche aufweisen, sind solche Bereiche zu vermeiden. Man kann das Problem durch entsprechende Einstellungen der Nebenbedingungen im Solver in Bezug auf die zu erwarten Modelltrajektorie verringern. Zuletzt, es wird vorgeschlagen den Fit zu wiederholen um die Stabilität der erreichten Lösung zu prüfen, da eventuell Zeit- und Iterationslimits im ersten Anpassungsvorgang überschritten wurden.