4. ODE unabhängig von t lösen(Problem mit RK4 + ABM Verfahren)

ODE unabhängig von t lösen(Problem mit RK4 + ABM Verfahren)

  • B

    Beim Runge-Kutta und anderen Lösungsverfahren wird ja meistens mit einer Zeitabhänigkeit gerechnet um den nächsten Wert zum nächsten Zeitschritt zu ermitteln. ODE45 in Matlab und andere Lösungsverfahren machen das ja auch so.

    Das sieht ja dann ungefähr so aus:

    k1 = f(tn,yn)
k2 = f(tn + h/2, yn + k1/h*2)
k3 = f(tn + h/2, yn + k2/h*2)
k4 = f(tn + h, yn + k3)

yn+1 = yn + h*(1/6 * k1 + 4/6 * k2 + 4/6 * k3 + 1/6 * k4)

    Das ganze Funktionier auch sehr gut für Funktionen, wo die Kraft von der Zeit abhängig ist und bekannt ist. z.B. Fn = F0 * sin(omega*tn), da kann ich ja den Zeitschrifft tn+h einsetzen und erhalte die richtige Kraft für die Berechnung.

    Konkret wäre das ein Feder-Masse-Dämpfersystem, welches durch die folgende ODE beschrieben wird.

    dy(1) = y(2);
dy(2) = (1/m) * (F(t) -yps*y(2)  -c*y(1)) ;

    Wie mach ich das jetzt bei Funktionen welche ich von außen mit einem Wert und einem Zeitschritt füttere? Ich bekomme quasi den Wert für Fn und h.

    Mein System wird nun von außen mit einer Kraft belastet, die sich unregelmäßig ändert, also der Zeitschritt h ist auch nicht konstant.

    dy(1) = y(2);
dy(2) = (1/m) * (F -yps*y(2)  -c*y(1)) ;

    Kann ich in diesem Fall die Zeit vernachlässigen und das RK4 so schreiben:

    k1 = f(tn,yn)
k2 = f(Fn, yn + k1/h*2)
k3 = f(Fn, yn + k2/h*2)
k4 = f(Fn, yn + k3)

yn+1 = yn + h*(1/6 * k1 + 4/6 * k2 + 4/6 * k3 + 1/6 * k4)

    Ich hab das Gefühl da fehlt dann irgendwas.

    Anschließen soll dann das genauere Ergebnis für die nächsten Schritte mittels Adams-Bashforth-Moulton ermittelt werden, dort brauch ich aber auch den "zukünftigen" Wert.

    // Adam Bashforth Step(Predictor)
w_ab = yn + h/24 * ( 55 * f(Fn, yn) - 19 * f(Fn-1,yn-1) + 37 * f(Fn-2, yn-2) - 9 * f(Fn-3, yn-3))

