Numerische Simulation eines einfachen Oszillators

Neue Frage »

hp32 Auf diesen Beitrag antworten »
Numerische Simulation eines einfachen Oszillators
Bei der numerischen Simulation eines einfachen Oszillators (Federpendel) ist die Simulation instabil, wenn alle neuen Zustände aus den momentanen Zuständen berechnet werden: (Variante A)
p_neu = p + dt * x * D;
x_neu = x + dt * p * mr;

Wird hingegen für x_neu gleich der neu berechnete Zustand p_neu verwendet, ist das ganze stabil: (Variante B)
p_neu = p + dt * x * D;
x_neu = x + dt * p_neu * mr;


Im folgenden Beispiel ist die Schwingungsdauer 2*pi*sqrt(D/m) = 4.44s bzw. 444 Zeitschritte.

In Variante A weicht die maximale Amplitude nach 100s bereits um 14% ab. Die Frequenz stimmt überhaupt nicht.
Bei Variante B bleibt alles stabil.

Ist die gewählte Anzahl Zeitschritte einfach zu gering? Mit 4444 Zeitschritten sieht es besser aus in Variante A.

Welches Verfahren würde sich besser für solche Aufgabenstellungen eignen, d.h. ohne Umstellung, Differenzierung usw. der einzelnen Gleichungen?

Gibt es Verfahren die automatisch den Zeitschritt anpassen?



Nachfolgend der verwendete Code:

code:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
<html>
<head>
</head>
<body>
<script type="text/javascript" >

var dt = 0.01; // Zeitschritt
var t_total = 100;
var n = t_total/dt;

var m = 2, mr = 1/m, D = -1;
var tp = 2 * Math.PI *Math.sqrt(-D/m);

console.log("Schwingungsdauer: "+tp+" entspricht "+(tp/dt)+ " Zeitschritten");

var x,p,x_max,nullstelle;

///////////////////////////////////////////////////////////////
// Variante A

// Anfangswerte
x = 1, p = 0, x_max = x, nullstelle=undefined;

for ( var i = 0; i < n; i++ )
{
    var p_neu = p + dt * x * D;
    var x_neu = x + dt * p * mr;

    if( x_neu * x <= 0 )
    {
       if ( nullstelle !== undefined) 
          console.log("Variante A: Abweichung Schwingungsdauer + " + tp/(2*(i-nullstelle)*dt));
       nullstelle = i;
    }

    x = x_neu; 
    p = p_neu;

    if ( x_max < x )
    {
       x_max = x;
       console.log("Variante A: Neues Maximum: x["+i+"]="+ x_max);
    }
}

///////////////////////////////////////////////////////////////
// Variante B

// Anfangswerte
x = 1, p = 0, x_max = x, nullstelle=undefined;

for ( var i = 0; i < n; i++ )
{
    var p_neu = p + dt * x * D;
    var x_neu = x + dt * p_neu * m; // hier wird im gleichen Zeitschritt der neue Wert von p für die Berechnung von x verwendet

    if( x_neu * x <= 0 )
    {
        if ( nullstelle !== undefined)
            console.log("Variante B: Abweichung Schwingungsdauer + " + tp/(2*(i-nullstelle)*dt));   
        nullstelle = i;
    }

    x = x_neu;
    p = p_neu;

    if ( x_max < x )
    {
        x_max = x;
        console.log("Variante B: Neues Maximum: x["+i+"]="+ x_max);
    }
}

</script>
</body>
Ehos Auf diesen Beitrag antworten »

Deine numerische Rechnung kann nicht funktionieren, weil deine beiden Gleichungen unrichtig sind, nämlich




Richtig muss es heißen




-----------------
Zur Herleitung der richtigen Gleichungen geht man von der Schwindungsgleichung des Federschwingers aus. Dabei ist M die Masse und D die Federkonstante. Diese Gleichung 2.Ordnung muss man in 2 Gleichungen 1.Ordnung umwandeln. Dazu führt man die Geschwindigkeit ein, welche die 1.Ableitung des Ortes x nach der Zeit t ist. Nochmaliges Differenzieren nach der Zeit liefert die Beschleunigung . Setzt man beides in die ursprüngliche Gleichung ein, hat man wie gewünscht 2 Gleichungen 1.Ordnung




Obwohl es nicht notwendig ist, hast du die Geschwindigkeit v durch den Impuls ausgedrückt, also bzw. differenziert . Einsetzen liefert




Darin ersetzen wir die Differenzialquotienten durch die Differenzenquotienten





Umstellen nach bzw. liefert die beiden richtigen Gleichungen, die ich anfangs angegeben habe.

Zur Probe deiner numerischen Rechnungen benutze die formelmäßige Lösung , wobei sind C, gewisse Konstanten sind, die von der Anfangsgeschwindigkeit und dem Anfangsort abhängen. Da die Kreisfrequenz der Sinuskurve lautet , ist der (zeitliche) Abstand zweier Nullstellen (also die halbe Periodendauer) gerade . Das kann man als bequemes Kriterium dafür nutzen, ob die Schrittweise bei deiner Numerik kurz genug ist.
hp32 Auf diesen Beitrag antworten »

Besten Dank für deine ausführliche Antwor! Ich bin sicher, dass die Gleichungen stimmen, da das Ergebnis in Variante B mit der analytischen Lösung übereinstimmt.

Ich habe nicht erwähnt (es ist nur im JavaScript Code zu sehen) , dass ich D = -1 und mr = 1/m gesetzt habe.

Auf jeden Fall bin ich jetzt sicher, dass die Gleichungen stimmen, da du auf dieselben kommst.
Ehos Auf diesen Beitrag antworten »

Wenn du in deiner Rechnung mit einer negativen Federkonstante D rechnest und anstelle der "wahren" Masse M deren Kehrwert 1/M benutzt, dann stimmen deine Gleichungen tatsächlich mit meinen Gleichungen überein.

Folglich hast du irgendeinen Programmierfehler gemacht. Den solltest du aber selber finden.
hp32 Auf diesen Beitrag antworten »

Das Programm an sich ist fehlerfrei, d.h. es macht das, was ich es sollte. Nur das Ergebnis ist aus numerischer Sicht jedoch unbefriediegend in Variante A.
Neue Frage »
Antworten »



Verwandte Themen

Die Beliebtesten »
Die Größten »
Die Neuesten »