RÉSOLUTION ÉQUATION DIFFERENTIELLE DU PENDULE PAR LA METHODE DE RUNGE KUTTA ORDRE 4

#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 */