// Adam Moulton Step(Corrector)
yn+1 = yn + h/24 * (9 * f(Fn+1,w_ab) + 19 * f(Fn, yn) - 5 * f(Fn-1, yn-1) + f(Fn-2, yn-2))

    Hier ist die Krux bei "Fn+1" ich kenne den Wert, ich kenne auch h, aber in der Formulierung die ich überall finde für das Verfahren steht auch wieder (t(i)+h)

    Vielleicht hat ja jemand eine Idee wie ich das formulieren muss. 🙂

    0
  • ?

    Hallo,

    wenn ich das noch richtig im Kopf habe, dann kommen die Gewichte 1/6 und 4/6 gerade daher, dass du für deinen Zeitschritt ein Polynom wählst, für welches du dann deinen Schritt exakt ausführst. Ich gehe davon aus, dass dort auch eingeht, dass du gerade die Funktion an den Stellen tnt_n, tn+h/2t_n + h/2 und tn+ht_n+h auswertest.

    Wenn du jetzt also deine Funktion an anderen Stellen auswertest, würde ich vermuten, dass sich mindestens die Koeffizienten vor k1k_1 bis k4k_4 ändern.

    0
  • B

    Hallo!

    Danke für deine Antwort, aber weiter hilft mir das leider nicht. 😞

    Wie gesagt, ich habe zwar den exakten Zeitschritt, aber dieser geht quasi gar nicht in die Formel mit ein, da diese eigentlich unabhängig von t ist.

    Meinst du ich muss mich von dem "t+h" verabschieden und dafür etwas anders nehmen?

    0
  • ?

    Hallo,

    vielleicht habe ich dich falsch verstanden. Ich war davon ausgegangen, dass dir die Zeitschritte von "außen" vorgegeben werden. Dann würdest du mindestens die "quadratur Gewichte" im Runge-Kutta-Verfahren ändern müssen, damit das noch funktioniert.

    Vielleicht fasse ich einfach mal zusammen, was ich verstanden habe und wir sehen mal weiter.

    Du willst das ODE-System:

    dy_1dt(t_0)=y_2(t_0)\frac{\mathrm{d} y\_1}{\mathrm{d} t}(t\_0) = y\_2(t\_0)
    dy_2dt(t_0)=1m(F(t_0)−αy_2(t_0)−βy_1(t0))\frac{\mathrm{d} y\_2}{\mathrm{d} t}(t\_0) = \frac{1}{m} (F(t\_0) - \alpha y\_2(t\_0) - \beta y\_1(t_0))
    wobei m,α,βm, \alpha, \beta konstant sind, lösen.

    Das Problem ist jetzt, dass du F(t0)F(t_0) nur an diskreten Punkten \tilde t\_1, \dots \tilde t\_n kennst. Habe ich das soweit richtig verstanden?

    Dann wäre jetzt meine erste Frage, warum du den Wert von F(t0)F(t_0) für beliebige Werte von t0t_0 nicht aus den Werten von F(t~j)F(\tilde t_j) interpolierst und ein normales Runge-Kutta Verfahren verwendest.

    0
  • B

    Ja, das genau ist das Problem.

    Ich habe ein System, welches sich mit Hilfe der ODE beschreiben lässt. Dieses wird von außen, in einem anderen Program bearbeitet, ausgewertet. Mir stehen dann die Daten für Fn{ F }_{ n } und hn{h}_{n} des jeweiligen Berechnungsschrittes n zur Verfügung. Damit kann ich das ODE dann wieder lösen, Übergebe mein Ergebnis des ODE an das andere Programme und erhalte die Werte für n+1 usw.

    Das Runge Kutta soll auch nur das Startbrett für das Adambs-Bashforth-Moulton verfahren sein, welches ja die eigentlichen Berechnungen dann ermittelt. Dort wird aber auch an einer Stelle auf tn+1{t}_{n+1} zugegriffen.

    α,β und &m; sind dann noch große Matritzen eines Systems, aber es funktioniert schon im "Kleinen" nicht.

    Das mit dem Interpolieren hatte ich bereits versucht. Da hatte ich die Idee den Gradienten aus den letzten beiden Werten zu bilden, sodass ich dann für mein RK4:

    F_grad = F(i) - F(i-1);
    s1 = f(F(i),            y(:,i));
    s2 = f(F(i)+(F_grad)/2, y(:,i) + s1 * h /2);
    s3 = f(F(i)+(F_grad)/2, y(:,i) + s2 * h /2);
    s4 = f(F(i)+F_grad,     y(:,i) + s3 * h);

    y(:,i+1) = y(:,i) + (s1 + s2+s2 + s3+s3 + s4)  *h/ 6;

    machen könnte. Das klappt auch für einen ein Massen-Schwinger-Dämpfer-System ganz gut, wenn ich eine einfache Sinusschwingung teste. Jedoch nicht so wirklich in meinem späteren System oder komplizierteren Kraftverläufen. Ich glaube, dass schon kleine Unstimmigkeiten da große Auswirkungen haben, also würde ich gerne soviele Fehlerquellen wie möglich aussschließen.

    Was meinst du denn mit "normalem" RK-Verfahren? 🙂

    Danke schon mal für deine Hilfe. Ich steh da irgendwie echt auf dem Schlauch.

    0
  • ?

    Ich glaube, so langsam verstehe ich. Wie sieht denn Matlab-Funktion f aus?

    0
  • C

    Der Threadtitel macht keinen Sinn, du willst die zeitabhängige! ODE
    d/dt x = g(F(t),x) = f(t,x)
    lösen, wobei f, wenn ich das richtig verstanden habe, durch die stückweise konstante Kraft unstetig von der Zeit abhängt.

    Wenn du ein Mehrschrittverfahren anwenden willst, brauchst du wohl oder übel einen konstanten Zeitschritt mit h < h_n.

    Beim Runge-Kutta Verfahren kannst du als Zeitschritt natürlich h_n oder Teilschritte davon wählen und dann in dieser Zeit mit dem konstanten F rechnen.
    Was sinnvoll ist, kommt darauf an, wie unstetig f ist.

    0
  • B

    @Progchild:

    Die Funktion f ist das entsprechende ODE:

    dy_i(:,1)=y_i(:,2);
