Anfangswertprobleme/Methoden zur Lösung von Anfangswertproblemen
Problemstellung
Fast immer sind im Maschinenbau Anfangswertprobleme über Differentialgleichungen im Zeitbereich definiert - die unabhängige Koordinate ist deshalb hier immer t. Für die gesuchte Funktion
suchen wir die Lösung zum Anfangswertproblem
im Zeitintervall 0 < t < T.
Numerische Löser für Anfangswertprobleme gibt es meist nur für Differentialgleichungen erster Ordnung (erste Ableitung nach der Zeit). Da Bewegungsgleichungen mechanischer Systeme fast immer Differentialgleichungen zweiter Ordnung sind (die Beschleunigung und damit die d'Alembert'sche Trägheitskrafte sind zweite Ableitungen des Weges nach der Zeit), brauchen wir einen Trick: Wir führen die Geschwindigkeit explizit als Koordinate ein, z.B. für die Bewegungsgleichung
Nun sind
die Zustandsgrößen des Systems.
Mit
schreiben wir die Differentialgleichung zweiter Ordnung als zwei Differentialgleichungen erster Ordnung:
Beispiel 1: Wachstum einer Population
Eine Fliegenpopulation besteht zum Zeitpunkt t=0 aus der Anzahl von n0 Fliegen.
In einem Zeitintervall stirbt ein Anteil λS von Fliegen, ein Anteil λG von Fliegen wird geboren, das Anfangswertproblem lautet:
Die analytische Lösung dazu ist
Beispiel 2: Wachstum einer Population mit Selbstvergiftung
Für sehr viele Fliegen n führen dies zu "Wachstum mit Selbstvergiftung". Die Stoffwechselrückstände r der Population in einem abgeschlossenen System (z.B. Erde) führen zum Aufbau von Hemmstoffen, die negativ auf die Geburtenrate λG wirken, hier mit dem Faktor (1-r).
Dieses Modell (die Bewegungsgleichung) für r setzt für die Änderungsgeschwindigkeit einen Prozess an, bei dem die Hemmstoffe
- proportional mit μ zur Anzahl der Fliegen n steigen und
- proportional zu r abgebaut werden (zerfallen):
Auch wenn es einfach aussieht: das IVP hat keine analytische Lösung mehr. Wir müssen es als Anfangswertproblem numerisch lösen. Und das geht so:
Zum Zeitpunkt t0 =0 wissen wir
Einsetzen in die Bewegungsgleichung liefert:
Die Änderungsgeschwindigkeit der Zustandsgrößen des Modells (n,r) ist also für t=t0 bekannt. Nun ersetzten wir - wie beim Verfahren der Finiten Differenzen - den Differentialquotienten durch den Differenzenquotienten:
Und analog für r.
/* equations of motion */
eom: ['diff(n,tau)=((1-r)-%lambda[S])*n,
'diff(r,tau)= mu*n- rho*r];
/* parameter */
par: [%lambda[S]=0.5, rho=1/10, mu=1/100];
/* solve IVP */
results: rk(subst(par,makelist(rhs(eom[i]),i,1,2)),[n,r],[1,0],[tau,0,50,0.1])$
results: [tau = makelist(results[i][1],i,1,length(results)),
n = makelist(results[i][2],i,1,length(results)),
r = makelist(results[i][3],i,1,length(results))]$
/* plot results */
plot2d([[discrete,subst(results,tau),subst(results,n)],
[discrete,subst(results,tau),subst(results,r)]],
[legend, "# Pop.", "Rückst."],
[xlabel, "τ→"],[ylabel, "n,r→"],
[style, [lines,3]])$
Klimamodelle machen genau das - nur mit deutlich komplizierteren Modellen.
Links
- Phasendiagramme: eine Methode zum Ausdeuten der Lösung von nunmerisch gelösten AWP.
Untergeordnete Seiten