Lösen von DGL in C++

sothis_

aus der "dünnen" gsl dokumentation mit ein paar kommentaren meinerseits :):

#include <stdio.h>
     #include <gsl/gsl_math.h>
     #include <gsl/gsl_deriv.h>

// das hier ist die funktion die numerisch differenziert werden soll
// in dem fall f(x) = x^{3/2}, die parameterliste bleibt, der rückgabe-
// wert kommt in deinem fall natürlich von einer trigonometrischen funktion ;)

     double f (double x, void * params)
     {
       return pow (x, 1.5);
     }

     int
     main (void)
     {
       gsl_function F;
       double result, abserr;
// der funktionszeiger der obigen funktion
       F.function = &f;
       F.params = 0;

       printf ("f(x) = x^(3/2)\n");
// anstieg an stelle 2.0     
       gsl_deriv_central (&F, 2.0, 1e-8, &result, &abserr);
       printf ("x = 2.0\n");
       printf ("f'(x) = %.10f +/- %.10f\n", result, abserr);
       printf ("exact = %.10f\n\n", 1.5 * sqrt(2.0));
// anstieg an stelle 0.0      
       gsl_deriv_forward (&F, 0.0, 1e-8, &result, &abserr);
       printf ("x = 0.0\n");
       printf ("f'(x) = %.10f +/- %.10f\n", result, abserr);
       printf ("exact = %.10f\n", 0.0);

       return 0;
     }

calli82

wow, vielen dank.
wo hast du die doku mit dem beispiel gefunden? mir liegt nur eine doku vor, wo die nackten funktionen beschrieben werden und wo ich absolut nix mit anfangen kann...
wenn die doku so ist, ist sie ja sogar brauchbar

lieben dank für deine hilfe!

sothis_

calli82 schrieb:

wo hast du die doku mit dem beispiel gefunden?

http://www.gnu.org/software/gsl/manual/html_node/

im grunde exisitiert da sogar für jede noch so kleine einzelne funktion ein programm beispiel

calli82

nur mal so als mathematische Frage (bin da wirklich absoluter Anfänger auf diesem Gebiet, ein bisschen schmunzeln ist also erlaubt):

wenn ich jetzt das System von DGL auflöse, dann müsste ich doch das Integral bilden der beiden rechten Seiten und dann habe ich einen speziellen skalaren Wert heraus, der (in dem Fall die x,y - Position eines Objektes im 2-dimensionalen Raum) beschreibt, oder? Wieso benutzt du in diesem Beispiel eine Ableitung zur Lösung?
und welche rolle spielen hierfür z.B. die verwendeten "Einheiten" der Variablen, beispielsweise die Geschwindigkeit (Pixel pro Sekunde, ...)

Kennt jemand ein "Mathe fitt" Programm zum Erlernen höherer Mathematik, wo man sich viele Dinge am besten nochmal anschaulich ansehen und lernen kann (bin zwar studierter Informatiker, aber mit nicht so dollen Mathekenntnissen)... am liebsten wäre mir was multimediales...

sothis_

calli82 schrieb:

wenn ich jetzt das System von DGL auflöse, dann müsste ich doch das Integral bilden der beiden rechten Seiten und dann habe ich einen speziellen skalaren Wert heraus, der (in dem Fall die x,y - Position eines Objektes im 2-dimensionalen Raum) beschreibt, oder?

da frag mal lieber einen anderen eierkopf, der sich mit dynamischen systemen auskennt

calli82 schrieb:

Wieso benutzt du in diesem Beispiel eine Ableitung zur Lösung?

es war nur ein beispiel wie man eine funktion numerisch differenzieren kann und sollte als anstoß dienen wie man mit funktionen in gsl umgeht. eine lösung eines systems von differentialgleichungen ist das nicht

calli82

habe mal eine mathematisch methodische Frage, vielleicht kann die einer beantworten.

Gegeben sei eine DGL der Form
dy/dt = y*...
So, nun will ich sozusagen den Wert von y zu einem Zeitpunkt t ausrechnen.
Nach meinem mathematischen Kenntnisstand muss ich dann die rechte Seite integrieren. Was kommt dann raus, ein unbestimmtes Integral und eine Konstante c (unendlich viele Möglichkeiten), die dann zuaddiert wird.

Wie erreiche ich, dass ich y zu einem zeitpunkt t rauskriege.
Reicht das nicht, wenn ich die Größen auf der rechten Seite einfach einsetze und bekomme ich nicht für dy/dt schon einen plausiblen, richtigen Wert raus, ohne zu integrieren?

Besten Dank und sry für die dämichen Fragen.

Calli82

Theston

Wie willst du die Größen auf der rechten Seite einsetzen, wenn du eine unbekannte Größe y(t) hast?

Daniel E.

calli82 schrieb:

Gegeben sei eine DGL der Form
dy/dt = y*...
So, nun will ich sozusagen den Wert von y zu einem Zeitpunkt t ausrechnen.
Nach meinem mathematischen Kenntnisstand muss ich dann die rechte Seite integrieren. Was kommt dann raus, ein unbestimmtes Integral und eine Konstante c (unendlich viele Möglichkeiten), die dann zuaddiert wird.

Wie erreiche ich, dass ich y zu einem zeitpunkt t rauskriege.
Reicht das nicht, wenn ich die Größen auf der rechten Seite einfach einsetze und bekomme ich nicht für dy/dt schon einen plausiblen, richtigen Wert raus, ohne zu integrieren?

Du hast ein Anfangswertproblem, also einen Wert y(0) gegeben? Ja, dann kannst Du das so ähnlich machen. Aus y(0) folgt y'(0) und damit kannst Du dann y(x) für x nahe 0 ganz gut abschätzen. So etwas nennt sich Einschrittverfahren zur Lösung von DGLs: die einfachste Idee -- die Du im wesentlichen nacherfunden hast -- nennt sich Euler-Cauchy-Verfahren, kompliziertere Verfahren, die dafür die DGL auch besser lösen, findest Du zB. unter Heun oder Runge-Kutta.

calli82

ok, also ich bin sehr verwirrt und je länger ich drüber nachdenken, desto mehr verwirrt mich das =).

ich versuchs nochmal, unmathematisch, qualitativ zu formulieren.

Eine dgl der form dy/dt = y*... ist gegeben. diese drückt aus, dass eine größe y zu einem zeitpunkt t einen bestimmten wert annimmt. diesen muss ich nun für ein bestimmtes ti (meinetwegen alle 5 zeiteinheiten) in abhängigkeit der anderen einflussgrößen aktualisieren. Daraus kann ich dann eine Simulation bauen. dafür brauche ich das.
Die anderen Größen, die in die dgl einfließen, werden teilweise durch nicht-dynamische gleichungen für jeden zyklus errechnet.

das bedeutet, ich müsste bei der dgl also folgendes machen:
int(dy/dt) = euler_int(y*...) bzw. y(t) = euler_int(y*...).

Ich bekomme dann nach auflösen ein approximatives Ergebnis für y(t) heraus, welches mir dann einen skalaren wert der größe y liefert ja? nur warum gibt es dann kein "exaktes ergebnis". Die Gleichung kann ich ja nicht vorhersagen oder?

wie kann ich das mit gsl machen? hier sind funktionen für systeme von differentialgleichungen gegeben. thereotisch und praktisch kann ich das doch einfach die bibliothek machen lassen... oder?

lieben dank,

matthias

Taurin

Das klingt bei dir noch etwas durcheinander. Dich interessieren Gleichungen der Form

dy/dt = y * f(y,t),

wobei f(y,t) nicht von irgendwelchen Ableitungen abhängt. Beispiel:

dy/dt = y * t    (1)

Auf beiden Seiten integrieren klappt nicht, denn dann steht da

integral dy/dt dt = integral y * t dt

Mit der rechten Seite kann man nix anfangen. (1) umformen bringt uns weiter:

1/y * dy/dt = y * t
integral 1/y * dy/dt dt = integral t dt
integral 1/y dy = integral t dt
ln y = t^2 / 2 + C               (2)
y(t) = exp( t^2 / 2 + C) = K * exp( t^2 / 2 )   (3)

In (2) habe ich beide Integrationskonstanten zu einer zusammengefasst, damit es sich etwas einfacher schreibt; in (3) habe ich C durch K ersetzt, wieder als "schreibhilfe".

Ist ein Anfangswert gegeben, kann man daraus das K bestimmen. Z.b. y(0) = 2 -> K = 2.

Um analytische Lösungen von DGLs zu bekommen, gibt es eine ganze Reihe von Kochrezepten, die je nach der Gestalt der Gleichung (oder des Gleichungssystems) einsetzbar sind. Oft lassen sich die Lösungen aber gar nicht mehr analytisch aufschreiben. Bei den interessanten DGLs ist das sogar eher die Regel als die Ausnahme. Deswegen greift man zu numerischen Verfahren. Das einfachste ist das hier schon angesprochende Eulerverfahren:

Ich kenne y(0), daher kenne ich auch y'(0). Am Punkt (0,y(0)) nähert man die Funktion durch eine Gerade an und geht einen kleinen Zeitschritt delta_t weiter. Dann kennt man eine Näherung y(t_1) an der Stelle t_1 = 0 + delta_t. Damit bekomme ich eine Näherung y'(t_1) und kann wieder ein delta_t weitergehen. Mit jedem Schritt macht man einen kleinen, zusätzlichen Fehler. Mit jedem Schritt wird der Fehler größer.
Wenn du die Schrittweite halbierst, halbiert sich auch der Fehler.

Wenn es dich interessiert, kannst du dir mal das Skript http://www.tu-harburg.de/ins/lehre/material/numersim.pdf anschauen. Das ist aber keine leichte Kost, auch wenn es auf den ersten Blick schlimmer aussieht, als es ist. Vor allem der Anfang sollte noch ganz gut verständlich sein. Das Euler-Verfahren (siehe oben) findest Du auf Seite 23 noch mal anständig erklärt.

calli82

Hallo,

super. vielen dank für die hilfe und infos. hier bekommt echt superschnelle antwort und sehr viel hochwertigen support! Absolut klasse!

meine dynamischen gleichungen sind immer in der form
dy/dt = y... gegeben, sodass ich vorbereitend auf die Lösung der DGL mit GSL also erst einmal die ganzen Gleichungen so umformen muss, dass ich sie in der Form y(t) = exp(...) habe (Exponentialansatz). Dann kann ich erst sinnvoll damit rumrechnen!

Die mathematischen Hintergründe bleiben mir zwar immernoch etwas verschlossen (z.B. warum ist es notwendig, die Funktion z.B. wie bei Euler durch eine Gerade anzunähern und kann nicht einfach t in die umgeformte Gleichung y(t) einsetzen und erhalte dann den wert [grinsen über die dummheit der fragen erlaubt]), jedoch kann ich fürs erste damit eine Lösung bauen.

Also bleibt für mich folgende Reihenfolge:
1.) Dynamische Gleichungen vernünftig umformen.
2.) Anfangswerte berechnen.
3.) DGL nummerisch lösen lassen durch GSL (z.B. Runge-Kutta, Heun oder Euler), möglichst kleine delta_t.

Merci,

Calli

Taurin

Du schreibst immer dy/dt = y*.... Kannst Du mal genauer sagen, was .... ist?
Steht da eine Funktion, die nur von t abhängt? Oder kann da auch ein y stehen?
(Je genauer Du die Pünktchen kennst, desto besser)

Hast Du einen Anfwangswert y(t_0) = y_0?

Zu Euler: Das Näherung durch eine Gerade macht man gerade deswegen, weil man davon ausgeht, dass die rechte Seite der Gleichung ...
- ...unbekannt ist oder
- ...so hässlich ist, dass man nicht umformen kann.

Nur für den Fall, dass bei dir was durcheinander geht: Wenn Du eine analytische Lösung hast, brauchst Du keine Numerik mehr!

Mal zur Übung und zum besseren Verständnis, löse mal folgende Gleichungen (numerisch oder analytisch):

dy/dt = y*sin(t)
dy/dt = y - t * exp(y)

(jeweils mit dem Anfangswert y(0) = 1)

Das Eulerverfahren lässt sich auch in C leicht hinschreiben. Hilft dir bestimmt auch beim Verständnis.

calli82

ok, ich gebe jetzt mal ein Beispiel meiner dynamischen Gleichung:

d_omega/d_t = (alpha + 1/2*pi)*c_o*f_o + alpha*omega-gamma*omega^3
wie du siehst, hängt t von verschiedenen Faktoren ab.

Einen Anfangswert für omega habe ich, jedenfalls kann ich ihn aus einer bestimmten Konfiguration errechnen.

Also das mit der Näherrung verstehe ich nicht (gibt es ein Beispiel hierzu?). obere Funktion umzuformen ist für mich schon ziemlich fies, werds aber wohl machen müssen. Am Ende steht dann in etwa omega(t) = e^(...).

Mir ist an dieser Stelle einfach noch nicht klar, warum man Euler oder dergleichen anwenden muss. Wenn ich die obere Gleichung umgeformt habe, kann ich doch einfach ein t einsetzen, warum dann die Krücke?

Die von dir gegebenen Aufgaben werde ich in einer ruhigen Stunde bearbeiten. Danke hierfür.

Taurin

Wir kommen der Sache näher. Wenn deine rechte Seite unabhängig von t ist und alpha, f_o,c_o alle konstant sind, ist es einfach.

Dann steht da sowas wie

domega/dt = p(omega),

wobei p(omega) einfach ein Polynom mit konstanten Koeffizienten ist.
Dann erhält man nach teilen durch p(omega)

integral 1/p(omega) d_omega = integral 1 dt
integral 1/p(omega) d_omega = t + C

Die Linke Seite kann man integrieren und erhält einen Term, der nur von omega abhängt. ( Das vorzuführen ist mir jetzt zu mühsam, das kannst du selber machen ). Jetzt gibt es zwei Fälle:

a) Die Gleichung lässt sich nach omega umstellen. Dann haben wir gewonnen und es kommt etwas raus wie

omega(t) = ... irgendwas, was von C und t abhängt ...

b) Die Gleichung ist zu kompliziert und wir können nicht nach omega umstellen.
Wir haben aber trotzdem eine implizite Lösung gewonne. Für gegebens t kann man einen numerischen Gleichungslöser anschmeissen, der dazu das entsprechende omega ausrechnet.
Problem: Mit etwas Pech ist die Gleichung nicht eindeutig lösbar.

Im Fall b) kannst du auch (statt der Umformung) direkt einen numerischen Löser für die DGL benutzen. Aus dem Bauch heraus würde ich das vorziehen, da es sehr gute Verfahren dafür gibt, z.b. wird Runge-Kutta sehr oft benutzt.
(Zum Gleichungslösen gibt es auch gute Verfahren, aber dafür muss man ja zuerst einen kleinen Umweg gehen. Und Umwege sind doof.)

edit:

Warum brauche ich Euler und Co?

Weil man viele Gleichung *nicht* anständig umformen kann. Versuch dich mal an den Aufgaben von mir oben. Wenn ich mich recht ans Grundstudium erinnere, dürfte es bei einer von denen nicht klappen.

calli82

ja, exakt. Das Lösen über Möglichkeit (a) suche ich!

Um es nochmal zu verdeutlichen: ich habe einen vorgeschriebenen berechnungszyklus. alle hirsebrei zeiteinheiten errechnen sich die werte wie alpha und f_o, die in der Gleichung vorkommen, neu durch andere formeln.

Die Werte setze ich dann in die DGL-Formel ein und erhalte somit die "Konstanten" so wie du sie eben geschrieben hast, sodass ich irgendwann nach Auflösen omega(t) = irgendwas von t und c abhängiges stehen habe. dabei setz ich dann halt die ts im abstand von hirsebrei ein.

