Gelöste Aufgaben/Kw24: Unterschied zwischen den Versionen

Aus numpedia
Zur Navigation springen Zur Suche springen
Keine Bearbeitungszusammenfassung
 
(2 dazwischenliegende Versionen desselben Benutzers werden nicht angezeigt)
Zeile 10: Zeile 10:


==Aufgabenstellung==
==Aufgabenstellung==
SOME TEXT
<onlyinclude>
[[Datei:Kw24-01.png|alternativtext=|links|mini|180x180px|Lageplan]]
Ein Mechanismus besteht aus einem starren Stab ''AC'' (Länge ''ℓ<sub>1</sub>'', Massenmoment ''J'') und einer Feder ''BC'' (Länge ''ℓ<sub>2</sub>'' im entspannten Zustand, Steifigkeit ''k'').
Ein Mechanismus besteht aus einem starren Stab ''AC'' (Länge ''ℓ<sub>1</sub>'', Massenmoment ''J'') und einer Feder ''BC'' (Länge ''ℓ<sub>2</sub>'' im entspannten Zustand, Steifigkeit ''k'').


Die Körper sind drehbar gelagert, der Abstand ''AB'' ist ''ℓ<sub>3</sub>.''
Die Körper sind drehbar gelagert, der Abstand ''AB'' ist ''ℓ<sub>3</sub>.''


Gesucht sind Lösungen als Anfangswertproblem, nichtlineare Schwingungen und Aussagen über die Stabilität der Bewegungen.</onlyinclude>
<onlyinclude>
[[Datei:Kw24-01.png|alternativtext=|links|mini|180x180px|Lageplan]]
Hier geht es um die verschiedenen Lösungstypen des Systems - in Abhängigkeit von den Anfangsbedingungen. Gesucht sind Lösungen als Anfangswertproblem, nichtlineare Schwingungen und Aussagen über die Stabilität der Bewegungen.</onlyinclude>


== Lösung mit Maxima ==
== Lösung mit Maxima ==
Lorem Ipsum ....


==tmp==
<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Header
{{MyCodeBlock|title=Header
|text=Text
|text=
Gesucht ist die Lösung zu dem nichtlinearen Anfangswertproblem.
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/*********************************************************/
/* MAXIMA script                                        */
/* version: wxMaxima 15.08.2                            */
/* author: Andreas Baumgart                              */
/* last updated: 2017-09-30                              */
/* ref: Kw24 (TM-C, Labor 6)                            */
/* description: finds the solution for                  */
/*              the nonlinear IVP                        */
/*********************************************************/
</syntaxhighlight>
</syntaxhighlight>
}}
}}


==tmp==
 
<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Declarations
{{MyCodeBlock|title=Declarations
|text=Text
|text=
Wir wählen die Längenverhältnisse
 
::<math>\displaystyle {{\ell}_{2}}={{\ell}_{1}}, {{\ell}_{3}}=\sqrt{2}\cdot {{\ell}_{1}}</math>
 
so dass die Feder für φ''(t) = 45''° entspannt ist.
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/*********************************************************/
/* declare variational variables - see 6.3 Identifiers  */
declare("δW", alphabetic); /* total virtual work        */
declare( "φ", alphabetic); /* coordinate of rotatin      */
declare("δφ", alphabetic); /* variation of coordinate phi*/
declare("δw", alphabetic); /* variation of coordinate w  */
 
/* parameter */
assume(l[1]>0, k>0, J>0);
params: [l[2] = l[1], l[3] = sqrt(2)*l[1]];
</syntaxhighlight>
</syntaxhighlight>
}}
}}


