Anfangswertprobleme/Methoden zur Lösung von Anfangswertproblemen/Integration der Differentialbeziehung (IVP)

Was ohne Computer geht ...

Ohne Computer interessieren uns meist nur Anfangswertprobleme, bei denen wir eine analytische Lösung des Anfangswertproblems angeben können.

Integration der Bewegungsgleichung

Die Bewegungsdifferentialgleichung eines Körpers im Erdschwerefeld sei - ohne Einflüsse der Luft - für die Koordinaten Wurf-Weite w und Wurf-Höhe h:

$\begin{array}{ll}\ddot{h}(t)& =0\\ \ddot{w}(t)& =-g\end{array}$

Diese Differentialgleichung können wir - wie auch die Biegedifferentialgleichung des Euler-Bernoulli-Balkens - integrieren, und finden

$\begin{array}{l}h(t)=-\frac{g}{2}\cdot t^2+c_{h,1}\cdot t+c_{h,0},\\ w(t)=c_{w,1}\cdot t+c_{w,0}\end{array}$

mit den Integrationskonstanten c_h1, c_h0, c_w1, c_w0.

Die Integrationsbedingungen können wir an "Randbedingungen" zu Zeitpunkten t₁, t₂ anpassen, meist aber interessiert uns die Lösung, wenn der Anfangszustand für vollständig durch Anfangsweg und Anfangsgeschwindigkeit beschrieben ist.

Die Lösungen sehen dann so aus:

