Lösung von retardierten Differentialgleichungen zur Modellierung von … Murmeltieren?

Mit der Veröffentlichung der Version 6.3 der Software COMSOL Multiphysics^® ist es nun möglich, frühere Lösungen aus einer zeitabhängigen Studie im aktuellen Zeitschritt zu verwenden. Das heißt, dass die Lösung von einem früheren Zeitpunkt verwendet werden kann, um die Lösung zum aktuellen Zeitpunkt zu beeinflussen. Dies ermöglicht zahlreiche Anwendungen, unter anderem für die Modellierung von retardierten Differentialgleichungen. Sehen wir uns diese neue Funktionalität anhand eines etwas ungewöhnlichen Beispiels an.

Gedanken auf der Bergwiese

Vor einiger Zeit hatte ich die Gelegenheit, ein wenig in den Bergen zu wandern, und entdeckte ein Murmeltier. Nachdem ich ein paar Schnappschüsse gemacht hatte, begann ich über den Lebenszyklus dieses liebenswerten Pflanzenfressers nachzudenken. Meine Gedanken wanderten schnell zu einer Gleichung, mit der ich die Populationsveränderung von Murmeltieren beschreiben könnte.

Das flauschige Murmeltier, das diesen Blog-Beitrag inspiriert hat.

Ich dachte an meine Fachbücher über numerische Methoden aus dem Studium und erinnerte mich an alle möglichen Gleichungen, die nützlich sein könnten, darunter die Lotka-Volterra-Gleichungen und die logistische Differentialgleichung. In Anlehnung an diese Gleichungen stellen wir zwei gewöhnliche Differentialgleichungen (Ordinary Differential Equations, ODEs) auf: eine für die Änderungsrate der Anzahl der Murmeltiere auf der Wiese, M(t), und eine weitere für die Änderungsrate der Menge an essbarer Biomasse, B(t).

Die erste Gleichung, für die Änderungsrate der Murmeltierpopulation, enthält Terme für die Fortpflanzung, R(t), und die Mortalität, D(t), welche wir in Kürze genauer besprechen werden. Die zweite Gleichung hat drei Konstanten: G_0 ist die Wachstumsrate, B_{max} ist die maximale Masse an Vegetation, die auf der Wiese wachsen kann, und E_0 ist der tägliche Verbrauch der Murmeltiere.

Für die Gleichung für die Murmeltierpopulation benötigen wir einen Term für die Fortpflanzung, der eine lineare Abhängigkeit von der aktuellen Population berücksichtigt:

Die Sterberate der Murmeltiere sollte andererseits mit den verfügbaren Ressourcen in Zusammenhang stehen. Daher definieren wir als Zwischenschritt die Futterfläche pro Tag:

wobei A_0 die Größe der Wiese ist. Wir gehen davon aus, dass die Wahrscheinlichkeit, einem Prädator zu begegnen, mit der Größe der Futterfläche steigt. Wir gehen außerdem von einer festen Anzahl von Prädatoren pro Gebiet aus, P(t) = P_0. Die Sterberate der Murmeltiere ist also proportional zur aktuellen Anzahl der Murmeltiere, M(t), dem täglichen Futtergebiet jedes Murmeltiers und der Prädatorendichte. Dies ergibt eine Gleichung für die Sterberate, die von beiden Unbekannten abhängt:

Wir gehen außerdem davon aus, dass die gesamte Murmeltierkolonie auf einmal aus dem Winterschlaf erwacht, und wir lösen diese Gleichungen, um die Populationsveränderung im Sommer vorherzusagen, woraufhin sich die Kolonie für den Winter wieder in ihren Bau zurückzieht. Wir werden auch eine Stoppbedingung einbauen, um den Löser anzuhalten, wenn die Anzahl der Murmeltiere oder die Menge der Vegetation unter einen kritischen Wert fällt, was auf einen Zusammenbruch der Kolonie oder der Umwelt hindeutet.

Der Screenshot zeigt, wie man eine Reihe von ODEs modelliert. Verschiedene ODEs können unterschiedliche Einheiten haben. Ein Feature Stop Condition wurde ebenfalls zum Feature Time-Dependent Solver hinzugefügt.

