Probleme avec Runge kutta

Signaler
Messages postés
273
Date d'inscription
samedi 5 juillet 2003
Statut
Membre
Dernière intervention
31 mars 2015
-
Messages postés
1
Date d'inscription
samedi 16 décembre 2006
Statut
Membre
Dernière intervention
13 janvier 2008
-
Salut tout le monde,
Je voudrais savoir si quelqu'un peut me donner un coup de main sur un programme concernant l'algorithme Runge Kutta ordre 4.
Voilà, en fait, j'ai fais deux programmes permettant de résoudre l'équation differentielle du pendule amortie avec mon algo rk4 et les données sont sauvegardées dans un fichier.
Cependant, les données qui sont sauvegardées par les 2 progs sont differentes. Je ne vois pas du tout où est mon erreur. Je voulais savoir si quelqu'un peut me dire
où est ce que ça cloche; si c'est dans mon 1er prog? ou dans le 2ème ?
Si ca vient de mon 1er prog comment faire pour arranger ce problème. Je développe avec Dev-Cpp.
Merci.

Voici mon 1er programme:


#include <stdio.h>
#include <math.h>
#define NMAX 65536
#define Y0 2.0
#define pi 3.1415926535
#define OM 2.0/3.0
#define TAU (2*pi)/(128.0*(2.0/3.0))
double q,b,w;
double f(double x, double y, double z)
{
return y;
}
double g(double x, double y, double z)
{
return -q*y-sin(x)+b*cos(z);
}
double h(double x, double y, double z)
{
return OM;
}


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,k1y,k1z,k2x,k2y,k2z,k3x,k3y,k3z,k4x,k4y,k4z;
double xi,yi,zi;
int np;


xi=*x;
yi=*y;
zi=*z;


k1x = TAU * f(xi,yi,zi);
k1y = TAU * g(xi,yi,zi);
k1z = TAU * h(xi,yi,zi);
k2x = TAU * f(xi+k1x/2.0,yi+k1y/2.0,zi+k1z/2.0);
k2y = TAU * g(xi+k1x/2.0,yi+k1y/2.0,zi+k1z/2.0);
k2z = TAU * h(xi+k1x/2.0,yi+k1y/2.0,zi+k1z/2.0);
k3x = TAU * f(xi+k2x/2.0,yi+k2y/2.0,zi+k2z/2.0);
k3y = TAU * g(xi+k2x/2.0,yi+k2y/2.0,zi+k2z/2.0);
k3z = TAU * h(xi+k2x/2.0,yi+k2y/2.0,zi+k2z/2.0);
k4x = TAU * f(xi+k3x,yi+k3y,zi+k3z);
k4y = TAU * g(xi+k3x,yi+k3y,zi+k3z);
k4z = TAU * h(xi+k3x,yi+k3y,zi+k3z);


*x = *x + (k1x+ 2.0*(k2x+k3x)+k4x)/6.0;
*y = *y + (k1y+ 2.0*(k2y+k3y)+k4y)/6.0;
*z = *z + (k1z+ 2.0*(k2z+k3z)+k4z)/6.0;
np = *x/(2.0*pi)+0.5;
*x = *x-2.0*pi*np;
}
int main(int argc, char *argv[])
{
double a1,a2,a3;
int i,n;
extern double q,b,w;
FILE *fichier;
fichier=fopen("data.txt","wt");
void rk4(a1, a2, a3);
a1=0.0;
a2=Y0;
a3=OM;
q=0.5;
b=1.15;
w=2.0/3.0;
for(i=1;i<=n;i++){
rk4(&a1,&a2,&a3);
fprintf(fichier,"%16.8lf %16.8lf\n",a1,a2);
}
fclose(fichier);
system("PAUSE");
return 0;
}


Mon second prog:


#include <stdio.h>
#include <stdlib.h>
#include <math.h>


#define g1(y1,y2,t) (y2)
#define g2(y1,y2,t) (-q*y2-sin(y1)+b*cos(w*t))
#define NMAX 65536
double y[2][NMAX];
double q,b,w;


int main(int argc, char *argv[])
{
int i,n,np,l;
double pi,h,t,y1,y2;
double dk11,dk21,dk12,dk22,dk13,dk23,dk14,dk24;
FILE *fichier;
fichier=fopen("data.txt","wt");
n = NMAX;
pi = 3.1415926535;
w = 2.0/3.0;
q = 0.5;
b = 1.15;
h = 2*pi/(l*w);
l=128.0;
y[0][0] = 0.0;
y[1][0] = 2.0;
for (i = 0; i < n-1; ++i)
{
t = h*(i+1);
y1 = y[0][i];
y2 = y[1][i];
dk11 = h*g1(y1,y2,t);
dk21 = h*g2(y1,y2,t);
dk12 = h*g1((y1+dk11/2.0),(y2+dk21/2.0),(t+h/2.0));
dk22 = h*g2((y1+dk11/2.0),(y2+dk21/2.0),(t+h/2.0));
dk13 = h*g1((y1+dk12/2.0),(y2+dk22/2.0),(t+h/2.0));
dk23 = h*g2((y1+dk12/2.0),(y2+dk22/2.0),(t+h/2.0));
dk14 = h*g1((y1+dk13),(y2+dk23),(t+h));
dk24 = h*g2((y1+dk13),(y2+dk23),(t+h));
y[0][i+1] = y[0][i]+(dk11+2.0*(dk12+dk13)+dk14)/6.0;
y[1][i+1] = y[1][i]+(dk21+2.0*(dk22+dk23)+dk24)/6.0;
np = y[0][i+1]/(2.0*pi)+0.5;
y[0][i+1] = y[0][i+1]-2.0*pi*np;
fprintf(fichier,"%16.8lf %16.8lf\n", y[0][i],y[1][i]);
}
fclose(fichier);
system("PAUSE");
return 0;
}