dy_i(:,2)= M_sys\(-Gamma_sys*y_i(:,2)-K_sys*y_i(:,1)+F_system);

    M_sys, Gamma_sys und K_sys sind Matrizen, F_system ist der Kraftvektor mit dem das System belastet wird. Dieser wird zu jedem Zeitpunkt i neu berechnet, für eben jenen Zeitpunkt.

    Das ganze Zeug funktioniert ja auch, es kommen Ergebnisse raus usw. die Frage ist halt, wie richtig sie sind.

    @C14:

    Hmmm, naja ich habe halt zum Berechnungszeitpunkt die nötigen Daten zur Verfügung. Ich habe eine Beispieldatei mit einem Kraftverlauf und den zugehörigen Zeitschritten. Damit "füttere" ich gerade von außen mein Programm. Aber ich hab halt kein klassisches F(t) was ich mir im Vorfeld berechnen könnte wie bei einer einfachen sinus-Schwingung.
    Daher kommen ja meine Probleme.

    So sieht ein Ausschnitt meines Kraftverlaufs aus:

    http://abload.de/img/lalat4u7b.jpg

    Im hinteren Bereich schwingt es sehr stark um 0,

    http://abload.de/img/lalavdka2.jpg

    Aber wie gesagt, wenn das Program dann im Einstz ist, habe ich diese Daten nicht sondern immer nur die aktuellen für den momentanen Berechnungsschritt, den Zeitschritt und alle alten Daten.

    Das Runge Kutta ist nicht mal so wichtig, da es quasi "nur" Startwerte für das ABM-Verfahren bereitstellt.

    http://abload.de/img/lalayvkb7.jpg

    Das Bild zeigt mal eine Rechnung, Türkis ist die Kraft, rot ist das ungedämpfte System und blau ist mit leichter Dämpfung.

    Die nächste Frage wäre, woher kommen die ganzen Schwingungen. 😃

    0
  • ?

    Hallo,

    ich war davon ausgegangen, dass du

    dydt=B(t,y)\frac{\mathrm{d}y}{\mathrm{d}t} = B(t, y)
    mit
    y:R→Rdy: \mathbb{R} \to \mathbb{R}^d,
    B:R×Rd→RdB: \mathbb{R} \times \mathbb{R}^d \to \mathbb{R}^d
    lösen möchtest.

    Wobei

    B(t0,y)=[y_2(t_0)1m(F(t_0)−αy_2(t_0)−βy_1(t0))]B(t_0, y) = \begin{bmatrix} y\_2(t\_0) \\ \frac{1}{m} (F(t\_0) - \alpha y\_2(t\_0) - \beta y\_1(t_0)) \end{bmatrix}

    .

    Du schreibst, dass du die Zeitschritte in ein externes Programm weiter gibst, um neue Werte für FF zu bekommen. Das bedeutet doch aber, dass dein FF vielleicht nicht von tt, aber sehr wohl von yy abhängt. (Habe ich das richtig verstanden?) Die rechte Seite deiner Gleichung ist also

    B(t0,y)=[y_2(t_0)1m(F(t_0,y(t_0))−αy_2(t_0)−βy_1(t_0))]B(t_0, y) = \begin{bmatrix} y\_2(t\_0) \\ \frac{1}{m} (F(t\_0, \mathbf{y(t\_0)} ) - \alpha y\_2(t\_0) - \beta y\_1(t\_0)) \end{bmatrix}

    .

    Nennen wir mal dein externes Programm GG. Dann liefert dir
    G(t_0,y(t_0))≈F(t_0+h~,y(t_0+h~))G(t\_0, y(t\_0) ) \approx F(t\_0 + \tilde h, y(t\_0+ \tilde h)),
    wobei h~\tilde h die von deinem externen Programm ermittelte Schrittweite ist.

    Damit ist

    B(t_0, y) \approx \begin{bmatrix} y\_2(t\_0) \\ \frac{1}{m} (G(t\_0 - \tilde h, y(t\_0 - \tilde h) ) \- \alpha y\_2(t\_0) - \beta y\_1(t\_0)) \end{bmatrix}

    .

    Angnommen, dass die Approximation GG gut genug ist, kannst du dies nun in z.B. das RK4 Verfahren (analog jeden anderen ODE-Löser) einsetzen:

    yn+1=y_n+h6(k_1+2k_2+2k_3+k4)y_{n+1} = y\_n + \tfrac{h}{6}(k\_1 + 2 k\_2 + 2 k\_3 + k_4)
    tn+1=tn+ht_{n+1} = t_n + h

    k_1=B(t_n,yn)k\_1 = B(t\_n, y_n)
    k_2=B(t_n+h2,y_n+h2k_1)k\_2 = B(t\_n + \tfrac{h}{2}, y\_n + \tfrac{h}{2} k\_1)
    k_3=B(t_n+h2,y_n+h2k_2)k\_3 = B(t\_n + \tfrac{h}{2}, y\_n + \tfrac{h}{2} k\_2)
    k_4=B(t_n+h2,y_n+hk_3)k\_4 = B(t\_n + \tfrac{h}{2}, y\_n + h k\_3)

    In diesem Fall siehst du, dass du BB an den Stellen tnt_n, tn+h2t_n + \tfrac{h}{2}, tn+ht_n + h benötigst. Das bedeutet, dass du dein Verfahren entweder zwingen musst, zweimal hintereinander den gleichen Zeitschritt auszuführen oder du musst irgendwie interpolieren.

    Falls dein Problem gutartig genug ist, dann könnte das funktionieren. Andernfalls solltest du überlegen, ob du die beiden Löser eher über eine Fixpunkt-Iteration o.Ä. koppelst.

    0
Anmelden zum Antworten