Problème sur une boucle en C++
Résolu/Fermé
turbochris
Messages postés
4
Date d'inscription
dimanche 10 mars 2013
Statut
Membre
Dernière intervention
10 mars 2013
-
10 mars 2013 à 10:55
turbochris Messages postés 4 Date d'inscription dimanche 10 mars 2013 Statut Membre Dernière intervention 10 mars 2013 - 10 mars 2013 à 13:49
turbochris Messages postés 4 Date d'inscription dimanche 10 mars 2013 Statut Membre Dernière intervention 10 mars 2013 - 10 mars 2013 à 13:49
6 réponses
KX
Messages postés
16752
Date d'inscription
samedi 31 mai 2008
Statut
Modérateur
Dernière intervention
31 août 2024
3 019
10 mars 2013 à 11:01
10 mars 2013 à 11:01
Sur ce site, il n'y a pas de devin, il faut que tu nous donnes le code qui pose problème parce que tes explications ne sont pas suffisantes pour t'aider...
turbochris
Messages postés
4
Date d'inscription
dimanche 10 mars 2013
Statut
Membre
Dernière intervention
10 mars 2013
10 mars 2013 à 12:15
10 mars 2013 à 12:15
#include <cstdlib>
#include <ctime>
#include <iomanip>
#include <iostream>
#include <math.h>
#include <fstream>
#include <string.h>
#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
#include <fcntl.h>
#define pi 3.1415926535897932384626433832795
#define mu 1.0
#define sig 1.0
#define mu2 0.15
#define ki 1.0
#define alfa2 0.4
#define kip 1.0
#define k0 0.01
#define lam 0.1
#define lamp 0.0
#define dlam 0.05
#define dv 0.025
#define dr 0.01
#define vfin 0.25
//#define rfin 0.25
#define lamfin 4.0
#define kmax 39000
#define ktrans 70000
#define hh 0.05
#define dim 2
#define dip 2
#define val 1
const double c41[4] = {0.0, 0.5*hh, 0.0, 0.5*hh};
const double c42[4] = {0.5*hh, 0.5*hh, hh, 0.0};
const double c43[4] = {hh/6.0, hh/3.0, hh/3.0, hh/6.0};
//const int nbv = int(ceil(vfin/dv1));
const double alfa22=alfa2*alfa2;
int lk;
double v1, v2, r1, r2;
double vx[dip][kmax], vxp[dip][kmax];
double ver[dip][kmax];
// declaration of functions
double sqr(double u);
void initu(double vx[dip][kmax], double vxp[dip][kmax], double ver[dip][kmax]);
double f2(double xt, double xx[dim], double xp[dim]);
double f4(double xt, double xx[dim], double xp[dim]);
void update_rk4_x(long int nt, double &t, double xx[dim], double xp[dim], double vx[dip][kmax], double vxp[dip][kmax], double ver[dip][kmax]);
double Abs(double u);
using namespace std;
int main ()
{
srand(time(NULL));
//declaration of variable
double t, xx[dim], xp[dim];
double er[dim];
double kxx, lam0;
double rfin;
ofstream fichx ("synchro2a.dat");
if ( fichx.is_open() ) // opening the datas file
{
r1=-0.93; r2=r1; rfin=-0.90;
do
{
v1=0.05; v2=v1;
xx[0]=0.0; xp[0]=0.0;
xx[1]=0.01; xp[1]=0.05;
er[0] = xx[1]-xx[0];
er[1] = xp[1]-xp[0];
t=0.;
update_rk4_x(ktrans,t,xx,xp,vx,vxp,ver); // transient
initu(vx,vxp,ver);
update_rk4_x(kmax,t,xx,xp,vx,vxp,ver);
fichx << ver[0][kmax-1] << " " << ver[1][kmax-1] << "\n";
cout << " r1 = " << r1 << " r2 = " << r2 << " "<< ver[0][kmax-1]<<"\n";
r1 += dr;
r2 += dr;
} while (r1<rfin);
cin >> lk;
fichx.close();
} // if no error when opening fichx
return(0);
} // end of main program
double sqr(double u)
{
return(u*u);
}
void initu(double vx[dip][kmax], double vxp[dip][kmax], double ver[dip][kmax])
{
for (int i=0; i<dip; ++i)
{
for (long int k=0; k<kmax; ++k)
{
vx[i][k] = 0.0; vxp[i][k] = 0.0;
ver[i][k] = 0.0;
}
}
}
double f2(double xt, double xx[dim], double xp[dim])
{
double result, arg1, pterm, a, b;
arg1 = 2.*pi*(xx[0]-v1*xt);
a = -mu*sig*mu2*(xp[0]-v1)-mu*ki*xx[0]*alfa22-mu*kip*k0*xx[0]*sqr(xx[0]);
a -= mu*lam*(xx[0]-xx[1])+mu*lamp*(xx[0]-xx[1])*sqr(xx[0]-xx[1]);
b = -(mu*val*(1./(2.*pi))*(sqr(1.-r1*r1))*sin(arg1))/sqr(1.+r1*r1+2.*r1*cos(arg1));
pterm = a+b;
result = pterm;
return (result);
}
double f4(double xt, double xx[dim], double xp[dim])
{
double result, pterm, arg2, a, b;
arg2 = 2.*pi*(xx[1]-v2*xt);
a = -mu2*(xp[1]-v2)-xx[1]*alfa22-k0*xx[1]*sqr(xx[1]);
a -=lam*(xx[1]-xx[0])+lamp*(xx[1]-xx[0])*sqr(xx[1]-xx[0]);
b = -1.*val/(2.*pi)*(sqr(1.-r2*r2))*sin(arg2)/sqr(1.+r2*r2+2.*r2*cos(arg2));
pterm = a+b;
result = pterm;
return (result);
}
void update_rk4_x(long int nt, double & t, double xx[dim], double xp[dim], double vx[dip][kmax], double vxp[dip][kmax], double ver[dip][kmax])
{
double kr4p[dim], kr4x[dim], xnew[dim], xrec[dim], xpnew[dim], xprec[dim], er[dim], som[dim];
som[0] = 0.0;
som[1] = 0.0;
for (long int k=0; k<nt; ++k)
{
for (int i=0; i<dim; ++i)
{
xpnew[i]=xp[i]; xprec[i]=xp[i];
xnew[i]=xx[i]; xrec[i]=xx[i];
}
for (int i4=0; i4<4; ++i4)
{
t += c41[i4];
kr4x[0]=xpnew[0];
kr4x[1]=xpnew[1];
kr4p[0]=f2(t,xnew,xpnew);
kr4p[1]=f4(t,xnew,xpnew);
xnew[0] = xrec[0]+kr4x[0]*c42[i4]; xp[0] += kr4p[0]*c43[i4]; xnew[1] = xrec[1]+kr4x[1]*c42[i4]; xx[0] += kr4x[0]*c43[i4];
xpnew[0] = xprec[0]+kr4p[0]*c42[i4]; xp[1] += kr4p[1]*c43[i4]; xpnew[1] = xprec[1]+kr4p[1]*c42[i4]; xx[1] += kr4x[1]*c43[i4];
} // end of i4 loop
er[0] = Abs(xx[1]-xx[0]);
er[1] = Abs(xp[1]-xp[0]);
if (k<kmax)
for (int i=0; i<dip; ++i)
{
som[i] += er[i];
vx[i][k] = xx[i];
vxp[i][k] = xp[i];
ver[i][k] = som[i];
}
} // end of time loop
} // ******* end of updating x **********
double Abs(double u)
{
if (u>=0)
return(u);
if (u<0)
return(-u);
}
#include <ctime>
#include <iomanip>
#include <iostream>
#include <math.h>
#include <fstream>
#include <string.h>
#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
#include <fcntl.h>
#define pi 3.1415926535897932384626433832795
#define mu 1.0
#define sig 1.0
#define mu2 0.15
#define ki 1.0
#define alfa2 0.4
#define kip 1.0
#define k0 0.01
#define lam 0.1
#define lamp 0.0
#define dlam 0.05
#define dv 0.025
#define dr 0.01
#define vfin 0.25
//#define rfin 0.25
#define lamfin 4.0
#define kmax 39000
#define ktrans 70000
#define hh 0.05
#define dim 2
#define dip 2
#define val 1
const double c41[4] = {0.0, 0.5*hh, 0.0, 0.5*hh};
const double c42[4] = {0.5*hh, 0.5*hh, hh, 0.0};
const double c43[4] = {hh/6.0, hh/3.0, hh/3.0, hh/6.0};
//const int nbv = int(ceil(vfin/dv1));
const double alfa22=alfa2*alfa2;
int lk;
double v1, v2, r1, r2;
double vx[dip][kmax], vxp[dip][kmax];
double ver[dip][kmax];
// declaration of functions
double sqr(double u);
void initu(double vx[dip][kmax], double vxp[dip][kmax], double ver[dip][kmax]);
double f2(double xt, double xx[dim], double xp[dim]);
double f4(double xt, double xx[dim], double xp[dim]);
void update_rk4_x(long int nt, double &t, double xx[dim], double xp[dim], double vx[dip][kmax], double vxp[dip][kmax], double ver[dip][kmax]);
double Abs(double u);
using namespace std;
int main ()
{
srand(time(NULL));
//declaration of variable
double t, xx[dim], xp[dim];
double er[dim];
double kxx, lam0;
double rfin;
ofstream fichx ("synchro2a.dat");
if ( fichx.is_open() ) // opening the datas file
{
r1=-0.93; r2=r1; rfin=-0.90;
do
{
v1=0.05; v2=v1;
xx[0]=0.0; xp[0]=0.0;
xx[1]=0.01; xp[1]=0.05;
er[0] = xx[1]-xx[0];
er[1] = xp[1]-xp[0];
t=0.;
update_rk4_x(ktrans,t,xx,xp,vx,vxp,ver); // transient
initu(vx,vxp,ver);
update_rk4_x(kmax,t,xx,xp,vx,vxp,ver);
fichx << ver[0][kmax-1] << " " << ver[1][kmax-1] << "\n";
cout << " r1 = " << r1 << " r2 = " << r2 << " "<< ver[0][kmax-1]<<"\n";
r1 += dr;
r2 += dr;
} while (r1<rfin);
cin >> lk;
fichx.close();
} // if no error when opening fichx
return(0);
} // end of main program
double sqr(double u)
{
return(u*u);
}
void initu(double vx[dip][kmax], double vxp[dip][kmax], double ver[dip][kmax])
{
for (int i=0; i<dip; ++i)
{
for (long int k=0; k<kmax; ++k)
{
vx[i][k] = 0.0; vxp[i][k] = 0.0;
ver[i][k] = 0.0;
}
}
}
double f2(double xt, double xx[dim], double xp[dim])
{
double result, arg1, pterm, a, b;
arg1 = 2.*pi*(xx[0]-v1*xt);
a = -mu*sig*mu2*(xp[0]-v1)-mu*ki*xx[0]*alfa22-mu*kip*k0*xx[0]*sqr(xx[0]);
a -= mu*lam*(xx[0]-xx[1])+mu*lamp*(xx[0]-xx[1])*sqr(xx[0]-xx[1]);
b = -(mu*val*(1./(2.*pi))*(sqr(1.-r1*r1))*sin(arg1))/sqr(1.+r1*r1+2.*r1*cos(arg1));
pterm = a+b;
result = pterm;
return (result);
}
double f4(double xt, double xx[dim], double xp[dim])
{
double result, pterm, arg2, a, b;
arg2 = 2.*pi*(xx[1]-v2*xt);
a = -mu2*(xp[1]-v2)-xx[1]*alfa22-k0*xx[1]*sqr(xx[1]);
a -=lam*(xx[1]-xx[0])+lamp*(xx[1]-xx[0])*sqr(xx[1]-xx[0]);
b = -1.*val/(2.*pi)*(sqr(1.-r2*r2))*sin(arg2)/sqr(1.+r2*r2+2.*r2*cos(arg2));
pterm = a+b;
result = pterm;
return (result);
}
void update_rk4_x(long int nt, double & t, double xx[dim], double xp[dim], double vx[dip][kmax], double vxp[dip][kmax], double ver[dip][kmax])
{
double kr4p[dim], kr4x[dim], xnew[dim], xrec[dim], xpnew[dim], xprec[dim], er[dim], som[dim];
som[0] = 0.0;
som[1] = 0.0;
for (long int k=0; k<nt; ++k)
{
for (int i=0; i<dim; ++i)
{
xpnew[i]=xp[i]; xprec[i]=xp[i];
xnew[i]=xx[i]; xrec[i]=xx[i];
}
for (int i4=0; i4<4; ++i4)
{
t += c41[i4];
kr4x[0]=xpnew[0];
kr4x[1]=xpnew[1];
kr4p[0]=f2(t,xnew,xpnew);
kr4p[1]=f4(t,xnew,xpnew);
xnew[0] = xrec[0]+kr4x[0]*c42[i4]; xp[0] += kr4p[0]*c43[i4]; xnew[1] = xrec[1]+kr4x[1]*c42[i4]; xx[0] += kr4x[0]*c43[i4];
xpnew[0] = xprec[0]+kr4p[0]*c42[i4]; xp[1] += kr4p[1]*c43[i4]; xpnew[1] = xprec[1]+kr4p[1]*c42[i4]; xx[1] += kr4x[1]*c43[i4];
} // end of i4 loop
er[0] = Abs(xx[1]-xx[0]);
er[1] = Abs(xp[1]-xp[0]);
if (k<kmax)
for (int i=0; i<dip; ++i)
{
som[i] += er[i];
vx[i][k] = xx[i];
vxp[i][k] = xp[i];
ver[i][k] = som[i];
}
} // end of time loop
} // ******* end of updating x **********
double Abs(double u)
{
if (u>=0)
return(u);
if (u<0)
return(-u);
}
KX
Messages postés
16752
Date d'inscription
samedi 31 mai 2008
Statut
Modérateur
Dernière intervention
31 août 2024
3 019
10 mars 2013 à 12:17
10 mars 2013 à 12:17
Donc ta boucle c'est ça :
Mais quel est le problème ? Ça a l'air de bien fonctionner...
r1=-0.93; r2=r1; rfin=-0.90; do { ... r1 += dr; r2 += dr; } while (r1<rfin);
Mais quel est le problème ? Ça a l'air de bien fonctionner...
turbochris
Messages postés
4
Date d'inscription
dimanche 10 mars 2013
Statut
Membre
Dernière intervention
10 mars 2013
10 mars 2013 à 12:47
10 mars 2013 à 12:47
quand tu initialise r à -0.94; les résultats affichés par ver[0][kmax-1] changent
Vous n’avez pas trouvé la réponse que vous recherchez ?
Posez votre question
KX
Messages postés
16752
Date d'inscription
samedi 31 mai 2008
Statut
Modérateur
Dernière intervention
31 août 2024
3 019
10 mars 2013 à 13:31
10 mars 2013 à 13:31
Il faut faire attention aux approximations sur les double !
Quand tu initialises r1=-0.93, la valeur "réelle" de r1 est -0.93000000000000005
Mais quand tu initialises r1=-0.94, la valeur de r1 est -0.93999999999999995 à quoi tu rajoutes 0.01, ce qui donne -0.92999999999999994.
La valeur de r1 n'est donc jamais vraiment égale à -0.93 et son approximation peut-être différente selon la manière dont elle est calculée (-0.94 + 0.01 ≠ -0.93)
Remarque : il en est de même pour r2.
Donc comme tes calculs utilisent intensivement r1 et r2 (dans update_rk4_x qui appelle souvent f2 et f4), les erreurs d'approximations s'accumulent et cela entraîne une bifurcation du résultat final.
Il faudrait donc revoir tout tes calculs pour diminuer l'impact de ces erreurs, par exemple en manipulant des types plus précis (comme long double).
Au passage, il serait bon de revoir ta structure de code en général, car l'utilisation de variables globales comme ça, c'est vraiment très moche !
Quand tu initialises r1=-0.93, la valeur "réelle" de r1 est -0.93000000000000005
Mais quand tu initialises r1=-0.94, la valeur de r1 est -0.93999999999999995 à quoi tu rajoutes 0.01, ce qui donne -0.92999999999999994.
La valeur de r1 n'est donc jamais vraiment égale à -0.93 et son approximation peut-être différente selon la manière dont elle est calculée (-0.94 + 0.01 ≠ -0.93)
Remarque : il en est de même pour r2.
Donc comme tes calculs utilisent intensivement r1 et r2 (dans update_rk4_x qui appelle souvent f2 et f4), les erreurs d'approximations s'accumulent et cela entraîne une bifurcation du résultat final.
Il faudrait donc revoir tout tes calculs pour diminuer l'impact de ces erreurs, par exemple en manipulant des types plus précis (comme long double).
Au passage, il serait bon de revoir ta structure de code en général, car l'utilisation de variables globales comme ça, c'est vraiment très moche !
turbochris
Messages postés
4
Date d'inscription
dimanche 10 mars 2013
Statut
Membre
Dernière intervention
10 mars 2013
10 mars 2013 à 13:49
10 mars 2013 à 13:49
merci pr ton aide