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

$\displaystyle q(t){\text{ und der Abkürzung }}{\frac {d}{dt}}(.):={\dot {(.)}}$

suchen wir die Lösung zum Anfangswertproblem ${\dot {q}}=f(q,t){\text{ mit dem Anfangswert }}q(0)=q_{0}$

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

${\ddot {u}}(t)=f(u,{\dot {u}},t)$

Nun sind

$u,{\dot {u}}$

die Zustandsgrößen des Systems.

Mit

$v(t)={\dot {u}}(t)$

schreiben wir die Differentialgleichung zweiter Ordnung als zwei Differentialgleichungen erster Ordnung:

$\underbrace {\left({\begin{array}{c}{\dot {v}}\\{\dot {u}}\end{array}}\right)} _{\displaystyle {\underline {\dot {q}}}}=\underbrace {\left({\begin{array}{c}f(u,v,t)\\v\end{array}}\right)} _{\displaystyle {\underline {f}}({\underline {q}},t)}$

Beispiel 1: Wachstum einer Population

Eine Fliegenpopulation besteht zum Zeitpunkt t=0 aus der Anzahl von n₀ Fliegen.

In einem Zeitintervall stirbt ein Anteil λ_S von Fliegen, ein Anteil λ_G von Fliegen wird geboren, das Anfangswertproblem lautet:

${\dot {n}}=\underbrace {(\lambda _{G}-\lambda _{S})} _{\displaystyle =:\lambda _{0}}\cdot n\;{\text{ mit dem Anfangswert }}n(0)=n_{0}$

Die analytische Lösung dazu ist

$\displaystyle n(t)=n_{0}\cdot e^{\lambda _{0}\cdot t}$

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):

${\begin{array}{ll}{\dot {n}}&=\left(\lambda _{G}\cdot (1-r)-\lambda _{S}\right)\cdot n\\{\dot {r}}&=\mu \cdot n-\varrho \cdot r\end{array}}\;{\text{ mit den Anfangswerten }}n(0)=n_{0},r(0)=0$

Implementierung in Maxima

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 t₀ =0 wissen wir

$\displaystyle \left({\begin{array}{c}n(0)\\r(0)\end{array}}\right)=\left({\begin{array}{c}n_{0}\\0\end{array}}\right)$

Einsetzen in die Bewegungsgleichung liefert:

$\displaystyle \left.\left({\begin{array}{c}{\dot {n}}\\{\dot {r}}\end{array}}\right)\right|_{t=0}=\left({\begin{array}{c}\lambda _{G}\cdot n_{0}\\\mu \cdot n_{0}\end{array}}\right)$

Die Änderungsgeschwindigkeit der Zustandsgrößen des Modells (n,r) ist also für t=t₀ bekannt. Nun ersetzten wir - wie beim Verfahren der Finiten Differenzen - den Differentialquotienten durch den Differenzenquotienten:

$n(\underbrace {t_{0}+\Delta t} _{\displaystyle :=t_{1}})=n(t_{0})+{\dot {n}}|_{t=t_{0}}\cdot \Delta t$

Und analog für _r. 

Integrationsschritt beim Lösen eines AWPs.

So berechnen wir aus dem Anfangspunkt den Zustand des Systems zu einem Folgepunkt, den Zeitpunkt t₁. Und für t₁ den Folgepunkt t₂'. Wenn wir diese Schritte oft genug machen, erhalten wir den approximierten Verlauf der Zustandsgrößen über der Zeit. Im Maxima-Modell machen wir das mit der dimensionslosen Zeit:

$\displaystyle \tau =\lambda _{G}\cdot t$

Für die gewählte Dauer von 50 Fliegen-Generationen erhalten wir:

Anzahl der Fliegen und Hemmstoffe über der Zeit.

/* wxMaxima 21.05.2                */
/* equations of motion */
eom: ['diff(n,τ)=((1-r)-λ[S])*n,
      'diff(r,τ)= μ*n- ρ*r];

/* parameter */
par: [λ[S]=0.5, ρ=1/10, μ=1/100];

/* solve IVP */
results: rk(subst(par,makelist(rhs(eom[i]),i,1,2)),[n,r],[1,0],[τ,0,50,0.1])$
results: [τ = 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,τ),subst(results,n)],
        [discrete,subst(results,τ),subst(results,r)]],
        [legend, "# Pop.", "Rückst."],
        [xlabel, "τ→"],[ylabel, "n,r→"],
        [style, [lines,3]])$

Implementierung in Matlab

Ganz analog erfolgt die Implementierung in Matlab. Nur dass wir hier zwei Dateien brauchen:

1. die Skriptdatei und
2. die Funktion, die dydt zurückliefert.

Anzahl der Fliegen und Hemmstoffe über der Zeit.

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% Modul Numerische Methoden der Mechanik                            %
% IVP - Introduction                                                %
% Autor: Andreas Baumgart                                           %
% Last Update: 2022-05-17                                           %
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

cd('C:\Users\abs384\OneDrive - haw-hamburg.de\2022-SS\Numerische Methoden der Mechanik\Unterricht\Matlab\KW21') % hier den richtigen Pfad eingeben!
addpath("./Functions");

%% parameter
lambda=0.5; rho=1/10; mu=1/100;

tspan = [0 50];
y0 = [1 0];

%% solve ODE
[t,y] = ode45(@(t,y) flies(t,y,lambda,rho,mu), tspan, y0);

%% plot
plot(t,y(:,1),t,y(:,2));
xlabel('Time τ');
ylabel('solutions n,r');
legend('n','r');

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% Function dydt(y,t)                                                %
%                                                                   %
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
function dydt = flies(t,y,lambda,rho,mu)
%RETURN rhs of ode
%   self-poisoning of population of flies

n=y(1);r=y(2);

dydt = [((1-r)-lambda)*n; ...
         mu*n- rho*r];

end

Klimamodelle machen genau das - nur mit deutlich komplizierteren Modellen.

Links

• Phasendiagramme: eine Methode zum Ausdeuten der Lösung von numerisch gelösten AWP.

Untergeordnete Seiten

1. Euler-Streckenzug-Verfahren
2. Integration der Differentialbeziehung (IVP)
3. Löser für steife Differentialgleichungen
4. Runge-Kutta-Verfahren 4.ter Ordnung