Autor Beitrag
LePomme
Hält's aus hier
Beiträge: 5



BeitragVerfasst: Sa 02.07.05 23:42 
Servus!

Ich möchte ein Programm schreiben, das die harmonische Schwingung simuliert,
also die DGl s(t)=-D/m*a(t) numerisch löst.
Den Teil habe ich (mit umfangreicher Unterstützung) auch geschafft.

Der Code, der den nächsten Schritt der DGl (mit der Schrittweite h) berechnet, lautet:

ausblenden volle Höhe Delphi-Quelltext
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:
  procedure TOsz.nextStep(var s: Single; var v: Single);
  var
    Fa, Fb, Fc, Fd: Array[0..1of Single;
    omega, z1, z2: Single;
  begin
    omega := - D / m;
    z1 := s;
    z2 := v;

    //############ Klassisches Runge-Kutta-Verfahren ############
    //# w1 = f(v[i], t[i])                                      #
    //# w2 = f(v[i] + h * w1 / 2, t + h / 2)                    #
    //# w3 = f(v[i] + h * w2 / 2, t + h / 2)                    #
    //# w4 = f(v[i] + h * w3, t + h)                            #
    //# v[i+1] = v[i] + h / 6 * (w1 + (2 * w2) + (2 * w3) + w4) #
    //###########################################################

    Fa[0] := z2;
    Fa[1] := omega * z1;

    Fb[0] := z2 + (h / 2) * Fa[1];
    Fb[1] := omega * (z1 + (h / 2) * Fa[0]);

    Fc[0] := z2 + (h / 2) * Fb[1];
    Fc[1] := omega * (z1 + (h / 2) * Fa[0]);

    Fd[0] := z2 + h * Fc[1];
    Fd[1] := omega * (z1 + h * Fc[0]);

    s := s + (h / 6) * (Fa[0] + (2 * Fb[0]) + (2 * Fc[0]) + Fd[0]);
    v := v + (h / 6) * (Fa[1] + (2 * Fb[1]) + (2 * Fc[1]) + Fd[1]);
  end;


Jetzt wollte ich nicht nur ein, sondern mehrere Oszillatoren abhängig voneinander, also gekoppelt, schwingen lassen.
Die untenstehende Funktion wird von einem Timer aufgerufen und das Ergebnis ausgegeben (habe ich hier weggelassen).


ausblenden volle Höhe Delphi-Quelltext
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:
  for i := 0 to AnzahlOsz - 1 do
  begin
    BuffS[i] := OszS[i];  //OszS -> Auslenkung des i-ten Schwingers
    BuffV[i] := OszV[i];  //OszV -> Geschwindigkeit
  end;
  
  if AnzahlOsz > 1 then         //AnzahlOsz: Konstanter Wert für die Anzahl der Schwinger
  begin

    //Falls mehrere Osz. vorliegen, werden die "Federn" nicht mehr zum 0-Wert gespannt,
    //sondern zu den Nachbarn.

    //Für einen ersten Test hielt ich es für unsinnig mit dem Pythagoras rumzufuchteln, 
    //daher werden im Folgenden nur Differenzen berechnet.

    //Die Buffer (BuffS und BuffV) sollen verhindern, dass beim einmaligen Durchlaufen der
    //Reihe nicht schon direkt Werte bis nach Hinten weitergegeben werden
    //Beispiel: 
    // 1 0 0 0 0 0     <- Jeder soll den halben Wert des Vorgängers annehmen
    //Ohne Buffer:
    // 1 1/2 1/4 1/8 usw.   (direkt beim ersten Durchlauf werden alle Werte gebildet)

    //Mit Buffer:           (Ausfüllen in mehreren Schritten)
    // 1 1/2 0 0 0 0
    // 1 1/2 1/4 0 0 0
    // 1 1/2 1/4 1/8 0 0
    // usw ... 
    

    for i := 0 to AnzahlOsz - 1 do
    begin
      if (i = 0then         //Beim ersten Oszillator: berücksichtige nur den nächsten
      begin
        OszS[i] := OszS[i] - BuffS[i + 1];
        Osz.nextStep(OszS[i], OszV[i]);
        OszS[i] := OszS[i] + BuffS[i + 1];
      end;

      if (i > 0and (i < AnzahlOsz - 1then     //Alles dazwischen
      begin
        OszS[i] := 2 * OszS[i] - BuffS[i - 1] - BuffS[i + 1];
        Osz.nextStep(OszS[i], OszV[i]);
        OszS[i] := 0.5 * (OszS[i] + BuffS[i - 1] + BuffS[i + 1]);
      end;

      if (i = AnzahlOsz - 1then        //Nur den vorherigen berücksichtigen
      begin
        OszS[i] := OszS[i] - BuffS[i - 1];
        Osz.nextStep(OszS[i], OszV[i]);
        OszS[i] := OszS[i] + BuffS[i - 1];
      end;
    end;
  end
  else
    Osz.nextStep(OszS[0], OszV[0]);      //Falls nur ein Schwinger vorliegt


Mein Problem ist, dass alles (bei mehr als einem Schwinger) ausser Kontrolle zu geraten scheint.
Zunächst läuft alles normal, aber nach einer Weile verschieben sich alles OszS-Werte immer weiter in den
negativen Bereich.

An der nextStep-Methode liegt es meines Erachtens nach nicht, denn ein einzelner Oszillator macht
nicht so nen Quatsch.

Habt ihr eine Idee?

Und Danke schonmal im Vorraus!
delfiphan
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic star
Beiträge: 2684
Erhaltene Danke: 32



BeitragVerfasst: So 03.07.05 00:22 
Erst mal hi und :welcome: im DF!

Du musst lediglich die Aufgabe mathematisch korrekt formulieren, dann erhältst du ein DGl-System. Dieses ist, wie das Problem selbst, mehrdimensional und gekoppelt; kann aber ohne weiteres mit Runge-Kutta gelöst werden. Deine ganze Geschichte mit dem "Die Buffer sollen verhindern, dass beim einmaligen Durchlaufen der Reihe nicht schon direkt Werte bis nach Hinten weitergegeben werden" schnall ich zwar nicht, ist aber auch nicht nötig. Denn, wenn du das physikalische Problem mathematisch korrekt formulierst, kannst du es mit Runge-Kutta lösen, ohne dich um die "Informationsverbreitung" kümmern zu müssen.

Übrigens: Zum Problem existiert (wenn ich mich nicht täusche) eine analytische Lösung. Du kannst das DGl-System nämlich entkoppeln, indem du ein Eigenwertproblem löst (dieses kann auch numerisch gelöst werden).
LePomme Threadstarter
Hält's aus hier
Beiträge: 5



BeitragVerfasst: So 03.07.05 00:37 
Danke für die schnelle Antwort :)

Wenn man die DGl analytisch löst erhält man eine Sinus-Funktion.
Aber das ist mir natürlich zu einfach ;)

Im Prinzip will ich versuchen nur anhand der Schwingungsgleichung eine Welle zu simulieren, was ja teilweise auch geklappt hat. Nur verschieben sich die Werte nach einer Weile - wie gesagt - immer weiter in den negativen Bereich.
delfiphan
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic star
Beiträge: 2684
Erhaltene Danke: 32



BeitragVerfasst: So 03.07.05 01:00 
Die analytische Lösung ergibt nicht "nur einen Sinus". Wenn du das gekoppelte System korrekt auflöst, kriegst du auch alle Effekte mit, die du in der Realität beobachtest. Also auch "Wellen". Wenn du mehrere gekoppelte Pendel hast, musst du das DGl-System auch entsprechend formulieren.

Beispiel mit 3 gekoppelten "Schwingern":

Seien im folgenden u1,u2,u3 die Auslenkungen von der Ruhelage. Dann gilt (Federkräfte summieren):
m1*ü1=k*(u2-u1)
m2*ü2=-k*(u2-u1)+k*(u3-u2)
m3*ü3=k*(u3-u2)

Das ist jetzt auch "ohne Pythagoras", gilt also nur für kleine Auslenkungen. Die 3 Gleichungen lassen sich komfortabel in Matrix-/Vektornotation schreiben:
M*ü = k*A*u

Dabei ist:
    /  1  -1   0 \
A = | -1   2  -1 |
    \  0  -1   1 /
    / u1 \
u = | u2 |
    \ u3 /
und M die Diagonalmatrix mit den Einträgen m1,m2,m3.

Da das Problem linear ist, sollte es auch ein Koordinatensystem geben, "aus dessen Sicht" die Gleichungen entkoppelt sind (man also also die Gleichungen einzeln und unabhängig voneinander lösen kann). Eine Koordinatentransformation geschieht mittels Multiplikation mit einer Transformationsmatrix T. Man setzt also u=T*y, wobei T eine noch unbekannte Matrix ist:
M*T*ÿ = k*A*T*y
bzw.
M*ÿ = k*T^-1*A*T*y

Das System ist genau dann entkoppelt, wenn die rechte Seite die Variablen nicht vermischt, also, wenn "k*T^-1*A*T" eine Diagonalmatrix ist.
Man muss also eine Matrix T und eine Diagonlamatrix D finden, sodass gerade gilt: T^-1*A*T = D bzw. A*T = T*D.
Etwas umformuliert: Man löse das Eigenwertproblem A*t=t*d (t = Vektor, d = Skalar). Die 3 Eigenvektoren und Eigenwerte schreibst du in die Matrizen T und D und erhältst:
    / -1  1  1 \
T = |  0  1 -2 |
    \  1  1  1 /
    /  1  0  0 \
D = |  0  0  0 |
    \  0  0  3 /

Das Problem ist jetzt auf M*ÿ = k*D*y reduziert, lautet also:
m1*ÿ1=1*k*y1
m2*ÿ2=0*k*y2
m3*ÿ3=3*k*y3

Das sind jetzt 3 voneinander völlig unabhängige Differentialgleichungen. Diese kannst du von Hand lösen, da sie jetzt entkoppelt sind (Sinus/Cosinus als Lösung bzw. eine lineare Funktion für 2. Gleichung). Wenn du die Lösung y1, y2, y3 hast, kannst du diese dann zurücktransformieren, um u zu erhalten:
u = T*y

Also:
u1=-y1+y2+y3
u2=y2-2*y3
u3=y1+y2+y3

Das wäre die analytische Lösung.

Falls du das ganze numerisch lösen willst, geht's aber wie gesagt ganz direkt mit Runge-Kutta. Einfach die Ausgangsgleichung M*ü = k*A*u einsetzen und fertig. Wo Zahlen waren, sind neu Vektoren. Die Formulierung bleibt aber gleich.

PS: Sorry für eventuelle Fehler
LePomme Threadstarter
Hält's aus hier
Beiträge: 5



BeitragVerfasst: So 03.07.05 12:28 
Wow!

Danke für die ausführliche Lösung! :D

Nur muss ich mir ein bisschen davon noch erarbeiten - wir sind noch net so weit im MatheLk :)
delfiphan
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic star
Beiträge: 2684
Erhaltene Danke: 32



BeitragVerfasst: So 03.07.05 14:59 
Wie gesagt, wenn du es numerisch lösen willst, brauchst du dich nicht um Koordinatentransformationen zu kümmern. Du musst lediglich dein Runge-Kutta in n-Dimensionen umschreiben. Im Forum geistert auch irgendwo die Funktion RKM4 (Runge-Kutta 4. Ordnung) herum. Die macht genau das..
Die mehr-dimensionale Funktion f (Beispiel für 4 Schwinger) muss dann folgendes zurückgeben:

(k*(u2-u1))/m1
(-k*(u2-u1)+k*(u3-u2))/m2
(-k*(u3-u2)+k*(u4-u3))/m3
(-k*(u4-u3))/m4
(keine Garantie auf Korrektheit. Hab das jetzt nur mal so hingeschrieben. Die Vorzeichen sind evtl. verkehrt rum)

Die Funktion f kannst du eigentlich genau so in Delphi implementieren...

Sieht vielleicht alles etwas kompliziert aus, aber im Grunde alles analog zum 1D-Fall.

In diesem Thread ging es um etwas ganz Ähnliches (dort findest du auch die Funktion RKM4).

//Edit: Sorry, die Funktion f für die numerische Lösung ist etwas komplizierter, denn die DGln sind von 2. Ordnung. 1 Gleichung 2. Ordnung <=> 2 Gleichungen 1. Ordnung. Mit Runge-Kutta kannst du nur DGLn 1. Ordnung lösen; also wird die Funktion f für 4 Schwinger 8-Dimensional. (Aber die 4 zusätzlichen Gleichungen sind mehr oder weniger Dummy-Einträge für die Ableitungen. Wenn du willst kann ich dir das auch erklären... Im Grunde nur eine kleine Substitution...)
LePomme Threadstarter
Hält's aus hier
Beiträge: 5



BeitragVerfasst: Mo 11.07.05 19:15 
Ich hab mittlerweile den Fehler gefunden!! :D

Ich hab das Programm einfach mit einer zu großen Schrittweite laufen lassen (peinlich! ;) )
Unter 0,001 läuft das alles glatt.

Aber trotzdem vielen Dank für die Tips :)
delfiphan
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic star
Beiträge: 2684
Erhaltene Danke: 32



BeitragVerfasst: Mo 11.07.05 23:02 
Ich frage mich, ob deine Buffer-Methode auch physikalisch korrekte Ergebnisse liefert, oder ob's einfach okay aussieht ;)
LePomme Threadstarter
Hält's aus hier
Beiträge: 5



BeitragVerfasst: Di 09.08.05 23:22 
Also es sieht anständig aus und die Geschichte mit der Reflexion am offenen und geschlossenen Ende funktioniert auch wunderbar :)