Probleme avec Runge kutta

Jarod1980 Messages postés 273 Date d'inscription samedi 5 juillet 2003 Statut Membre Dernière intervention 31 mars 2015 - 1 mars 2005 à 16:36
cs_spirit31 Messages postés 1 Date d'inscription samedi 16 décembre 2006 Statut Membre Dernière intervention 13 janvier 2008 - 13 janv. 2008 à 20:06
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

11 réponses

cs_JCDjcd Messages postés 1138 Date d'inscription mardi 10 juin 2003 Statut Membre Dernière intervention 25 janvier 2009 4
1 mars 2005 à 18:05
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é ?
0
Jarod1980 Messages postés 273 Date d'inscription samedi 5 juillet 2003 Statut Membre Dernière intervention 31 mars 2015 2
1 mars 2005 à 18:28
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
0
cs_JCDjcd Messages postés 1138 Date d'inscription mardi 10 juin 2003 Statut Membre Dernière intervention 25 janvier 2009 4
1 mars 2005 à 18:36
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é ?
0
Jarod1980 Messages postés 273 Date d'inscription samedi 5 juillet 2003 Statut Membre Dernière intervention 31 mars 2015 2
1 mars 2005 à 18:50
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
0

Vous n’avez pas trouvé la réponse que vous recherchez ?

Posez votre question
cs_JCDjcd Messages postés 1138 Date d'inscription mardi 10 juin 2003 Statut Membre Dernière intervention 25 janvier 2009 4
1 mars 2005 à 22:29
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é ?
0
Jarod1980 Messages postés 273 Date d'inscription samedi 5 juillet 2003 Statut Membre Dernière intervention 31 mars 2015 2
2 mars 2005 à 10:23
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
0
cs_JCDjcd Messages postés 1138 Date d'inscription mardi 10 juin 2003 Statut Membre Dernière intervention 25 janvier 2009 4
2 mars 2005 à 13:55
a partir de quelle ligne de code, les deux algo. sont differents ?

Pourquoi faire simple quand on peut faire compliqué ?
0
Jarod1980 Messages postés 273 Date d'inscription samedi 5 juillet 2003 Statut Membre Dernière intervention 31 mars 2015 2
10 mars 2005 à 10:45
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
0
cs_JCDjcd Messages postés 1138 Date d'inscription mardi 10 juin 2003 Statut Membre Dernière intervention 25 janvier 2009 4
12 mars 2005 à 18:46
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é ?
0
Jarod1980 Messages postés 273 Date d'inscription samedi 5 juillet 2003 Statut Membre Dernière intervention 31 mars 2015 2
13 mars 2005 à 13:02
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
0
cs_spirit31 Messages postés 1 Date d'inscription samedi 16 décembre 2006 Statut Membre Dernière intervention 13 janvier 2008
13 janv. 2008 à 20:06
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
0
Rejoignez-nous