==tmp==
<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Kinematics
{{MyCodeBlock|title=Kinematics
|text=Text
|text=
Der Mechanismus hat genau einen Freiheitsgrad, den wir entweder durch die Federlängung ''w'' oder den Drehwinkel ''φ'' erfassen. Den Zusammen zwischen den beiden stellen wir über
 
::<math>{{\ell}_{1}^{2}}\cdot {{\mathrm{sin}\left( \mathrm{\varphi}\left( t\right) \right) }^{2}}+{{\left( {{\ell}_{3}}-{{\ell}_{1}}\cdot \mathrm{cos}\left( \mathrm{\varphi}\left( t\right) \right) \right) }^{2}}={{\left( {{\ell}_{2}}+w\right) }^{2}}</math>
 
her. Nun steht uns frei, ob wir diese Beziehung nach ''w'' oder ''φ'' auflösen. Wir wählen ''φ'' als Minimal-Koordinate und erhalten
 
::<math>w = \displaystyle \frac{{{\ell}_{1}}\cdot \sqrt{12-{{2}^{\frac{7}{2}}}\cdot \cos(\varphi) }-2\cdot {{\ell}_{1}}}{2}</math>
 
Nach den Regeln der [[Sources/Lexikon/Variationsmethoden|Variationsmethoden]] erhalten wir daraus:
 
::<math>\delta w = \displaystyle \frac{{2^{\frac{3}{2}}}\cdot {{\ell}_{1}}\cdot \sin( \varphi ) \cdot \delta \varphi}{\sqrt{12-{{2}^{\frac{7}{2}}}\cdot \cos(\varphi)}}</math>.
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/*********************************************************/
/**** kinematics ***/
 
/* coordinates */
coords : [φ, δφ];
 
kin : (l[1]*sin(φ(t)))^2+(l[3]-l[1]*cos(φ(t)))^2=(l[2]+w)^2;
kin : trigsimp(ratsimp(solve(subst(params, kin),w)[2]));
kin : [kin,
      δw = subst([epsilon=0],
                diff(subst([φ(t)=φ(t)+epsilon*δφ],
                                  subst(kin,w)),epsilon))];
</syntaxhighlight>
</syntaxhighlight>
}}
}}


==tmp==
 
<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Equilibrium Conditions
{{MyCodeBlock|title=Equilibrium Conditions
|text=Text
|text=
Mit dem Prinzip der virtuellen Verrückungen schreiben wir
 
::<math>\displaystyle \delta W=-\left( \frac{{{d}^{2}}}{d\,{{t}^{2}}}\cdot \mathrm{\varphi}\left( t\right) \right) \cdot J\cdot \delta \varphi-k\cdot w\cdot \delta w</math>
 
und erhalten die Bewegungsgleichung
 
::<math>\displaystyle J \cdot \ddot{\varphi} + \sqrt{2}\cdot l_1^2 \cdot k \cdot \sin\varphi - \frac{{{2}^{\frac{3}{2}}}\cdot l_1^2 \cdot k \cdot \sin\varphi}{\sqrt{12-{{2}^{\frac{7}{2}}}\cdot \cos\varphi}} = 0</math>
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/*********************************************************/
/* equilibrium conditions */
PvV : δW = -J*'diff(φ(t),t,2)*δφ - k*w*δw;
 
/* dgl for φ ist the coefficient of δφ in δW */
dgl[1] : expand(coeff(subst([δW = 0], subst(kin, -PvV)),δφ));
</syntaxhighlight>
</syntaxhighlight>
}}
}}


==tmp==
 