/* equations of motion */
eom: ['diff(w(t),t,2)= 0,
      'diff(h(t),t,2)= -g];
 
/* solve eom */
wurf : solve(integrate(integrate(eom,t),t),[h(t),w(t)])[1];
wurf: subst([%c1 = V[w,0], %c3 = W[0],
             %c2 = V[h,0], %c4 = H[0]],wurf);
 
/* parameter */
par: [H[0] = 2*m, W[0] = 0*m,
      V[h,0] = 10*m/s, V[w,0] = 3*m/s,
      g=10*m/s^2];
 
/* plot results */
toPlot : append([parametric],
  subst([t=t*s],subst(par,subst(wurf, [w(t),h(t)]/m))),
  [[t,0,2.2]]);
/* Wurf-Parabel */
plot2d(toPlot, [legend, "wurf"],
  [style, [lines,3]],
  [x,0,8], [xlabel, "Weite w/m →"],
  [y,0,8], [ylabel, "Höhe h/m →"])$
/* Weg-Zeit-Diagramm */
plot2d([toPlot[2], toPlot[3]],[t,0,2.2],
  [legend, "w(t)", "h(t)"],
  [style, [lines,3]],
  [xlabel, "Zeit t/s →"],
  [ylabel, "Weg h/m, w/m →"])$
/* Geschwindigkeits-Zeit-Diagramm */
plot2d(diff([toPlot[2], toPlot[3]],t),[t,0,2],
  [legend, "dw/dt", "dh/dt"],
  [style, [lines,3]],
  [xlabel, "Zeit t/s →"],
  [ylabel, "Geschwindigkeit dh/dt/(m/s), dw/dt/(m/s) →"])$

Wurfparabel

Weg-Zeit-Diagramm

Geschwindigkeits-Zeit-Diagramm.

Diese Lösungsvariante spielt in Ingenieursproblemen praktisch keine Rolle, weil die Bewegungsgleichungen sich nicht geschlossen integrieren lassen.

Der e^λt-Ansatz

Meist sind die Zustandsgrößen (Auslenkung, Geschwindigkeit) der Körper unseres Systems elastisch miteinander oder an eine Umgebung gekoppelt, die Koordinaten u kommen dann in der Bewegungsgleichungen in verschiedenen Zeitableitungen (Geschwindigkeit, Beschleunigung) vor. Dann sehen die Bewegungsgleichungen so aus:

$\underset{\_}{\underset{\_}{M}}\cdot \underset{\_}{\ddot{u}}+\underset{\_}{\underset{\_}{K}}\cdot \underset{\_}{u}=\underset{\_}{f}$

mit

$\begin{array}{ll}\underset{\_}{u}& :\text{ Koordinaten des Systems}\\ \underset{\_}{\underset{\_}{M}}& :\text{ Massenmatrix}\\ \underset{\_}{\underset{\_}{K}}& :\text{ Steifigkeitsmatrix und}\\ \underset{\_}{f}& :\text{ eingeprägte Lasten auf das System}\end{array}$.

Für diese Systeme interessiert uns sehr oft nur die homogene Lösung (f=0). Für dieses System können wir die Lösung mit den Ansatz

$\underset{\_}{u}(t)=\widehat{\underset{\_}{u}}\cdot e^{\lambda \cdot t}$

ansetzen (erraten). Dieser Ansatz führt auf ein Eigenwertproblem mit den Eigenwerten λ und den Eigenvektoren u.

Zwei-Massen-Schwinger.

Das schauen wir uns an dem schwingungsfähige System links an. Es besteht aus zwei Körper der Masse m sowie Federn der Steifigkeit k1=3/2 k und k2=k. Gesucht ist die Schwingung des Systems in Zeitbereich, wenn die untere Masse durch einen Schlag eine Anfangsgeschwindigkeit erfährt.

Beispiel: Zwei-Massen-Schwinger

Koordinaten der Bewegung.

Zur Lösung mit dem Prinzip der virtuellen Verrückungen und den Prinzip von d'Alembert brauchen wir zunächst Koordinaten der Bewegung:

Die Gleichgewichtsbedingung lautet

$\begin{array}{ll}\delta W& =\underset{:=\delta W^a}{\underbrace{-m\cdot {\ddot{w}}_1\cdot \delta w_1-m\cdot {\ddot{w}}_2\cdot \delta w_2}}-\underset{:=\delta \mathrm{\Pi }}{\underbrace{k_1\cdot w_1\cdot \delta w_1+k_2\cdot (w_2-w_1)\cdot (\delta w_2-\delta w_1)}}\\ & \stackrel{!}{=}0\end{array}$

Aufgelöst nach den Koeffizienten der virtuellen Verrückungen folgt die gewöhnliche, lineare Differentialgleichung:

$\left(\begin{array}{cc}m& 0\\ 0& m\end{array}\right)\cdot \left(\begin{array}{c}{\ddot{w}}_1\\ {\ddot{w}}_2\end{array}\right)+\left(\begin{array}{cc}\frac{5}{2}k& -k\\ -k& k\end{array}\right)\cdot \left(\begin{array}{c}{w}_1\\ w_2\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)$

Dimensionslos machen wir die Bewegungsgleichung mit 

$\tau =\frac{t}{T}\text{ und }T=\sqrt{\frac{m}{k}}$

Die Bezugszeit T ist das 2π-fache der Schwingungsperiode eines Ein-Masse Schwingers mit Masse m und Federsteifigkeit k.

Der e^λτ-Ansatz mit

$\left(\begin{array}{c}{w}_1\\ w_2\end{array}\right)=\underset{\_}{\hat{q}}\cdot e^{\lambda \cdot \tau }$

liefert die nun dimensionslosen Eigenwerte

${\lambda }_1=\frac{i}{\sqrt{2}},{\lambda }_2=\sqrt{3}\cdot i\text{ mit }i=\sqrt{-1}$.

Aus dem charakteristischen Polynom 4.ter Ordnung kommen eigentlich 4 Nullstellen und damit Eigenwerte des Systems. Hier unterscheiden sich jeweils paarweise zwei Eigenwerte nur durch das Vorzeichen, was beim Quadrieren in der zweiten Zeitableitung allerdings zu gleichen Eigenwerten führt, nämlich

$\begin{array}{ll}\text{zu }{\lambda }_1=\frac{i}{\sqrt{2}}& :\hat{q}=\left(\begin{array}{c}1\\ 2\end{array}\right),\\ \text{zu }{\lambda }_2=\sqrt{3}\cdot i& :\hat{q}=\left(\begin{array}{c}1\\ -\frac{1}{2}\end{array}\right).\end{array}$

Die Lösung lautet also mit den Integrationskonstanten C_i

$\left(\begin{array}{c}{w}_1(t)\\ w_2(t)\end{array}\right)=C_1\cdot {\underset{\_}{\hat{q}}}_1\cdot e^{{\lambda }_1\cdot \tau }+C_2\cdot {\underset{\_}{\hat{q}}}_2\cdot e^{{\lambda }_2\cdot \tau }$

Die Schreibweise sieht einfach aus, verbirgt aber das Wesen der Lösung: die λs sind komplex-wertig und nur durch "Zufall" (wie bei allen Bewegungsgleichungen ohne erste Zeitableitung) sind die Eigenvektoren rein reell-wertig. Dafür dürfen die Cjs nun auch komplex-wertig sein, also

$C_1=c_{1,R}+i\cdot c_{1,I}\text{ , }{C}_2=c_{2,R}+i\cdot c_{2,I}$

und wir brauchen vier Anfangsbedingungen zu ihrer Bestimmung.

Die sind:

$\begin{array}{rccl}1)& w_1(0)& =& 0\\ 2)& w_2(0)& =& 0\\ 3)& {\dot{w}}_1(0)& =& 0\\ 4)& {\dot{w}}_2(0)& =& V_0\end{array}$