Gleichungen dieser Art können in COMSOL Multiphysics^® mithilfe des Interfaces Global ODEs and DAEs gelöst werden, wie im Screenshot dargestellt ist. Das Interface definiert die beiden Variablen, die wir lösen, M(t) und B(t), sowie die Gleichung für jede in Bezug auf ihre Zeitableitungen und ihre Anfangswerte, M_\text{init}=M(t=0) und B_\text{init}=B(t=0). Wir können diese Gleichungen lösen, um vorherzusagen, wie sich die Population im Laufe des Sommers entwickelt.

Lösungen der Gleichungen für die Dynamik der Vegetationsbiomasse (grün) und der Murmeltierpopulation (braun), mit einer Anmerkung zur Population am Ende des Sommers. Die Murmeltierpopulation steigt zunächst steil an und fällt dann aufgrund erhöhter Prädatorenaktivität.

Die Ergebnisse zeigen, dass die Murmeltierpopulation zunächst stark ansteigt und die verfügbare Nahrung schnell verbraucht. Dies führt zu einer Korrektur, also einem Rückgang der Population aufgrund erhöhter Prädatorenaktivität. Nach einiger Zeit erreicht die Murmeltierpopulation ein Gleichgewicht (obwohl es durchaus schwierig sein kann, festzustellen, ob es sich hierbei um ein stabiles Gleichgewicht handelt).

Retardierte Differentialgleichungen zur Berücksichtigung zusätzlicher Faktoren

Bisher haben wir plausible Ergebnisse, aber mindestens zwei wichtige Faktoren sollten noch berücksichtigt werden:

1. Murmeltiere beginnen sich nach dem Ende des Winterschlafs zu paaren, pflanzen sich aber nicht sofort fort; es gibt eine Tragzeit.
2. Wenn die Murmeltierpopulation über einen ausreichend langen Zeitraum wächst, werden mehr Prädatoren auf die Wiese gelockt.

Im Folgenden wird gezeigt, wie diese beiden Faktoren mathematisch beschrieben und in das Modell integriert werden können. Um den ersten Faktor zu berücksichtigen, legen wir die Fortpflanzungsrate als Funktion der Population zum Zeitpunkt der Empfängnis fest und berücksichtigen die Tragzeit, T_g, in der die Fortpflanzungsrate gleich Null ist. Unsere Gleichung für die Fortpflanzungsrate lautet also:

wobei H(t-T_g) die Heaviside-Funktion ist. Die Rate wird auch durch den Bruch \frac{M\left( T_g \right)}{M_\text{init}} reduziert, also den gesamten Bevölkerungsrückgang während der Tragzeit. Wir gehen davon aus, dass die älteren Murmeltiere nicht mehr gejagt werden, sobald die Jungtiere mit der Nahrungssuche beginnen. Daher bleibt dieser Bruch nach t=T_g unverändert.

Wir haben nun eine retardierte Differentialgleichung, die in die Software eingegeben werden kann, indem ein if()-Ausdruck mit dem Operator at(time,variable) kombiniert wird. Wir können den Fortpflanzungsterm schreiben als:

R0*if(t > Tg,at(t-Tg,M)*at(Tg,M)/Minit,0)

Das erste Argument für den Operator at() ist der vorherige Zeitpunkt, zu dem wir das zweite Argument auswerten möchten. Um diese Funktion zu aktivieren, bei der während der Lösung auf vorherige Zeitpunkte zurückgegriffen wird, müssen wir die Löser-Einstellungen wie im folgenden Screenshot gezeigt ändern. Die Verwendung von Strict Time Stepping mit einem Ausgabezeitschritt, der nicht größer als die Verzögerungszeit ist, wird sichergestellt, dass wir den Beginn der Geburtensaison erfassen. Ansonsten wird das Modell auf die gleiche Weise wie zuvor gelöst.

Der Screenshot zeigt, wie die Verwendung von Zeitoperatoren während der Lösung aktiviert wird.

Nun betrachten wir die Population unter Einbeziehung dieser nicht konstanten Fortpflanzungsrate. Wir beobachten zunächst einen Populationsrückgang in dem Zeitraum, in dem keine Murmeltiere geboren werden, und dann einen Anstieg der Population, der sich allmählich abflacht.