<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Equilibrium Points and Stability
{{MyCodeBlock|title=Equilibrium Points
|text=Text
|text=
Die Bewegungsgleichung machen wir dimensionslos, indem wir charakteristische Größen aus der linearisierten Bewegungsgleichung verwenden.
 
Die drei Gleichgewichts-Lagen bzgl. derer wir die Bewegungsgleichung linearisieren können, sind
 
::<math>\begin{array}{ll}\varphi_{0,1} &= \displaystyle -\frac{\pi}{4}\\\varphi_{0,2} &= 0\\\varphi_{0,3} &= \displaystyle +\frac{\pi}{4}\end{array}</math>.
 
Eine Taylor-Reihenentwicklung des Rückstellmoments bzgl. φ''<sub>0,3</sub>'' ,also
 
::<math>\varphi = \varphi_{0,3} + \phi</math>
 
liefert die linearisierte Bewegungsgleichung
 
::<math>\displaystyle\ddot{\phi}+\omega_0^2\cdot\phi = 0 \text{ mit } \omega_0^2 = \frac{{{\ell}_{1}^{2}}\cdot k}{J}</math>.
 
Dabei ist ''ω<sub>0</sub>'' die Eigenkreisfreuquenz der linearisierten Bewegungsgleichung.
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/*********************************************************/
/* Zeros of restoring function */
equilib : solve(subst(['diff(φ(t),t,2)=0],dgl[1]),φ(t));
 
/*********************************************************/
/* derive the eigenfrequency                            */
/* of linearized ODE at stable point of equilob.        */
LinODE : taylor (subst([φ(t)=psi],rhs(dgl[1]/J)), psi, %pi/4,1);
ref : omega[0]^2 = coeff(expand(LinODE),psi);
ref : [t[B] = 1/(subst(solve(ref,omega[0])[2],omega[0])/(2*%pi))];
</syntaxhighlight>
</syntaxhighlight>
}}
}}


[[Datei:Kw24-11.png|left|mini|Rückstellmoment der Feder]]


<!-------------------------------------------------------------------------------->{{MyCodeBlock|title=Make Equations of Motion Dimensionless
|text=
Eine Bezugszeit für unsere nichtlineare Bewegungsgleichung kommt nun aus


==tmp==
::<math>\displaystyle \omega_0 = \frac{2\pi}{t_{Bez}}</math>
<!-------------------------------------------------------------------------------->
 
{{MyCodeBlock|title=Make Equations of Motion Dimensionless
Damit lautet die dimensionslose Bewegungsgleichung
|text=Text
 
::<math>\displaystyle \frac{d^2}{d\tau^2}\Phi + \frac{{{2}^{\frac{5}{2}}}\cdot {{\pi }^{2}}\cdot \sqrt{12-{{2}^{\frac{7}{2}}}\cdot \mathrm{cos}\left( \Phi\right) }\cdot \mathrm{sin}\left( \Phi\right) -{{2}^{\frac{7}{2}}}\cdot {{\pi }^{2}}\cdot \mathrm{sin}\left( \Phi\right) }{\sqrt{12-{{2}^{\frac{7}{2}}}\cdot \mathrm{cos}\left( \Phi\right) }} = 0 \text{ mit } \tau = t/t_{Bez}</math>
 
[[Datei:Kw24-11.png|mini|Rückstellmoment der Feder|alternativtext=]]Den nichtlinearen (zweiten) Summand können wir als dimensionsloses "Rückstellmoment auffassen und über dem Drehwinkel auftragen:
 
Nullstellen des Rückstellmoments zeigen die Gleichgewichtslagen des Systems an. Die Steigung dieses Rückstellmoments in den Gleichgewichtslagen deutet dabei auf die Stabilität der linearisierten Bewegung (oder Bewegung im Kleinen) hin. So ist die Lösung ''Φ'' = 0 eine instabile Gleichgewichtslage, die beiden anderen sind stabil.
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/*********************************************************/
/* make dgl dimensionless employing the eigenfrequency  */
 
dgl[2] : subst([diff(φ(t),t,2) = 'diff(Phi,t,2)/t[B]^2,
                    φ(t)      =      Phi ] , dgl[1]);
dgl[2] : solve(subst(ref,dgl[2]), 'diff(Phi,t,2));
 
/* equilibrium angles für Φ */
plot2d(-rhs(dgl[2][1]), [Phi,-1,1],
                        [legend, "Rückstell'moment'"],
                        [xlabel, "Φ/1->"], [ylabel, "M/1->"]);
</syntaxhighlight>
</syntaxhighlight>
}}
}}


==tmp==
 
<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Numerical Solution
{{MyCodeBlock|title=Numerical Solution
|text=Text
|text=Für die numerische Lösung müssen wir die Bewegungsgleichung
 
::<math>\ddot{\Phi} = f(\Phi)</math>
 
umschreiben auf eine Differentialgleichung erster Ordnung:
 
::<math>\left(\begin{array}{c}\Phi\\\Omega\end{array}\right) = \left(\begin{array}{c}\Omega\\f(\Phi)\end{array}\right)</math>
 
Diese lösen wir mit einer fertigen Routine - hier einem Runge-Kutta-Verfahren - als Anfangswertproblem. Anders als bei linearisierten Bewegungsgleichungen hängt die Charakteristik der Lösung von den Anfangsbedingungen ab.
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/*********************************************************/
/* numerical solution of IVP */
fpprintprec:3$
lhs(dgl[2][1]) = float(expand(rhs(dgl[2][1])));
times : subst([t0 = 0, tmax = 6, dt = 0.01],
                    [t, t0, tmax, dt]);
dgl1stOrder : [Omega,float(expand(rhs(dgl[2][1])))];
stateVabs : [Phi,Omega];
initiVals : [[1,0],[-1,0],
            [0,0.001],[0,3],
            [0,8], [0,10]];
/* solution of IVP (ivs) */
for i: 1 thru length(initiVals) do
  ivs[i] : rk(dgl1stOrder, stateVabs, initiVals[i], times)$
</syntaxhighlight>
</syntaxhighlight>
}}
}}




