kramann.info

2.3 Programmierung des Eulerschen Integrationsverfahrens mit C++ für das Seerosenmodell

Die Steigung dA/dt hat die Einheit m²/s. Entsprechend

Die Steigung dA/dt hat die Einheit m²/s
k muss dann die Einheit 1/s haben
Es wird angenommen, dass die Seerosenfläche zu Beginn 20qm gross ist und sich nach einer Woche zu 40qm verdoppelt hat. Dann gilt nach der analythischen Lösung, die wir später noch betrachten werden: k = 1/872542s
Als Schrittweite Δt wollen wir 60s nehmen.
Unser konkretes Modell lautet dann:

Bild 2.3-1: Simulationsmodell-Gleichungen

Ein einfaches C/C++ - Programm liefert uns nun die tägliche Seerosenfläche ab dem Zeitpunkt t₀, an dem sie 20qm betrug, für einen Zeitraum von 30 Tagen:

#include<iostream>

using namespace std;


int main(void)
{
    double  A = 20.0;
    double dt = 60.0;
    double  k = 1.0/872542.0;

    double sekunden_pro_tag = 86400.0; 
    double tage       = 30.0;
    long   schritte   = (long)((tage*sekunden_pro_tag)/dt);
    long   i          = 0;

    long   schritte_pro_tag = (long)(sekunden_pro_tag/dt);

    double t = 0.0;

    do
    {
        if(i%schritte_pro_tag==0)
            cout<<"Tag Nr."<<i/schritte_pro_tag<<" Flaeche="<<A<<"qm"<<endl;        

        A *=k*dt+1.0; 

        i++;
    } while(i<=schritte);
}

Code 2.3-1: seerose.cpp

seerose.zip (C++ Programm)

Auf der Konsole wird nach dem Compilieren und dem Start des Programmes folgendes ausgegeben:

Bild 2.3-2: Konsolenausgabe mit Schrittweite dt=60.0s

Offenbar liefert der Integrator nach sieben Tagen ein leicht falsches Ergebnis. Satt genau 40qm werden nur 39,999qm errechnet. Es sei gleich gesagt, dass diese Differentialgleichung sehr gutmütig ist. Bei schwingenden Systemen mit sehr harten federnden Elementen (steife/stiff Differentialgleichungssysteme/Differential-Equations), würde der Fehler sehr viel größer ausfallen.

Eine Änderung der Schrittweite von dt=60s auf dt=1s liefert sofort das gewünschte Ergebnis:

Bild 2.3-2: Konsolenausgabe mit Schrittweite dt=1.0s