Résolution équation differentielle du pendule par la methode de runge kutta ordre 4

Commenter

Contenu du snippet

Ce programme permet de résoudre numeriquement une équation differentielle par la methode Runge Kutta ordre 4. Cette methode de résolution numerique donne une meilleur precision que l'algorithme d'Euler.
Vous pouvez changer dans le programme le pas TAU et les conditions initiales. Les résultats sont enregistrés dans le fichier output.
Vous pouvez vous amusez à changer les paramètres X0, MU, A et OM.
Ici dans le code source du programme XO, MU, A et OM definissent un cas chaostique pour le pendule (essayer de tracer à l'aide d'un tableur l'espace des phases xy).

Source / Exemple : 


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI acos(-1.0)
#define A 3.0
#define OM 1.0
#define TAU 0.01
#define MU 0.3
#define X0 1.0
#define PI 2*acos(-1.0)

/* Fonction Runge kutta ordre 4 */
void rk4(double *x, double *y , double *z)
{
	double f(double x, double y, double z);
	double g(double x, double y, double z);
	double h(double x, double y, double z);
        double k1x,k2x,k1y,k2y,xi,yi,xi2,yi2,xi3,yi3;
	double k3x,k3y,k4x,k4y;


    	
  • / Condition pour empêcher le pendule de faire un tour sur lui même */

		if(*x>PI)(*x)=*x-DPI;
		else if (*x<-PI)(*x)=DPI+*x;

        /* ALGORITHME RUNGE KUTTA 4 */
	xi = *x;
	yi = *y;
	zi = *z;

	k1x = TAU * f(xi,yi,zi);
	k1y = TAU * g(xi,yi,zi);
	k1z = TAU * h(xi,yi,zi);

	xi = *x + 0.5*k1x;
	yi = *y + 0.5*k1y;
	zi = *z + 0.5*k1z;

	k2x = TAU * f(xi,yi,zi);
	k2y = TAU * g(xi,yi,zi);
	k2z = TAU * h(xi,yi,zi);

	xi2 = *x + 0.5*k2x;
	yi2 = *y + 0.5*k2y;
	zi2 = *z + 0.5*k2z;

	k3x = TAU * f(xi2,yi2,zi2);
	k3y = TAU * g(xi2,yi2,zi2);
	k3z = TAU * h(xi2,yi2,zi2);

	xi3 = *x + k3x;
	yi3 = *y + k3y;
	zi3 = *z + k3z;

	k4x = TAU * f(xi3,yi3,zi3);
	k4y = TAU * g(xi3,yi3,zi3);
	k4z = TAU * h(xi3,yi3,zi3);

    /* INTEGRATION */


        
  • x = *x + (1.0/6.0) * (k1x + k4x) + (1.0/3.0) * (k2x + k3x);

    •     
  • y = *y + (1.0/6.0) * (k1y + k4y) + (1.0/3.0) * (k2y + k3y);	

    •     
  • z = *z + (1.0/6.0) * (k1z + k4z) + (1.0/3.0) * (k2z + k3z);

}
double f(double x, double y,double z)
{
	return y;
}

double g(double x, double y,double z)
{
       return A*cos(z) - sin(x) - MU*y;
}
double h(double x, double y, double z)
{
       return OM;
}

/* PROGRAMME PRINCIPAL */
void main(void)
{
	int nbrIT, i;
	double cn1,cn2,cn3;
	FILE *output;
	output=fopen("donnees","w");

	printf("\nEntrer le nombre d'itérations:);
	scanf("%d", &nbrIT);

	/* Conditions initiales */
	cn1 = X0;
	cn2 = 0.0;
	cn3 = OM;
       
	for(i=1;i<=nbrIT;i++)
	{
		rk4(&cn1, &cn2, &cn3); /* On appel la fonction rk4 */
		fprintf(output,"%f %lf %lf %lf \n", i*TAU, cn1,cn2, cn3); /* On sauve périodiquement les résultats */
	}
	fclose(output); /* On ferme le fichier */
} /* fin programme */

Conclusion :


-

A voir également

S'inscrire maintenant

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.

Du même auteur (Jarod1980)