WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

 
Gast > Registrieren       Autologin?   

Partner:




Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

Runge-Kutta (Dormand/Prince) auf System von DGLn anwenden

 

amtmann

Gast


Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 02.09.2017, 11:48     Titel: Runge-Kutta (Dormand/Prince) auf System von DGLn anwenden
  Antworten mit Zitat      
Hallo Leute,

ich habe folgende DGL gegeben (und entsprechende AB y = y' = ... = y'''' = 0)
Code:
b4*y''''+b3*y'''+b2*y''+b1*y'+b0*y = a0*D


dabei sind bi und a0 bzw. D Konstanten. Die DGL habe ich nun in ein System erster Ordnung überführt

Code:

y1' = y2
y2' = y3
y3' = y4
y4' = 1/b4*(D*a0-y4*b3-y3*b2-y2*b1-y1*b0)
 


mit y1(0) = y2(0) = y3(0) = y4(0) = 0

... darauf würde ich nun gerne das Dormand/Prince-Verfahren anwenden. Da ich die ode-Suite leider nicht verwenden kann verwende ich hierzu ein älteres Skript (von http://www.mathstools.com/section/main/dormand_prince_method):

Code:
function [y, stats]=dopri54c(funcion,t0,t1,y0,tol,rpar)

t=t0;
y=y0;
h=tol^(1/5)/4;
step=0;
nrej=0;
fcall=1;
a4=[44/45 -56/15 32/9]';
a5=[19372/6561 -25360/2187 64448/6561 -212/729]';
a6=[9017/3168 -355/33 46732/5247 49/176 -5103/18656]';
a7=[35/384 0 500/1113 125/192 -2187/6784 11/84]';
e=[71/57600 -1/40 -71/16695 71/1920 -17253/339200 22/525]';
k1=feval(funcion,t,y,rpar);

while t < t1
    if t+h > t1
        h=t1-t;
    end
   
    k2=feval(funcion,t+h/5,y+h*k1/5,rpar);
    k3=feval(funcion,t+3*h/10,y+h*(3*k1+9*k2)/40,rpar);
    k4=feval(funcion,t+4*h/5,y+h*(a4(1)*k1+a4(2)*k2+a4(3)*k3),rpar);
    k5=feval(funcion,t+8*h/9,y+h*(a5(1)*k1+a5(2)*k2+a5(3)*k3+...
    a5(4)*k4),rpar);
    k6=feval(funcion,t+h,y+h*(a6(1)*k1+a6(2)*k2+a6(3)*k3+a6(4)*k4+...
    a6(5)*k5),rpar);
    yt=y+h*(a7(1)*k1+a7(3)*k3+a7(4)*k4+a7(5)*k5+a7(6)*k6);
    k2=feval(funcion,t+h,yt,rpar);
    est=norm(h*(e(1)*k1+e(2)*k2+e(3)*k3+e(4)*k4+e(5)*k5+...
    e(6)*k6),inf);
    fcall=fcall+6;

    if est < tol
        t=t+h;
        k1=k2;
        step=step+1;
        y=yt;
    else
        nrej=nrej+1;
    end

    h=.9*min((tol/(est+eps))^(1/5),10)*h;

end

if nargout > 1
    stats=[step nrej fcall];
end


bzw. mit der (hier beispielhaften Funktion y' = y)

Code:
function[y] = fun(x, y, rpar)
u=x;
end


das dann aufgerufen wird (für y(0) = 1) mittels
Code:

>> [y stad] = dopri54c('fun', 0, 20, 1, 0.0001,0);
 


... soweit verstehe ich das ja alles wie das für eine "einzelne" DGL funktioniert; aber ich schaffe es nicht so wirklich das nun auf mein DGL-System anzuwenden.
Irgendwie denke ich hier etwas zu umständlich; oder ich weiß einfach nicht wie ich "fun()" entsprechend abändern müsste.

Vielleicht kann mir hier jemand "unter die Arme greifen"?

Vielen Dank im Voraus,
amtmann


Harald
Forum-Meister

Forum-Meister


Beiträge: 24.448
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 02.09.2017, 15:42     Titel:
  Antworten mit Zitat      
Hallo,

die DGLen hast du ja schon aufgestellt. Du brauchst eigentlich nur noch aus y1 bis y4 dann y(1) bis y(4) zu machen.

Grüße,
Harald
Private Nachricht senden Benutzer-Profile anzeigen
 
amtmann

Gast


Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 02.09.2017, 17:05     Titel:
  Antworten mit Zitat      
Hallo Harald,

Zitat:
die DGLen hast du ja schon aufgestellt. Du brauchst eigentlich nur noch aus y1 bis y4 dann y(1) bis y(4) zu machen.


danke für Deine Antwort - bin mittlerweile auch dahintergekommen. Ich hatte nur einen Knoten im Kopf bzgl. der Übertragung des RK-Verfahrens auf DGL-Systeme ...

Eines, was mir noch einfällt ist etwas anderer Natur. Die DGL die ich hier Löse gilt prinzipiell als "steif". Allerdings bin ich mir ob ihrem "Steifheitsmaß" nicht so sicher, da sowohl ode45 als auch ode23s die DGL gut lösen können; ode45 braucht zwar geringfügig mehrere Schritte je Zeitschritt (41) als ode23s (11), was mich etwas wundert. Ich dachte, dass da sicher der Faktor 100 oder so dazwischenliegt.

Nun denn, ich habe mich nun etwas auch mit den Rosenbrock-Methoden (wie z.B. ode23s) beschäftigt, da mit diese für mein Problem etwas effizienter erscheinen, und bin über https://www.mathworks.com/help/pdf_doc/otherdocs/ode_suite.pdf auf den Abschnitt 3.1. gestolpert, der die in ode23s "verbauten" Gleichungen beschreibt.

Allerdings bin ich auf dem Gebiet von impliziten RK/RO-Methoden ein richtiger "Newbie". Allerdings möchte ich mir auch einen ähnlichen Löser wie ode23s schreiben, da ich die Matlab-Suite leider in einem "externen" Programm nicht verwenden kann/darf.

Ist es viel aufwändiger, ein derartiges Verfahren zu programmieren, als z.B. das vorhin thematisierte DP-Verfahren?
Gibt es evtl. hierzu auch ähnliche Codesnippets, an denen ich mich orientieren kann bzw., die ich auf meine Bedürfnisse anpassen kann?

Vielen Dank,
amtmann
 
Harald
Forum-Meister

Forum-Meister


Beiträge: 24.448
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 02.09.2017, 18:15     Titel:
  Antworten mit Zitat      
Hallo,

es ist sehr ungewöhnlich, dass man die fertigen Löser nicht verwenden kann / darf. Darf ich nach dem Grund dafür fragen?

Grüße,
Harald
Private Nachricht senden Benutzer-Profile anzeigen
 
amtmann

Gast


Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 02.09.2017, 20:06     Titel:
  Antworten mit Zitat      
Hallo,

Zitat:
es ist sehr ungewöhnlich, dass man die fertigen Löser nicht verwenden kann / darf. Darf ich nach dem Grund dafür fragen?


Sicherlich. Ich muss das o.g. Problem leider in einer Basic-ähnlichen Skriptsprache gelöst bekommen, die selbst keinerlei fertige Algorithmen bereitstellt und nur if/else/schleifen bzw. als "höchsten" Datentyp eindimensionale Arrays kennt.

Bisher habe ich das so gemacht, dass ich die ursprüngliche DGL auf analytischem Weg je Zeitschritt gelöst habe, wobei ich unterstelle, dass sich während eines Zeitschrittes die Koeffizienten der DGL nicht ändern. - Diese Näherung war aus ingenieurmäßiger Sicht hinreichend genau und hat auch gut funktioniert, allerdings treten mit fortschreitender Zeit irgendwann in den Exponenten der e-Funktionen große Zahlen auf, die leider dann einen Overflow verursacht haben - kurz: das Lösungsgleichungssystem war irgendwann immer sehr schlecht skaliert.

Ich versuche zunächst noch, meine Eingangsdaten anders zu skalieren, dass ich z.B. nicht in Tagen sondern in Monaten oder Halbjahren rechne, sodass besagte e-Funktionen relativ kleine Exponenten haben (3 Jahre im Vergleich zu 1000 Tagen macht hier schon einen unterschied).


Den "alternativsten" Weg wollte ich nun beschreiten durch das selbstständige Programmieren eines Zeitintegrationsverfahrens, um dieses dann "inkrementell" zu verwenden: Die DGL wird nun wiederum in Zeitschritte geteilt und auf jeden Zeitschritt wird das Verfahren angewendet. Die Werte die das Verfahren am Ende für y, y', y'' ... ermittelt hat, dienen dann als Anfangswerte für den folgenden Zeitschritt usw.

Um mich nicht den Rosenbrockmethoden (die ich - da ich nur "einfacher Bauingenieur" bin - noch nicht durchschaut habe) stellen zu müssen, habe ich zunächst versucht, ob sich denn auch z.B. mit ode45 oder auch ode23 auf obige DGL zur Lösung anwenden ließe.
Zu meiner Verwunderung hat beides recht unkompliziert geklappt (ode23 auch mit weniger Schritten) obwohl ich der Meinung war, dass die DGL selbst eine steife DGL ist und ich definitiv zu einem impliziten Verfahren ausweichen muss. (Die exakte Lösung ist nämlich von der Form y = S+ Summe_i_N C_i*exp(-bk*t), wobei i = 1, ..., 5 und bk = 0.1, 1, 10, 100, 1000, und S eine Konstante ist.)

Es zeigt sich wohl, dass ode45 mehrere Funktionsaufrufe als ode23s benötigt, dennoch habe ich - beim Vergleich eines Zeitschrittes der Matlab-Solver [die sicher 1000fach besser als eine eigene Implementierung sind] - den Unterschied zwischen 11 (ode23s und auch ode23) und 41 (ode45) nicht allzu groß empfunden.
- Vor allem da ich in meiner bisherigen Lösung auch je Zeitschritt ständig die Nullstellen des charakteristischen Polynoms (Bairstow-Verfahren) ermitteln und eine Gauss-Elimination durchführen muss und das wohl auch annähernd gleich viele Rechenschritte bedeutet ...

Ganz kurz zusammengefasst: Der Grund liegt in einer Skriptsprache, die die Matlabfunktionalität nicht bereitstellt bzw. in die ich diese leider nicht direkt einbinden kann. (Ich hätte zwar noch die Möglichkeit LUA zu verwenden, allerdings weiß ich auch hier nicht, ob man in LUA Matlab-Funktionen "einfach" einbinden kann, ohne diese selbst schreiben zu müssen ...).

Viele Grüße,
amtmann
 
Harald
Forum-Meister

Forum-Meister


Beiträge: 24.448
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 03.09.2017, 12:15     Titel:
  Antworten mit Zitat      
Hallo,

sorry, dass ich weiter bohre, aber so ersoarst du dir vielleicht viel Arbeit:
Warum musst du diese Skript-Sprache bzw. LUA verwenden?
Falls es darum geht, den Algorithmus auf spezielle Hardware aufzuspielen, können vielleicht auch die Codegenerierungsmöglichkeiten von MATLAB weiterhelfen?

Grüße,
Harald
Private Nachricht senden Benutzer-Profile anzeigen
 
Jan S
Moderator

Moderator


Beiträge: 11.057
Anmeldedatum: 08.07.10
Wohnort: Heidelberg
Version: 2009a, 2016b
     Beitrag Verfasst am: 16.09.2017, 16:23     Titel:
  Antworten mit Zitat      
Hallo amtmann,

Eiegene ODE-Solver in einer beliebigen Sprache zu schreiben, ist kein Hexenwerk. Alelrdings ist es nicht trivial die Funktionen zu Testen und die Ergebnisse zu validieren. Eine numerisch korrekte Integration benötigt unbedingt auch eine Analyse der Sensitivitäten, also der Empfindlichkeit der berechneten Trajektorie von einer Variation der Startwerte und der Parameter. Wenn eine Trajektorie hier sehr empfindlich ist (also "instabil"), ist die Integration eher eine Art komplizierter Zufallszahlen-Generator.

Das ODE45 4 mal so viele Funktionsauswertungen für ein steigfes System benötigt wie für ODE23 ist durchaus plausibel, hängt aber auch von der verwendeten Toleranz ab. Das kann darüber hinaus aber auch starke Auswirkungen auf die Sensitivität der Lösung haben.

Zitat:
Ganz kurz zusammengefasst: Der Grund liegt in einer Skriptsprache, die die Matlabfunktionalität nicht bereitstellt

Da dies hier ein Matlab-Forum ist, würde ich den Thread in die Kategorie Off-Topic verschieben, falls die Diskussion sich nicht auf Matlab bezieht.

Gruß, Jan
Private Nachricht senden Benutzer-Profile anzeigen
 
Neues Thema eröffnen Neue Antwort erstellen



Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:

Du kannst Beiträge in dieses Forum schreiben.
Du kannst auf Beiträge in diesem Forum antworten.
Du kannst deine Beiträge in diesem Forum nicht bearbeiten.
Du kannst deine Beiträge in diesem Forum nicht löschen.
Du kannst an Umfragen in diesem Forum nicht mitmachen.
Du kannst Dateien in diesem Forum posten
Du kannst Dateien in diesem Forum herunterladen
.





 Impressum  | Nutzungsbedingungen  | Datenschutz | FAQ | goMatlab RSS Button RSS

Hosted by:


Copyright © 2007 - 2024 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.