==tmp==
 
<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Post-Processing
{{MyCodeBlock|title=Post-Processing
|text=Text
|text=
Das schauen wir uns am Verlauf von Lösungen im Zeitbereich an:
 
<table class="wikitable" style="background-color:white;">
<tr>
<td>a)[[Datei:K24-21.png|rahmenlos]]</td>
<td>b)[[Datei:Kw24-24.png|rahmenlos]]</td>
</tr><tr>
<td>c)[[Datei:Kw24-22.png|rahmenlos]]</td>
<td>d)[[Datei:Kw24-25.png|rahmenlos]]</td>
</tr><tr>
<td>e)[[Datei:Kw24-23.png|rahmenlos]]</td>
<td>f)[[Datei:Kw24-26.png|rahmenlos]]</td>
</tr>
</table>
 
Lösung f) fällt aus dem Raster - hier gibt es keine periodische Lösung mehr: Bei einer anderen Skalierung der Φ-Achse erkennen wir, dass der Stab AC um den Punkt A rotiert:[[Datei:Kw24-30.png|mini|Lösung f) im Zeitbereich.|alternativtext=|ohne]]
Oft erkennt man das Charakteristische einer Lösung im Phasendiagramm. Hier wird die Geschwindigkeit ''Ω'' über dem Drehwinkel ''Φ'' aufgetragen:
[[Datei:Kw24-31.png|mini|Phasendiagramme der Lösungen|alternativtext=|ohne]]
Periodische Lösungen - die oft gesucht sind - erkennt man an einem geschlossenen Verlauf. Die hellblaue Lösung gehört zu den Anfangswerten von f), bei dem der Drehwinkel immer größer wird - sich also keine periodische Lösung ergibt.
 
Wir sehen, dass die Lösung e) (schwarz) eine "gerade noch" periodische Lösung ist - etwas mehr Anfangsgeschwindigkeit im System und die Lösung sieht wie in f) aus. Diese Lösung nennt man auch Separatrix - sie separiert zwei Charakteristiken von Lösungen.
 
Analog sehen Sie in der Lösung c) (grün) eine Separatrix zwischen der periodischen Lösung der kleinen Schwingung um die Gleichgewichtslage und solchen periodischen Lösungen, bei denen der Stab die Horizontale durchbricht: [[Datei:Kw24-32.png|mini|Phasendiagramme|alternativtext=|ohne]]
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
 
/*********************************************************/
/* plot results */
for i: 1 thru length(initiVals) do
(T[i] : makelist(ivs[i][j][1],j,1,length(ivs[i])),
P[i] : makelist(ivs[i][j][2],j,1,length(ivs[i])),
O[i] : makelist(ivs[i][j][3],j,1,length(ivs[i])),
plot2d([discrete, T[i], P[i]], [y,-%pi,%pi],
      [title, sconcat("start @: ",string(initiVals[i]))],
                  [xlabel,"τ/1->"], [ylabel,"Φ/1->"]));
 
/* phase plot */
curves : makelist([discrete,P[i],O[i]],i,1,length(initiVals))$
plot2d(curves, [legend, false], [x,-2,10],
                  [ylabel,"Ω/1->"], [xlabel,"Φ/1->"]);
plot2d(curves, [legend, false], [x,-2,2], [y,-5,5],
                  [ylabel,"Ω/1->"], [xlabel,"Φ/1->"]);
plot2d(curves, [legend, false], [x,-2,2], [y,-5,5], [svg_file, "Kw24-20.svg"],
                  [ylabel,"Ω/1->"], [xlabel,"Φ/1->"]);
</syntaxhighlight>
</syntaxhighlight>
}}
}}
<table class="wikitable" style="background-color:white; float: left;  margin-right:14px;">
<tr><td>[[Datei:K24-21.png|rahmenlos]]
</td><td>[[Datei:Kw24-24.png|rahmenlos]]
</td></tr><tr><td>[[Datei:Kw24-22.png|rahmenlos]]</td><td>[[Datei:Kw24-25.png|rahmenlos]]</td></tr><tr><td>[[Datei:Kw24-23.png|rahmenlos]]</td><td>[[Datei:Kw24-26.png|rahmenlos]]</td></tr></table>
[[Datei:Kw24-30.png|mini|Lösung im Zeitbereich]]
[[Datei:Kw24-31.png|mini|Phasendiagramme]]
[[Datei:Kw24-32.png|mini|Phasendiagramme]]


<hr/>
<hr/>

Aktuelle Version vom 29. März 2021, 06:04 Uhr


Aufgabenstellung

Ein Mechanismus besteht aus einem starren Stab AC (Länge 1, Massenmoment J) und einer Feder BC (Länge 2 im entspannten Zustand, Steifigkeit k).

Die Körper sind drehbar gelagert, der Abstand AB ist 3.


Lageplan

Hier geht es um die verschiedenen Lösungstypen des Systems - in Abhängigkeit von den Anfangsbedingungen. Gesucht sind Lösungen als Anfangswertproblem, nichtlineare Schwingungen und Aussagen über die Stabilität der Bewegungen.

Lösung mit Maxima

Header

Gesucht ist die Lösung zu dem nichtlinearen Anfangswertproblem.


/*********************************************************/
/* MAXIMA script                                         */
/* version: wxMaxima 15.08.2                             */
/* author: Andreas Baumgart                              */
/* last updated: 2017-09-30                              */
/* ref: Kw24 (TM-C, Labor 6)                             */
/* description: finds the solution for                   */
/*              the nonlinear IVP                        */
/*********************************************************/




Declarations

Wir wählen die Längenverhältnisse

so dass die Feder für φ(t) = 45° entspannt ist.


/*********************************************************/
/* declare variational variables - see 6.3 Identifiers   */
declare("δW", alphabetic); /* total virtual work         */
declare( "φ", alphabetic); /* coordinate of rotatin      */
declare("δφ", alphabetic); /* variation of coordinate phi*/
declare("δw", alphabetic); /* variation of coordinate w  */

/* parameter */
assume(l[1]>0, k>0, J>0);
params: [l[2] = l[1], l[3] = sqrt(2)*l[1]];




Kinematics

Der Mechanismus hat genau einen Freiheitsgrad, den wir entweder durch die Federlängung w oder den Drehwinkel φ erfassen. Den Zusammen zwischen den beiden stellen wir über

her. Nun steht uns frei, ob wir diese Beziehung nach w oder φ auflösen. Wir wählen φ als Minimal-Koordinate und erhalten

Nach den Regeln der Variationsmethoden erhalten wir daraus:

.

/*********************************************************/
/**** kinematics ***/

/* coordinates */
coords : [φ, δφ];

kin : (l[1]*sin(φ(t)))^2+(l[3]-l[1]*cos(φ(t)))^2=(l[2]+w)^2;
kin : trigsimp(ratsimp(solve(subst(params, kin),w)[2]));
kin : [kin,
       δw = subst([epsilon=0],
                 diff(subst([φ(t)=φ(t)+epsilon*δφ],
                                  subst(kin,w)),epsilon))];




Equilibrium Conditions

Mit dem Prinzip der virtuellen Verrückungen schreiben wir

und erhalten die Bewegungsgleichung


/*********************************************************/
/* equilibrium conditions */
PvV : δW = -J*'diff(φ(t),t,2)*δφ - k*w*δw;

/* dgl for φ ist the coefficient of δφ in δW */
dgl[1] : expand(coeff(subst([δW = 0], subst(kin, -PvV)),δφ));




Equilibrium Points

Die Bewegungsgleichung machen wir dimensionslos, indem wir charakteristische Größen aus der linearisierten Bewegungsgleichung verwenden.

Die drei Gleichgewichts-Lagen bzgl. derer wir die Bewegungsgleichung linearisieren können, sind

.

Eine Taylor-Reihenentwicklung des Rückstellmoments bzgl. φ0,3 ,also

liefert die linearisierte Bewegungsgleichung

.

Dabei ist ω0 die Eigenkreisfreuquenz der linearisierten Bewegungsgleichung.


/*********************************************************/
/* Zeros of restoring function */
equilib : solve(subst(['diff(φ(t),t,2)=0],dgl[1]),φ(t));

/*********************************************************/
/* derive the eigenfrequency                             */
/* of linearized ODE at stable point of equilob.         */
LinODE : taylor (subst([φ(t)=psi],rhs(dgl[1]/J)), psi, %pi/4,1);
ref : omega[0]^2 = coeff(expand(LinODE),psi);
ref : [t[B] = 1/(subst(solve(ref,omega[0])[2],omega[0])/(2*%pi))];




Make Equations of Motion Dimensionless

Eine Bezugszeit für unsere nichtlineare Bewegungsgleichung kommt nun aus

Damit lautet die dimensionslose Bewegungsgleichung

Rückstellmoment der Feder

Den nichtlinearen (zweiten) Summand können wir als dimensionsloses "Rückstellmoment auffassen und über dem Drehwinkel auftragen:

Nullstellen des Rückstellmoments zeigen die Gleichgewichtslagen des Systems an. Die Steigung dieses Rückstellmoments in den Gleichgewichtslagen deutet dabei auf die Stabilität der linearisierten Bewegung (oder Bewegung im Kleinen) hin. So ist die Lösung Φ = 0 eine instabile Gleichgewichtslage, die beiden anderen sind stabil.


/*********************************************************/
/* make dgl dimensionless employing the eigenfrequency   */

dgl[2] : subst([diff(φ(t),t,2) = 'diff(Phi,t,2)/t[B]^2,
                     φ(t)      =      Phi ] , dgl[1]);
dgl[2] : solve(subst(ref,dgl[2]), 'diff(Phi,t,2));

/* equilibrium angles für Φ */
plot2d(-rhs(dgl[2][1]), [Phi,-1,1],
                        [legend, "Rückstell'moment'"],
                        [xlabel, "Φ/1->"], [ylabel, "M/1->"]);




Numerical Solution

Für die numerische Lösung müssen wir die Bewegungsgleichung

umschreiben auf eine Differentialgleichung erster Ordnung:

Diese lösen wir mit einer fertigen Routine - hier einem Runge-Kutta-Verfahren - als Anfangswertproblem. Anders als bei linearisierten Bewegungsgleichungen hängt die Charakteristik der Lösung von den Anfangsbedingungen ab.


/*********************************************************/
/* numerical solution of IVP */
fpprintprec:3$
lhs(dgl[2][1]) = float(expand(rhs(dgl[2][1])));
		
times : subst([t0 = 0, tmax = 6, dt = 0.01],
                    [t, t0, tmax, dt]);
dgl1stOrder : [Omega,float(expand(rhs(dgl[2][1])))];
stateVabs : [Phi,Omega];
initiVals : [[1,0],[-1,0],
             [0,0.001],[0,3],
             [0,8], [0,10]];
/* solution of IVP (ivs) */
for i: 1 thru length(initiVals) do
   ivs[i] : rk(dgl1stOrder, stateVabs, initiVals[i], times)$





Post-Processing

Das schauen wir uns am Verlauf von Lösungen im Zeitbereich an:

a) b)
c) d)
e) f)

Lösung f) fällt aus dem Raster - hier gibt es keine periodische Lösung mehr: Bei einer anderen Skalierung der Φ-Achse erkennen wir, dass der Stab AC um den Punkt A rotiert:

Lösung f) im Zeitbereich.

Oft erkennt man das Charakteristische einer Lösung im Phasendiagramm. Hier wird die Geschwindigkeit Ω über dem Drehwinkel Φ aufgetragen:

Phasendiagramme der Lösungen

Periodische Lösungen - die oft gesucht sind - erkennt man an einem geschlossenen Verlauf. Die hellblaue Lösung gehört zu den Anfangswerten von f), bei dem der Drehwinkel immer größer wird - sich also keine periodische Lösung ergibt.

Wir sehen, dass die Lösung e) (schwarz) eine "gerade noch" periodische Lösung ist - etwas mehr Anfangsgeschwindigkeit im System und die Lösung sieht wie in f) aus. Diese Lösung nennt man auch Separatrix - sie separiert zwei Charakteristiken von Lösungen.

Analog sehen Sie in der Lösung c) (grün) eine Separatrix zwischen der periodischen Lösung der kleinen Schwingung um die Gleichgewichtslage und solchen periodischen Lösungen, bei denen der Stab die Horizontale durchbricht:

Phasendiagramme

/*********************************************************/
/* plot results */
for i: 1 thru length(initiVals) do
	(T[i] : makelist(ivs[i][j][1],j,1,length(ivs[i])),
	 P[i] : makelist(ivs[i][j][2],j,1,length(ivs[i])),
	 O[i] : makelist(ivs[i][j][3],j,1,length(ivs[i])),
	 plot2d([discrete, T[i], P[i]], [y,-%pi,%pi],
	      [title, sconcat("start @: ",string(initiVals[i]))],
	                   [xlabel,"τ/1->"], [ylabel,"Φ/1->"]));

/* phase plot */
curves : makelist([discrete,P[i],O[i]],i,1,length(initiVals))$
plot2d(curves, [legend, false], [x,-2,10],
                   [ylabel,"Ω/1->"], [xlabel,"Φ/1->"]);
plot2d(curves, [legend, false], [x,-2,2], [y,-5,5],
                   [ylabel,"Ω/1->"], [xlabel,"Φ/1->"]);
plot2d(curves, [legend, false], [x,-2,2], [y,-5,5], [svg_file, "Kw24-20.svg"],
                   [ylabel,"Ω/1->"], [xlabel,"Φ/1->"]);





Links

  • ...

Literature

  • ...