und daraus kommt

$\begin{array}{llc}{c}_{1,R}& =& 0,\\ c_{1,I}& =& -\frac{2^{\frac{3}{2}}}{5},\\ c_{2,R}& =& 0,\\ c_{2,I}& =& \frac{2}{5\cdot \sqrt{3}}\end{array}$.

/*******************************************************/
/* MAXIMA script                                       */
/* version: wxMaxima 15.08.2                           */
/* author: Andreas Baumgart                            */
/* last updated:                                       */
/* ref: Gross e.a Aufg. 1.6 (räumliches Gleichgewicht) */
/* description                                         */
/*                                                     */
/*******************************************************/

/* declare variational variables - see 6.3 Identifiers */
declare("δW", alphabetic);
declare("δΠ", alphabetic);
declare("δw", alphabetic);

/*******************************************************/
/* declare variational variables - see 6.3 Identifiers */
declare("δW", alphabetic);
declare("δΠ", alphabetic);
declare("δw", alphabetic);
assume(k>0,m>0);

/* values to make eom dimensionless */
dims: [W = V[0]*T, T=sqrt(m/k)];

/* coordinates */
q : [[w[1](t),w[2](t)],[δw[1],δw[2]]];
/*** equations of motion ***/
/* princ. of virt. work */
PvV: [δW^a = -m*q[2].diff(q[1],t,2),
      δΠ   = 3/2*k*w[1](t)*δw[1] + k*(w[2](t)-w[1](t))*(δw[2]-δw[1])];
eom: makelist(coeff(-expand(subst(PvV,δW^a-δΠ)),δw[i]),i,1,2);

/* system matrices */
M : funmake('matrix,makelist(makelist(coeff(eom[i],diff(w[j](t),t,2)),j,1,2),i,1,2));
K : funmake('matrix,makelist(makelist(coeff(eom[i],     w[j](t)     ),j,1,2),i,1,2));

/*** solve eigenvalue problem ***/

/* characteristic polynomial */
charPoly : determinant((%lambda/T)^2*M+K)=0;

/* eigenvalues */
eigenVals: subst(dims,solve(charPoly, %lambda));
eigenVals: [eigenVals[2],eigenVals[4]];

/* eigenvevtors */
Q : matrix([1],[W[2]]);
eigenVects: ratsimp(subst(dims,subst(solve(((
                       (%lambda/T)^2*M+K).Q)[1][1],W[2]),Q)));
eigenVects: ratsimp(makelist(subst(eigenVals[i], eigenVects),i,1,2));

/* time-domain-solution */
solution: W*sum(subst(eigenVals[i],(c[i,R]+%i*c[i,I])*eigenVects[i]*%e^(%lambda*tau)),i,1,2);
initVals: append(makelist(realpart(subst([tau=0],solution)[i][1])=0,i,1,2),
                 makelist(realpart(subst([tau=0],diff(solution,tau)))[i][1]/T=[0,V[0]][i],i,1,2));
coeffs: solve(subst(dims,initVals), [c[1,R],c[1,I],c[2,R],c[2,I]])[1];

/*** post-processing ***/
toPlot : makelist(realpart(subst(coeffs,solution[i][1]/W)),i,1,2);

/* Weg-Zeit-Diagramm */
plot2d(toPlot,[tau,0,2*2*%pi],
  [legend, "w[1]", "w[2]"],
  [style, [lines,3]],
  [xlabel, "Zeit τ →"],
  [ylabel, "Auslenkung w/W →"])$
/* Geschwindigkeit-Zeit-Diagramm */
plot2d(diff(toPlot,tau),[tau,0,2*2*%pi],
  [legend, "w[1]", "w[2]"],
  [style, [lines,3]],
  [xlabel, "Zeit τ →"],
  [ylabel, "Geschwindigkeit dw/dτ →"])$
/* Phasen-Diagramm */
plot2d([[parametric, subst(t,tau,toPlot[1]), subst(t,tau,diff(toPlot[1],tau)), [t,0,8*2*%pi]],
  [parametric, subst(t,tau,toPlot[2]), subst(t,tau,diff(toPlot[2],tau)), [t,0,2*2*%pi]]],
  [legend, "w[1]", "w[2]"],
  [style, [lines,3]],
  [xlabel, "Auslenkung w/W →"],
  [ylabel, "Geschwindigkeit dw/dτ →"])$

Aufgetragen in Diagrammen erhalten wir dann diese analytische Lösung des Anfangswertproblem.

✔ {{{title}}}:

{{{text}}}