Population von Murmeltieren (braun) und Menge an Vegetation (grün), wenn ein Verzögerungsterm in die Geburtenrate einfließt.

Schließlich wollen wir die Sterberate anpassen, um die Zunahme der Prädatoren zu berücksichtigen, die eintritt, wenn die durchschnittliche Anzahl der Murmeltiere über einen ausreichend langen Zeitraum ansteigt. Das heißt, wir müssen die zeitgemittelte Murmeltierpopulation der vorangegangenen Tage verfolgen.

Anstatt dieses Integral über alle vorherigen Zeitschritte über den Zeitraum T_p zu bewerten, können wir den Hauptsatz der Differentialrechnung anwenden und eine Gleichung für die Änderungsrate der durchschnittlichen Bevölkerung in den vorangegangenen Tagen aufstellen:

Diese Gleichung könnte einfach gelöst werden, indem unserem Modell eine weitere globale Gleichung hinzugefügt wird, aber sie würde ab t=0 ausgewertet werden, also verwenden wir einen weiteren Ausdruck if() zur Definition von M(t<0) = M_\text{init}. Wir verwenden den gleitenden Durchschnitt, um die Prädatorendichte P(t) = \left( P_0 /M_\text{init}\right) M_\text{ave}(t) zu erhöhen, was sich wiederum auf die Sterberate auswirkt. Löst man die Gleichung erneut, sieht man die Auswirkungen auf die Population am Ende des Sommers.

Gute Nachrichten! Die Population hat sich vergrößert und die Murmeltiere können sich für ihren langen Winterschlaf eingraben.

Die Ergebnisse zeigen die Murmeltierpopulation (braun) und die Biomasse (grün), wobei sowohl die Fortpflanzungsrate als auch die Sterberate von vorherigen Lösungen abhängen.

Experimentelle Beobachtungen bestätigen einen Anstieg der Murmeltierpopulation.

Eine kurze Anmerkung zur Modellierung des Bevölkerungsgleichgewichts

An dieser Stelle sollte betont werden, dass das hier vorgestellte Beispiel einige spannende neue Funktionen im Kontext eines leicht verständlichen Modells veranschaulichen soll. Tatsächlich können reale Populationsdynamikmodelle wesentlich komplexer sein und verwenden häufig einen kompartimentierten oder intervallbasierten Ansatz, bei dem verschiedene Variablen zur Verfolgung verschiedener Teile der Gesamtpopulation verwendet werden. Das SEIR-Modell aus der Epidemiologie (Susceptible, Exposed, Infectious, Recovered and immune) ist ein gutes Beispiel für diesen Ansatz und wird in unserem Blog-Beitrag „Modeling the Spread of COVID-19 with COMSOL Multiphysics^®“ diskutiert.

Dieser intervallbasierte Ansatz ist besonders für die Chemietechnik relevant, wo er häufig zur Bestimmung der Größe von Tröpfchen und Niederschlägen eingesetzt wird. Die folgenden Tutorial-Modelle veranschaulichen diese Art der Modellierung:

• Oil–Water Flow Through an Orifice — A Droplet Population Model
• Precipitation of Barium Sulfate
• Crystallization of Benzoic Acid in a Mixed Suspension, Mixed Product Removal Crystallizer

Schlussbemerkungen

Wir haben die Grundlagen der Implementierung und Lösung von retardierten Differentialgleichungen in COMSOL Multiphysics^® Version 6.3 gezeigt. Obwohl dieses einfache Modell einige sehr interessante Dynamiken erzeugen kann, ist es nur als Lernbeispiel gedacht. Wenn Sie dieses einfache Beispiel verstanden haben, sind Sie bereit, anspruchsvollere Modellierungsaufgaben anzugehen. Sie können retardierten Differentialgleichungen der Zeit etwa mit Differentialgleichungen im Raum kombinieren. Jedes Interface in der Software unterstützt diese neuen Operatoren. Unabhängig von der Kombination der Interfaces in Ihrer Multiphysik-Modellierung ist die Technik zum Abrufen früherer Lösungen genau dieselbe wie die hier gezeigte.

Möchten Sie das hier beschriebene Modell testen? Laden Sie die zugehörige MPH-Datei aus der Application Gallery herunter: