Mathc initiation/Fichiers c : c64ea


Sommaire

L'intégrale de flux par le théoreme de la divergence (Partie 1).
Intégrale triple (dydzdx) 1 L'intégrale étudiée (dydzdx)
/* --------------------------------- */
double simpson_dydzdx(

double (*P_f)(double x, double y, double z),


   
double (*Puy)(double x, double z),
double (*Pvy)(double x, double z),
   int    ny,
   
double (*Psz)(double x),
double (*Ptz)(double x),
   int    nz,

double ax,
double bx,
   int nx
)
{
   int i = 0;
double m = 0.;
double M = 0.;

 for(i = 0; i <= nx; i++)
 {
       if(i ==0 || i== nx){m = 1.;}
  else if(fmod(i,2) == 0){m = 2.;}
  else                   {m = 4.;}

  M += m * intz_dydzdx(      (*P_f),
                               


                             (*Puy),
                             (*Pvy),
                                ny,
                             
                             (*Psz),
                             (*Ptz),
                                nz,
  
                             (ax+i*(bx-ax)/nx));
 }

  return( ((bx -ax)*M) / (3*nx) );
}
/* --------------------------------- */
/* --------------------------------- */
double flux_dydzdx(

double (*P_M)(double x, double y, double z),
double (*P_N)(double x, double y, double z),
double (*P_P)(double x, double y, double z),

double (*P_uy)(double x, double z),
double (*P_vy)(double x, double z),
   int ny,

double (*P_sz)(double x),
double (*P_tz)(double x),
   int nz,

double ax,
double bx,
   int nx
)
{
   int i = 0;
double m = 0.;
double M = 0.;

 for(i = 0; i <= nx; i++)
 {
       if(i ==0 || i== nx){m = 1.;}
  else if(fmod(i,2) == 0){m = 2.;}
  else                   {m = 4.;}

  M += m * flux_z_dydzdx(    (*P_M),
                             (*P_N),
                             (*P_P),
                             
                             (*P_uy),
                             (*P_vy),
                                 ny,
                               
                             (*P_sz),
                             (*P_tz),
                                 nz,
                             
                             (ax+i*(bx-ax)/nx));
 }

  return( ((bx -ax)*M) / (3*nx) );
}
/* --------------------------------- */

Comparons les deux fonctions.

Dans les trois premières colonnes, il y a la fonction de référence pour calculer une intégrale triple par la méthode de Simpson. Dans les deuxièmes colonnes il y a les trois fonctions pour calculer l'intégrale de flux par le théorème de la divergence.


On peut remarquer qu'il y a des fonctions supplémentaires en entrée. (Voir (*P_M),(*P_N),(*P_P),). Il y a aussi le paramètre h pour calculer les dérivées partielles. (Voir fxyz_x());

Dans la troisième partie de la fonction étudiée, il y a le calcul des dérivées partielles (divergence) au lieu d'un simple appel à la fonction f.


En comparant ces fonctions aux fonctions de référence, on voit immédiatement l'analogie qu'il existe entre ces fonctions.


L'intégrale de flux par le théoreme de la divergence (Partie 2).
Intégrale triple (dydzdx) 2 L'intégrale étudiée (dydzdx)
/* --------------------------------- */
double intz_dydzdx(

double (*P_f)(double x, double y, double z),
  

 
double (*Puy)(double x, double z),
double (*Pvy)(double x, double z),
   int    ny,

double (*Psz)(double x),
double (*Ptz)(double x),
   int    nz,

double x
)
{
   int i = 0;
double m = 0.;
double M = 0.;

 for(i = 0; i <= nz; i++)
 {
       if(i ==0 || i== nz){m = 1.;}
  else if(fmod(i,2) == 0 ){m = 2.;}
  else                    {m = 4.;}

  M += m * inty_dydzdx(     (*P_f),
                           


                             (*Puy),
                             (*Pvy),
                                ny,
                               
 (((*Psz)(x))+i*(((*Ptz)(x))-((*Psz)(x)))/nz),

                                x);
 }

 return( ((((*Ptz)(x)) -((*Psz)(x)))*M) / (3*nz) );
}
/* --------------------------------- */
/* --------------------------------- */
double flux_z_dydzdx(

double (*P_M)(double x, double y, double z),
double (*P_N)(double x, double y, double z),
double (*P_P)(double x, double y, double z),

double (*P_uy)(double x, double z),
double (*P_vy)(double x, double z),
   int ny,
   
double (*P_sz)(double x),
double (*P_tz)(double x),
   int nz,

double x
)
{
   int i = 0;
double m = 0.;
double M = 0.;

 for(i = 0; i <= nz; i++)
 {
       if(i ==0 || i== nz){m = 1.;}
  else if(fmod(i,2) == 0 ){m = 2.;}
  else                    {m = 4.;}

  M += m * flux_y_dydzdx(    (*P_M),
                             (*P_N),
                             (*P_P),
                             
                             (*P_uy),
                             (*P_vy),
                                 ny,
                               
 (((*P_sz)(x))+i*(((*P_tz)(x))-((*P_sz)(x)))/nz),
                           
                                x);
 }

 return( ((((*P_tz)(x)) -((*P_sz)(x)))*M) / (3*nz) );
}
/* --------------------------------- */


L'intégrale de flux par le théoreme de la divergence (Partie 3).
Intégrale triple (dydzdx) 3 L'intégrale étudiée (dydzdx)
/* --------------------------------- */
double inty_dydzdx(

double (*P_f)(double x, double y, double z),



double (*Puy)(double x, double z),
double (*Pvy)(double x, double z),
   int    ny,

double z,

double x
)
{

pt3d = n;

   int i = 0;
double m = 0.;
double M = 0.;

 for(i = 0; i <= ny; i++)
 {
       if(i ==0 || i== ny){m = 1.;}
  else if(fmod(i,2) == 0 ){m = 2.;}
  else                    {m = 4.;}

  n.x = x;
  n.y = ((*P_uy)(x,z)) + 
        i*(((*P_vy)(x,z))-((*P_uy)(x,z)))/ny;
  n.z = z;

  M += m * (*P_f)( n.x,
                   n.y, 
                   n.z);
 }

  return(((((*Pvy)(x,z))-((*Puy)(x,z)))*M)
                        /(3*ny));
}
/* --------------------------------- */
/* ---------------------------------- */
double flux_y_dydzdx(

double (*P_M)(double x, double y, double z),
double (*P_N)(double x, double y, double z),
double (*P_P)(double x, double y, double z),

double (*P_uy)(double x, double z),
double (*P_vy)(double x, double z),
   int     ny,
   
double z,

double x
)
{

pt3d n;

   int i = 0;
double m = 0.;
double M = 0.;

 for(i = 0; i <= ny; i++)
 {
       if(i ==0 || i== ny){m = 1.;}
  else if(fmod(i,2) == 0 ){m = 2.;}
  else                    {m = 4.;}

 n.x = x;
 n.y = ((*P_uy)(x,z)) + 
           i*(((*P_vy)(x,z))-((*P_uy)(x,z)))/ny;
 n.z = z;

  M += m * (fxyz_x((*P_M),H,n) + 
            fxyz_y((*P_N),H,n) + 
            fxyz_z((*P_P),H,n)) ;
  }

  return(((((*P_vy)(x,z))-((*P_uy)(x,z)))*M)
                         /(3*ny));
}
/* --------------------------------- */