Runge-kutta generique

Description

Voila la methode de Runge-Kutta, qui sert a resoudre des equations differentielles avec un erreur d'ordre 4. Le programme montre que la methode est generique, i.e. qu'elle est appliquable a toute les equations differentielles, et non spécifiques a une partiuclierement.

Pour l'exemple on resout une equation differentielle qui regie une circuit RLC (avec R!=0), pour infomation c'est la meme equation qui est resolu dans ma source sur la methode d'Euler, donc s'y reporter (http://www.cppfrance.com/code.aspx?id=24944).

En capture d'ecran, la comparaison entre le resultat de la methode d'Euler, et cella de Runge-Kutta. (le methode d'Euler est moins bonne).

Les operations sur les vecteurs sont dans math.h et math.c

Source / Exemple :


#ifndef _ODE_H_
	#include "ODE.h"
#endif // _ODE_H_

//----------------------------------------------------------
//	RUNGE-KUTTA
//----------------------------------------------------------
void _AssertRungeKutta(BOOL b,char *text)
{
_ErrorIf(!b,text);
} // _AssertRungeKutta()
//----------------------------------------------------------
BOOL _RungeKuttaODE(
                    P_VECTND  yf,
                    double    x0,
                    P_VECTND  y0,
                    int       n,
                    double    h,
                    BOOL      (*derivFunc)(P_VECTND res,double x,P_VECTND y,void *param),
                    void      *param
                   )
{
P_VECTND  k1,k2,k3,k4,tmp;
int       i,dim;
double    x;

AssertGoodVectND(y0);
AssertRungeKutta(NULL != derivFunc,"RungeKuttaODE : <derivFunc> est nul");

dim = yf->dim;

k1  = CreateVectND(dim);
k2  = CreateVectND(dim);
k3  = CreateVectND(dim);
k4  = CreateVectND(dim);
tmp = CreateVectND(dim);

x   = x0;
_CopyVectND(yf,y0);
for(i=0;i<n;i++)
  {
  // k1 = f(x,y)
  derivFunc(k1,x,yf,param);
  // k2 = f(x+h/2,y+h/2*k1)
  _AddVectND(tmp,yf,k1,1.,0.5*h);
  derivFunc(k2,x+0.5*h,tmp,param);
  // k3 = f(x+h/2,y+h/2*k2)
  _AddVectND(tmp,yf,k2,1.,0.5*h);
  derivFunc(k3,x+0.5*h,tmp,param);
  // k4 = f(x+h,y+h*k3)
  _AddVectND(tmp,yf,k3,1.,h);
  derivFunc(k4,x+h,tmp,param);
  // y(n+1) = y(n) + h*(k1/6 + k2/3 + k3/3 + k4/6)
  // on essaye d'eviter ce bout de code (deux fois le meme argument)
  //  | AddVectND(yf,yf,k1,1.,h/6.);
  //  | AddVectND(yf,yf,k2,1.,h/3.);
  //  | AddVectND(yf,yf,k3,1.,h/3.);
  //  | AddVectND(yf,yf,k4,1.,h/6.);
  // car on peut faire autrement :
  _AddVectND(tmp,yf,k1,1.,h/6.);
  _AddVectND(yf,tmp,k2,1.,h/3.);
  _AddVectND(tmp,yf,k3,1.,h/3.);
  _AddVectND(yf,tmp,k4,1.,h/6.);  

  x += h;
  }

Free(k1);
Free(k2);
Free(k3);
Free(k4);
Free(tmp);
return TRUE;
} // _RungeKuttaODE()

Conclusion :


Voila maintenant comment appelle la fonction RungeKuttaODE pour le circuit RLC :
il vous faut ecrire un fonction qui lie (q',q'') a (q,q'), la voila pour un RLC :
//----------------------------------------------------------
BOOL df(P_VECTND res,double x,P_VECTND y,void *param)
{
P_RLC rlc;
rlc = (P_RLC) param;

Assert(res != y);

res->coord[0] = y->coord[1];
res->coord[1] = -(y->coord[0] + rlc->R * rlc->C * y->coord[1])/(rlc->L * rlc->C);
return TRUE;
} // df()
//----------------------------------------------------------

le "void*para" sert a enregistrer les donnees du systeme d'equations differentielles, ici c'est une structure contenant les valeurs de R, L, et C.

Codes Sources

A voir également

Vous n'êtes pas encore membre ?

inscrivez-vous, c'est gratuit et ça prend moins d'une minute !

Les membres obtiennent plus de réponses que les utilisateurs anonymes.

Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.

Le fait d'être membre vous permet d'avoir des options supplémentaires.