Ich möchte jetzt einfach nur durch ein kleines c++-programm mit gsl die entsprechenden gleichungen lösen lassen. Durch das Umstellen der Formel habe ich die Gleichung soweit vorbereitet, dass ich sie dem gsl via funktionspointer übergeben kann.

mein prof. schreibt mir in der aufgabenstellung einfach nur, dass "für das Lösen der dynamischen Gleichungen z.B. das Eulerverfahren benutzt werden kann". mehr nicht. ist ja egal, ich kann auch ein anderes verfahren nehmen.
das ist denke ich nun ein weg, den ich umsetzen kann (an soner aufgabe kann man viel lernen, wenn man masochisisch ist)

und genau da setzt jetzt meine Frage ein. Wenn ich doch sowieso omega(t) = irgendwas, was von c und t abhängt habe, warum gehe ich dann nicht einfach her und setze in die Gleichung für t meine zeiten ein und erhalte direkt den funktionswert, sondern benutze die krücke eines approximativen verfahrens. das verstehe ich noch nicht.

edit:
habe deinen edit erst jetzt gesehen. m.a.w. brauche ich das im falle der oben genannten gleichung gar nicht dinge wie euler & co, sondern nur bei etwas komplexeren funktionen wie sin und cos...
dann brauche ich ja nichtmal gsl dafür...

Taurin

edit: Hier stand Quatsch.

calli82

Hallo,

habe nun was implementiert (Eulerverfahren) und noch eine kleine Rückfrage:

Für jeden Zeitpunkt t errechne ich ein paar Parameter, die in die DGL mit einfließen. Also ändern sich jedesmal die Konstanten, von denen meine Größe abhängt.
Nun muss ich aber nicht jedesmal nochmal komplett neu anfangen und alle vorherigen Schritte berechnen, sondern muss zu jedem Zeitpunkt t eigentlich nur den vorherigen Zustand wissen, damit die Gleichung gelöst werden kann, oder?

Viele Grüße,

calli

Daniel E.

Ich denke, man sollte nicht "Konstanten" sagen, wenn sie offensichtlich nicht konstant sind.

Aber ansonsten: ja, eine große Klasse von DGLs kann man so numerisch approximieren.

calli82

wunderbar.

Wenn ich mir das Eulerverfahren zur Lösung von DGL angucke, wird die Lösung ja durch einen Polygonzug approximiert mit der Berechnungsvorschrift

...
h = 1/N;
...
T[k+1] = T[k] + h;
Y[k+1] = Y[k] + h * [f(x)];
...

Nun habe ich das Problem, dass ich zu Beginn meiner Simulation noch nicht weiß, wie lange diese läuft, da sie ein ereignisbezogenes Abbruchkriterium hat. Insofern kann ich N, welches mein T in gleich große Intervalle partitioniert, nicht wissen. Also ist es wohl "best practice", dass ich ein hohes N wähle und hoffe, dass es ausreicht.
Oder ist es auch möglich, sozusagen das N in jedem Schritt einfach zu ändern, ohne die ergebnisse zu verfälschen?

thx for ideas

Taurin

Die Schrittweite h hat nichts mit der Anzahl deiner Schritte zu tun. Du kannst von Anfang an ein bliebiges h wählen (möglichst klein).
Du kannst h auch während der Simulation ändern. Das hat aber nichts mit der Anzahl deiner Schritte zu tun (das ist sogar üblich - es gibt Lösungsverfahren, die automatisch die Schrittweite an die Funktion anpassen).

Du musst nur daran denken: Mit jedem Schritt wird der Fehler größer!

Das kannst du ausprobieren, indem du eine DGL analytisch löst und dann mit dem Eulerverfahren, z.b. y'(t) = y(t), y(0)=1
Guck dir den Abstand zwischen der echten Kurve und der numerisch berechneten an. Je größer die Zeit, desto größer der Fehler.