Jarod_Delaware
A voir également:

11 réponses

Messages postés
1138
Date d'inscription
mardi 10 juin 2003
Statut
Membre
Dernière intervention
25 janvier 2009
3
Tu as essaye de faire du pas-a-pas et voir a partir de quand les deux programmes ne donne plus les memes valeurs ?

Pourquoi faire simple quand on peut faire compliqué ?
Messages postés
273
Date d'inscription
samedi 5 juillet 2003
Statut
Membre
Dernière intervention
31 mars 2015
2
Salut JCDjcd,
Merci pour ta réponse. En fait, les valeurs diffèrent dès la première ligne. Par exemple les 4 premières valeurs imprimées avec le second programme donne:


<COLGROUP>
<COL style=\"WIDTH: 60pt\" span=2 width=80>

----

0.00000000,
2.00000000,
----

0.14752512,
2.00525627,
----

0.29502501,
1.99938879,
----

0.44170028,
1.98297502

par contre avec le premier prog j'ai:


<COLGROUP>
<COL style=\"WIDTH: 60pt\" span=2 width=80>

----

0.14684113 ,
1.98641098,
----

0.29221305,
1.96020017,
----

0.43521337,
1.92214182,
----

0.57500529,
1.87323367

En fait, j'ai l'impression que mes 2 conditions initiales ne figurent pas dans le fichier de données généré par le premier programme. Je ne comprend pas où est mon erreur de programmation. Je suis pourtant persuadé que l'algo dans le premier prog est bon. Si tu as une idée ou quelqu'un n'hésitez pas parce que moi je ne sais plus quoi faire.
Jarod_Delaware
Messages postés
1138
Date d'inscription
mardi 10 juin 2003
Statut
Membre
Dernière intervention
25 janvier 2009
3
Si ca foire des la premiere boucle il suffit de faire soi-meme le calul a la main et voir plus precisement ou dans le code ca bug, en plus le code n'est pas tres long. verifies d'abord que tu pars avec les memes donnees.

Pourquoi faire simple quand on peut faire compliqué ?
Messages postés
273
Date d'inscription
samedi 5 juillet 2003
Statut
Membre
Dernière intervention
31 mars 2015
2
Ok, je vais essayé de faire comme tu m'as dit. Par contre est ce que tu as une idée de là où ça foire dans mes codes?

Jarod_Delaware
Messages postés
1138
Date d'inscription
mardi 10 juin 2003
Statut
Membre
Dernière intervention
25 janvier 2009
3
pourquoi TAU et h sont differents ? dan la seconde version le "l" est mis a 128 APRES son utilisation juste au dessus

Pourquoi faire simple quand on peut faire compliqué ?
Messages postés
273
Date d'inscription
samedi 5 juillet 2003
Statut
Membre
Dernière intervention
31 mars 2015
2
Salut,
TAU et h ont la même valeur. C'est vrai que dans la seconde version j'ai définis l juste après son utilisation. J'ai corrigé cette erreur mais c'est toujours pareil. Je ne sais plus trop quoi faire. Si tu as une idée elle sera la bien venue.

Jarod_Delaware
Messages postés
1138
Date d'inscription
mardi 10 juin 2003
Statut
Membre
Dernière intervention
25 janvier 2009
3
a partir de quelle ligne de code, les deux algo. sont differents ?

Pourquoi faire simple quand on peut faire compliqué ?
Messages postés
273
Date d'inscription
samedi 5 juillet 2003
Statut
Membre
Dernière intervention
31 mars 2015
2
Salut,
En fait ce qui est bizarre c'est que si je met 1024 itérations les valeurs sont les mêmes, mais au dela de 1024 ça part en sucette. Si tu as une idée du pourquoi.

Jarod_Delaware
Messages postés
1138
Date d'inscription
mardi 10 juin 2003
Statut
Membre
Dernière intervention
25 janvier 2009
3
peut etre que la methode elle meme diverge, mais si tu recommences avec les valeurs en 1000 que tu mets en valeurs initiales, ca foire toujours ?

Pourquoi faire simple quand on peut faire compliqué ?
Messages postés
273
Date d'inscription
samedi 5 juillet 2003
Statut
Membre
Dernière intervention
31 mars 2015
2
Salut,


En fait, tu veux que j'essais de mettre comme conditions initiales des valeurs multiples de 1000 sur x et y? ou bien dans le nombre d'itérations ?
Ca change pas grand chose. En fait, je voulais avoir comme nombre d'itérations des multiples de 2^N afin d'effectuer une FFT sur mes points. J'aurais aimer savoir pourquoi les valeurs sont differentes au delà de 1024. J'ai essayé plusieurs trucs mais ça n'a rien donné.

Jarod_Delaware
Messages postés
1
Date d'inscription
samedi 16 décembre 2006
Statut
Membre
Dernière intervention
13 janvier 2008

salut je suis débutant en c++ et je serai trés reconnaissant si quelqu'un me donne un coup de mains pour programmer rk4 pour x=f(t,x,y) et y=g(t,x,y)
merci