Schrittweitensteuerung RKV < Matlab < Mathe-Software < Mathe < Vorhilfe
|
Zur Zeit schlage ich mich gerade mit einer Schrittweitensteuerung beim klassischen RKV rum. Hat jemand vielleicht sowas schon mal in Matlab für DGL-Systeme gemacht mit nem klassischen RKV 4. Ordnung?
Das Klassische RKV für einen Iterationsschritt hab ich wie folgt implementiert
| 1: |
| | 2: | function [t y] = rk4_system_schritt(f, t_0, h, y_0, i)
| | 3: |
| | 4: | %% RK4_SYSTEM_SCHRITT macht einen expliziten Runge-Kutta-Schritt 4.Ordnung
| | 5: | %
| | 6: | % f = handle auf die DGL
| | 7: | % t_0 = Startpunkt
| | 8: | % h = aktuelle Schrittweite
| | 9: | % y_0 = Startwerte
| | 10: | % i = Iterationsschritt
| | 11: |
| | 12: | m = size(y_0,1);
| | 13: |
| | 14: | if m == 1
| | 15: | y_0 = y_0';
| | 16: | end
| | 17: |
| | 18: | k1 = h*f(t_0, y_0);
| | 19: | k2 = h*f(t_0+h/2, y_0(:)+0.5*k1);
| | 20: | k3 = h*f(t_0+h/2, y_0(:)+0.5*k2);
| | 21: | k4 = h*f(t_0+h, y_0(:)+k3);
| | 22: | y = y_0(:) + (k1 + 2*k2 + 2*k3 + k4)/6;
| | 23: | t = t_0 + h;
| | 24: |
| | 25: | end
|
Wie kann ich zu einen vernüftigen Schrittweitensteuerung bei einem DGL-System kommen?
Schwierigkeiten macht mir dabei die Fehlerabschätzung. Bei einer einzelnen DGL würde ich wohl den Betrag des vorherigen Ergebnisses vom Betrag des aktuellen Eergebnisses abziehen und den Wert anschauen. Sollte dieser Wert zu groß sein würde ich die Schrittweite wohl halbieren und das ganze wieder versuchen. In umgekehrter Richtung die Schrittweite vergrößern, wenn der Wert sehr sehr klein ist.
Gruß Jojo
Ich habe diese Frage auch in folgenden Foren auf anderen Internetseiten gestellt:
![[]](/images/popup.gif) HIER
An die angegebene Literatur komme ich jedoch leider nicht und aus der Erklärung werde ich nicht so richtig schlüssig.
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig  | | Datum: | 15:09 Fr 10.07.2009 | | Autor: | jojo2502 |
Also ich habe soweit wohl eine Lösung gefunden, die ich dem Forum nicht vorenthalten möchte.
Zuerst der Code für einen RK-Schritt:
| 1: | function [t y] = rk4_system_schritt(f, t_0, h, y_0)
| | 2: | %% RK4_SYSTEM_SCHRITT macht einen expliziten Runge-Kutta-Schritt 4.Ordnung
| | 3: | %
| | 4: | % f = handle auf die DGL
| | 5: | % t_0 = Startpunkt
| | 6: | % h = aktuelle Schrittweite
| | 7: | % y_0 = Startwerte
| | 8: | % i = Iterationsschritt
| | 9: |
| | 10: | m = size(y_0,1);
| | 11: |
| | 12: | if m == 1
| | 13: | y_0 = y_0';
| | 14: | end
| | 15: |
| | 16: | k1 = h*f(t_0, y_0);
| | 17: | k2 = h*f(t_0+h/2, y_0(:)+0.5*k1);
| | 18: | k3 = h*f(t_0+h/2, y_0(:)+0.5*k2);
| | 19: | k4 = h*f(t_0+h, y_0(:)+k3);
| | 20: | y = y_0(:) + (k1 + 2*k2 + 2*k3 + k4)/6;
| | 21: | t = t_0 + h;
| | 22: |
| | 23: | end
|
und dann die Simulation
| 1: | function LinSim(y_a, ta, te, h_0, epsi, maxPos)
| | 2: |
| | 3: | % y_a = Startvektor
| | 4: | % ta = Startzeit
| | 5: | % te = Endzeit
| | 6: | % h_0 = Startschrittweite
| | 7: | % epsi = Genauigkeit
| | 8: |
| | 9: |
| | 10: | %initial conditions
| | 11: | T = ta;
| | 12: | Y_1 = y_a(1,1);
| | 13: | Y_2 = y_a(2,1);
| | 14: | t_akt = ta;
| | 15: | y_akt1 = y_a(1,1);
| | 16: | y_akt2 = y_a(2,1);
| | 17: | FS = 0;
| | 18: | epsi_stern = epsi / 50;
| | 19: | h = h_0;
| | 20: | MINMAX_H = [h_0, h_0];
| | 21: | i = 1;
| | 22: |
| | 23: | while(0<1)
| | 24: | [t_h, y_h] = rk4_system_schritt(@DGL,t_akt, h, [y_akt1;y_akt2]);
| | 25: | [t_hh, y_hh] = rk4_system_schritt(@DGL,t_akt, 0.5*h, [y_akt1;y_akt2]);
| | 26: | [t_hh, y_hh] = rk4_system_schritt(@DGL,t_hh, 0.5*h, y_hh);
| | 27: | fehl = (y_hh(1) - y_h(1))/15;
| | 28: |
| | 29: | % Schrittweitenkontrolle(darf nicht kleiner als 10^-13 sein)
| | 30: | if((h < 1e-13) && T(i) < te)
| | 31: | h = 1e-13;
| | 32: | % warning(i, 'Schrittweite < 1e-13');
| | 33: | end
| | 34: |
| | 35: | % erste und zweite Bedingung zusammengefasst
| | 36: | if(abs(fehl) <= h*epsi)
| | 37: | % Ende erreicht?
| | 38: | if(t_h < te)
| | 39: | % Eintragen der relevanten Werte
| | 40: | T = [T, t_h];
| | 41: | Y_1 = [Y_1, y_h(1)];
| | 42: | Y_2 = [Y_2, y_h(2)];
| | 43: | FS = [FS, fehl];
| | 44: | t_akt = t_h;
| | 45: | y_akt1 = y_h(1);
| | 46: | y_akt2 = y_h(2);
| | 47: | if(h < MINMAX_H(1))
| | 48: | MINMAX_H(1) = h;
| | 49: | elseif(h > MINMAX_H(2))
| | 50: | MINMAX_H(2) = h;
| | 51: | end
| | 52: | % Zusatz für zweite Bedingung
| | 53: | if(abs(fehl) < epsi_stern)
| | 54: | h = 2 * h;
| | 55: | end
| | 56: | i = i + 1;
| | 57: | % Abbruchbedingung
| | 58: | else
| | 59: | break
| | 60: | end
| | 61: | % dritte Bedingung
| | 62: | else
| | 63: | h = h/2;
| | 64: | end
| | 65: | end
| | 66: |
| | 67: | % letzter Schritt
| | 68: | last_h = abs(te - T(i));
| | 69: | [t_h, y_h] = rk4_system_schritt(@DGL,t_akt, last_h, [y_akt1;y_akt2]);
| | 70: | [t_hh, y_hh] = rk4_system_schritt(@DGL,t_akt, 0.5*last_h, [y_akt1;y_akt2]);
| | 71: | [t_hh, y_hh] = rk4_system_schritt(@DGL,t_hh, 0.5*last_h, y_hh);
| | 72: |
| | 73: | fehl = (y_hh(1) - y_h(1))/15;
| | 74: |
| | 75: | T = [T, t_h];
| | 76: | Y_1 = [Y_1, y_h(1)];
| | 77: | Y_2 = [Y_2, y_h(2)];
| | 78: | FS = [FS, fehl];
|
Das ganze gilt für ein DGL System...
@DGL ist ein handel auf ein Funktion..
Greetz jojo
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 10:20 So 12.07.2009 | | Autor: | matux |
$MATUXTEXT(ueberfaellige_frage)
|
|
|
|
