Runge Kutta bei unabhängiger Variable?

[raeuberhotzenplotz]

Hi!

Ich versuche eine, besser zwei - Differentialgleichung iterativ zu integrieren, und zwar nach der Zeit. Die DGL hängt aber nicht mehr direkt von der Zeit ab. (Lotka-Volterra-DGLn).

Das Euler-Verfahren ist da ja recht einfach, aber ich verstehe nicht, wie ich das Problem auf das Runge-Kutta-Verfahren (zunächst RK2) ummünzen kann.

Für das Euler verfahren benutze ich jetzt diese Schleife:

double w_punkt (const double w, const double k)
{ return w*(1-k); }

double k_punkt (const double w, const double k)
{ return (-k)*(1-w);}

int main () {

  double k=0.5;
  double w=0.5;
  double a=0.5;
  double h=0.001;
  double t=0; 

  while (t<100) {

    a +=h*w_punkt(w,k);
    k +=h*k_punkt(w,k);
    t+=h;
    w=a; }

Wenn ich das ausführe und die Werte in eine Datei ausgebe lasse und das ganze plotte schaut das ja schon ganz gut aus.

Damit das noch genauer wird, muss jetzt das Runge-Kutta verfahren her, und ich verstehe nicht ganze, wie ich das anstellen muss, wenn meine Integrationsvariable gar nicht mehr explizit in einer dieser Gleichungen drin steckt? Ich muss ja jetzt einen halben Integrationsschritt weiter, da die "Steigung" bestimmen und dann damit den ganzen Integrationsschritt gehen.
Integriere ich eine "normale" Dgl, zb. eine Parabelgleichung oder sowas, dann ist das ja recht einfach, aber die ist ja von der Integrationsvariable auch abhängig. Diese DGLn hier ja nicht.

Hat jemand vllt. irgendwo ein Beispiel für so eine Fall in Petto, oder einen Link, in der so ein ähnliches Problem behandelt wird? Oder kann mir jemand von der Grundidee her sagen, wie ich hier zur Lösung kommen könnte?

Vielen Dank für die Hilfe!

Grüße,

der raeuberhotzenplotz

[headmyshoulder]

Integrieren ist beim Lösen von DGLs nicht wortwörtlich gemeint. Man integriert zwar formal über dx/dt aber praktisch verkettet man bei Runge-Kutta-Verfahren geschickt mehrere Euler-Schritte.

Schau beispielsweise hier: http://de.wikipedia.org/wiki/Klassisches_Runge-Kutta-Verfahren.

Ich würd an deiner Stelle auch die die DGL in eine Funktion packen, vielleicht so

void lotka_volterra( double *x , double *dxdt )
{
  dxdt[0] = ...
  dxdt[1] = ...
}

Und die Runge-Kutta Solver in eine eigene Funktion packen. Dann kannst Du auch Solver und DGL austauschen. Das wird z.B. in den Numerical Recipes beschrieben. Oder Du nimmst was vorgefertigtes: odeint (um mal ein bisschen Werbung zu machen).

Ich hoffe dass hilft Dir weiter.

[raeuberhotzenplotz]

Hm, nein so ganz verstehe ich das nicht.

Ich verstehe glaube ich das ganze Prinzip nicht.

Vllt. kannst du mir ja mal kurz weiterhelfen, ich habe mir jetzt ein extrem einfaches Problem zurechtgestrickt.

die DGL x'(t)= t/x(t)

Nach Euler habe ich das ganze jetzt so gelöst :

#include <iostream>
#include <cmath>
using namespace std; 

double f(const double t, const double x) {
 return t/x;

}

double h=0.001;

double x=0.5; 
double y=0.5;
double t=0;

int main(){

while (t<5) {

  x+=h*f(t,x); 
  t+=h; 
}
cout << x << endl; 
}

Nach Runge Kutta 2. Ordnung sieht das ganze jetzt so aus:

#include <iostream>
#include <cmath>
using namespace std; 

double f(const double t, const double x) {
 return t/x;

}

double h=0.001;

double x=0.5; 
double y=0.5;
double t=0;

int main(){

while (t<5) {

  x+=h*f(t+(h/2),x+h/2*f(t,x)); 
  t+=h; 
}
cout << x << endl; 
}

Hier klappt alles wunderbar. Ich kann das Verfahren direkt aus dem Lehrbuch übernehmen.

Wenn mir aber die explizite Zeitabhängigkeit fehlt, wie hier:

x'[t]= 1/ x[t]

Funktioniert das Euler verfahren genauso:

#include <iostream>
#include <cmath>
using namespace std; 

double f(const double x) {
 return 1/x;

}

double h=0.001;

double x=0.5; 
double y=0.5;
double t=0;

int main(){

while (t<5) {

  x+=h*f(x); 
  t+=h; 
}
cout << x << endl; 
}

Wenn ich jetzt richtig liege, müsste das ganze doch jetzt so aussehen:

#include <iostream>
#include <cmath>
using namespace std; 

double f(const double x) {
 return 1/x;

}

double h=0.001;

double x=0.5; 
double y=0.5;
double t=0;

int main(){

while (t<5) {

  x+=h*f(x+h/2*f(x)); 
  t+=h; 
}
cout << x << endl; 
}

Vom Wert her könnte es stimmen.

Wenn das stimmt, kannst du mir noch sagen, ob mein Euler-Ansatz im ersten Beitrag richtig ist? Ich hab das gefühl, vor lauter Wald die Bäume nicht mehr zu sehen ..

Danke für die Hilfe!

[headmyshoulder]

raeuberhotzenplotz schrieb:

Vom Wert her könnte es stimmen.

Wenn das stimmt, kannst du mir noch sagen, ob mein Euler-Ansatz im ersten Beitrag richtig ist? Ich hab das gefühl, vor lauter Wald die Bäume nicht mehr zu sehen ..

Danke für die Hilfe!

Was Du über das 1D Beispiel schreibst ist fast richtig. Ich glaub Du musst nur dt/2 im zweiten Schritt bei dem RK 2.ter Ordnung nehmen. Also

while (t<5) {  
  x+=h/2*f(t+(h/2),x+h/2*f(t,x)); 
  t+=h; 
}

Auch dein Euler für das Lotka-Volterra-System ist richtig. Ich glaub Du hast das schon fast verstanden. Du musst halt bei dem RK2 und dem Lotka-Volterra einen Zwischenschritt explizit hinschreiben. Das geht in einer einzigen Zeile nicht, da w und k voneinander abhängen. Berechnen einfach erst f(x) und mach dann die einzelnen Schritte machen:

while( t < tmax )
{
  double dk = k_punkt( w , k );
  double dw = w_punkt( w , k );
  k = k + dt / 2 * dk;
  w = w + dt / 2 * dw;
  dk = k_punkt( w , k );
  dw = w_punkt( w , k );
  k = k + dt / 2 * dk;
  w = w + dt / 2 * dw;
}

Das würd für höhere Ordnung übrigens nur noch komplizierter. Ich würde Dir auf jedenfall empfehlen, hier zwischen der DGL und dem Runge-Kutta-Solver trennen. Du kannst das obige Beispiel auch so schreiben:

void einfache_dgl( double* x , double *dxdt , double t )
{
  dxdt[0] = t / x[0];
}

void euler( double *x , double t , double dt , int n , void (*dgl)( double* , double* , double ) )
{
  double dxdt[n];
  dgl( x , dxdt , t );
  for( int i=0 ; i<n ; ++i ) x[i] = x[i] + dt * dxdt[i];
}

void rk2( double *x , double t , double dt , int n , void (*dgl)( double* , double* , double ) )
{
  double dxdt[n];
  dgl( x , dxdt , t );
  for( int i=0 ; i<n ; ++i ) x[i] = x[i] + dt * dxdt[i];
  dgl( x , dxdt , t + dt/2 );
  for( int i=0 ; i<n ; ++i ) x[i] = x[i] + dt * dxdt[i];
}

und benutzt wird das dann so

double x[1] = 0.5;
euler( x , 0.5 , 0.5 , 1 , einfache_dgl );

oder für das Lotka Volterra System

double x[2] = ...;
double t = ...;
double dt = ...;
rk2( x , t , dt , 2 , lotka_volterra );

[raeuberhotzenplotz]

headmyshoulder schrieb:

Was Du über das 1D Beispiel schreibst ist fast richtig. Ich glaub Du musst nur dt/2 im zweiten Schritt bei dem RK 2.ter Ordnung nehmen.

Sicher? Soweit ich das verstanden habe nehme ich bei RK2 den "Steigungswert" des halben Iterationsschritt, und iteriere damit dann aber einen ganzen Schritt, also h.
Ich kann mich natürlich auch täuschen, aber bei h/2 kommt dann ein anderer Wert heraus, und ich hab den Wert mit Wolfram Alpha auch rausbekommen(also mit ganzem h).

Deinen Code :

void einfache_dgl( double* x , double *dxdt , double t )
{
  dxdt[0] = t / x[0];
}

void euler( double *x , double t , double dt , int n , void (*dgl)( double* , double* , double ) )
{
  double dxdt[n];
  dgl( x , dxdt , t );
  for( int i=0 ; i<n ; ++i ) x[i] = x[i] + dt * dxdt[i];
}

void rk2( double *x , double t , double dt , int n , void (*dgl)( double* , double* , double ) )
{
  double dxdt[n];
  dgl( x , dxdt , t );
  for( int i=0 ; i<n ; ++i ) x[i] = x[i] + dt * dxdt[i];
  dgl( x , dxdt , t + dt/2 );
  for( int i=0 ; i<n ; ++i ) x[i] = x[i] + dt * dxdt[i];
}

Schaffe ich nicht zu 100% nachzuvollziehen. Sind das dxdt[n], x[1] etc. nicht arrays? Speichere ich dann damit alle iterationsschritte in dem Array ab und lass mir am Schluss nur den letzten ausgeben?

habe das ganze jetzt so weit umgeschrieben bei mir:

#include <iostream>
#include <cmath>

using namespace std; 

double dw(const double w, const double k)
{ return w*(1-k); }

double dk (const double w, const double k)
{ return (-k)*(1-w);}

void euler ( double tmax, double w, double k, double t, double a, const double h) {
  while (t<tmax)  {

    a +=h*dw(w,k);
    k +=h*dk(w,k);
    w=a;
    t+=h;
  }
cout << "Iteration nach Euler:" << endl; 
cout << "W(t=100) =" << w << endl; 
cout << "K(t=100) =" << k << endl; 

}

void rk2 (double tmax, double w, double k, double t, double a, const double h) {

  while (t<tmax) {

   a +=h/2*dw(w,k); 
   k +=h/2*dk(w,k); 
   w=a; 
   a += h*dw(w,k);
   k += h*dk(w,k);
   w=a;
   t+=h;

  }
cout << "Iteration nach RK2:" << endl; 
cout << "W(t=100) =" << w << endl; 
cout << "K(t=100) =" << k << endl; 
   }

void rk4 (double tmax, double w, double k, double t, double a, const double h) {

  while (t<tmax) {
    double w1= w + h/2 *dw(w,k);
    double k1= k + h/2 *dk(w,k); 
    double Kw1= dw(w1,k1);
    double Kk1= dk(w1,k1);
    double w2= w + h/2*dw(w1,k1);
    double k2= k + h/2*dk(w1,k1);

int main () {

euler(100, 0.5, 0.5, 0, 0.5, 0.001);
rk2 (100, 0.5, 0.5, 0, 0.5, 0.001);
}

Es funktioniert auch soweit (Ohne den RK4 - Teil, damit kämpfe ich gerade). Leider weichen die Werte von euler und rk2 soweit voneinander ab, dass ich nicht mehr weiß, ob ich überhaupt noch im Bereich des richtigen bin.

Ich arbeite schon das ganze Wochenende nur an dem Mist und bekomme gar nichts hin ... sehr deprimierend.

Danke auf jeden Fall für deine Hilfe!

Grüße!

[headmyshoulder]

raeuberhotzenplotz schrieb:

Sicher? Soweit ich das verstanden habe nehme ich bei RK2 den "Steigungswert" des halben Iterationsschritt, und iteriere damit dann aber einen ganzen Schritt, also h.
Ich kann mich natürlich auch täuschen, aber bei h/2 kommt dann ein anderer Wert heraus, und ich hab den Wert mit Wolfram Alpha auch rausbekommen(also mit ganzem h).

Sry, Du hast recht. Da muss nur dt hin.

raeuberhotzenplotz schrieb:

Schaffe ich nicht zu 100% nachzuvollziehen. Sind das dxdt[n], x[1] etc. nicht arrays? Speichere ich dann damit alle iterationsschritte in dem Array ab und lass mir am Schluss nur den letzten ausgeben?

Ja, dxdt[n] sind arrays. Aber n ist die Dimension der DGL, also für das einfach Beispiel 1 und für Lotka-Volterra 2. In der Fkt. wird nur ein Schritt ausgeführt. Und man kann Euler und RK2 auch für andere DGLs benutzen. Darum gehts. Das solltest Du auch für RK4 machen. Stell Dir vor, Du willst morgen das Lorenz-System oder was anderes mit RK4 lösen. Dann kannst Du einfach die fertige RK4-Routine nehmen und musst nur eine neue Funktion schreiben, welche dx1/dt, dx2/dt, dx3/dt berechnet.

raeuberhotzenplotz schrieb:

Es funktioniert auch soweit (Ohne den RK4 - Teil, damit kämpfe ich gerade). Leider weichen die Werte von euler und rk2 soweit voneinander ab, dass ich nicht mehr weiß, ob ich überhaupt noch im Bereich des richtigen bin.

Jepp, euler und rk2 müssen was anderes liefern. Die Ergebnisse solten aber nicht zuweit auseinander sein. Wenn Du damit sehr viele Schritte machst, sollten die Kurven bei Lotka-Volterra fast übereinander liegen. Dass kannst Du ganz einfach plotten.