77891.cpp

   1 // vim: set tabstop=4 softtabstop=4 shiftwidth=4 expandtab:
   2 //
   3 // Trabajo Práctico II de Análisis Numérico I
   4 // Este programa resuelve un sistema de ecuaciones diferenciales
   5 // resultante de un problema físico de oscilación de líquidos.
   6 // Copyright (C) 2002 Leandro Lucarella <leandro@lucarella.com.ar>
   7 //
   8 // Este programa es Software Libre; usted puede redistribuirlo
   9 // y/o modificarlo bajo los términos de la "GNU General Public
  10 // License" como lo publica la "FSF Free Software Foundation",
  11 // o (a su elección) de cualquier versión posterior.
  12 //
  13 // Este programa es distribuido con la esperanza de que le será
  14 // útil, pero SIN NINGUNA GARANTIA; incluso sin la garantía
  15 // implícita por el MERCADEO o EJERCICIO DE ALGUN PROPOSITO en
  16 // particular. Vea la "GNU General Public License" para más
  17 // detalles.
  18 //
  19 // Usted debe haber recibido una copia de la "GNU General Public
  20 // License" junto con este programa, si no, escriba a la "FSF
  21 // Free Software Foundation, Inc.", 59 Temple Place - Suite 330,
  22 // Boston, MA  02111-1307, USA.
  23 //
  24 // $URL$
  25 // $Date$
  26 // $Rev$
  27 // $Author$
  28 //
  29
  30 #include <iostream.h>
  31 #include <stdlib.h>
  32 #include <math.h>
  33
  34 // Tipos de datos.
  35 typedef unsigned int Indice;
  36 typedef float Numero;
  37 struct Datos;
  38 typedef Numero (*Friccion)( Datos&, Numero, Numero );
  39 typedef Numero (*Funcion)( Datos&, Numero, Numero, Numero );
  40 typedef void (*Metodo)( Datos& );
  41 struct Datos {
  42     Numero   to; // Tiempo inicial
  43     Numero   tf; // Tiempo final
  44     Numero   k;  // Paso
  45     Numero   xo; // X inicial = X(to) = Z(to)
  46     Numero   yo; // Y inicial = Y(to) = X'(to) = Z'(to)
  47     Numero   D;  // Diametro del tubo
  48     Numero   n;  // Viscosidad cinemática del líquido
  49     Numero   g;  // Aceleración de la gravedad
  50     Numero   f;  // Factor de fricción
  51     Numero   L;  // Longitud del tubo
  52     Numero   Le; // Factor de pérdida de carga
  53     Numero   G;  // Factor geométrico
  54     Funcion  fx; // Función asociada a la derivada de x
  55     Funcion  fy; // Función asociada a la derivada de y
  56     Friccion fi; // Factor asociado a pérdidas por fricción
  57 };
  58
  59 // Factor de fricción nulo.
  60 Numero friccionNula( Datos&, Numero, Numero );
  61
  62 // Factor de fricción laminar.
  63 Numero friccionLaminar( Datos&, Numero, Numero );
  64
  65 // Factor de fricción turbulenta.
  66 Numero friccionTurbulenta( Datos&, Numero, Numero );
  67
  68 // Factor de fricción turbulenta (con depósitos).
  69 Numero friccionTurbulentaConDepositos( Datos&, Numero, Numero );
  70
  71 // Función asociada a la primera derivada (x'=z').
  72 Numero funcionX( Datos&, Numero, Numero, Numero );
  73
  74 // Función asociada a la segunda derivada (y'=x''=z'').
  75 Numero funcionY( Datos&, Numero, Numero, Numero );
  76
  77 // Calcula la ecuación diferencial por el método de Euler.
  78 void euler( Datos& );
  79
  80 // Calcula la ecuación diferencial por el método de Runge-Kutta de órden 4.
  81 void rk4( Datos& );
  82
  83 // Calcula la ecuación diferencial por el método de Nystrom.
  84 void nystrom( Datos& );
  85
  86 // Constantes.
  87 const Numero   DEFAULT_N  = 1000,        // Cantidad de pasos por defecto
  88                DEFAULT_to = 0.0,
  89                DEFAULT_tf = 440.0,
  90                DEFAULT_xo = 2.9866369,
  91                DEFAULT_yo = 0.0,
  92                DEFAULT_D  = 0.5,
  93                DEFAULT_n  = 1e-6,
  94                DEFAULT_g  = 9.8,
  95                DEFAULT_f  = 0.03,
  96                DEFAULT_L  = 892.0,
  97                DEFAULT_Le = 1070.4,
  98                DEFAULT_G  = 2.0;
  99 const Funcion  DEFAULT_fx = &funcionX,
 100                DEFAULT_fy = &funcionY;
 101 const Friccion DEFAULT_fi = &friccionNula;
 102
 103 int main( int argc, char* argv[] ) {
 104
 105     // Se fija que tenga los argumentos necesarios para correr.
 106     if ( argc < 2 ) {
 107         cerr << "Faltan argumentos. Modo de uso:" << endl;
 108         cerr << "\t" << argv[0] << " método fi G tf N to xo yo D n g f L Le" << endl;
 109         cerr << "Desde fi en adelante son opcionales. Los valores por defecto son:" << endl;
 110         cerr << "\tfi  = n" << endl;
 111         cerr << "\tG   = " << DEFAULT_G  << endl;
 112         cerr << "\ttf  = " << DEFAULT_tf << endl;
 113         cerr << "\tN   = " << DEFAULT_N  << endl;
 114         cerr << "\tto  = " << DEFAULT_to << endl;
 115         cerr << "\txo  = " << DEFAULT_xo << endl;
 116         cerr << "\tyo  = " << DEFAULT_yo << endl;
 117         cerr << "\tD   = " << DEFAULT_D  << endl;
 118         cerr << "\tn   = " << DEFAULT_n  << endl;
 119         cerr << "\tg   = " << DEFAULT_g  << endl;
 120         cerr << "\tf   = " << DEFAULT_f  << endl;
 121         cerr << "\tL   = " << DEFAULT_L  << endl;
 122         cerr << "\tLe  = " << DEFAULT_Le << endl;
 123         return EXIT_FAILURE;
 124     }
 125
 126     // Selecciona el método deseado.
 127     Metodo metodo = NULL;
 128     switch ( char( argv[1][0] ) ) {
 129         case 'e':
 130             metodo = &euler;
 131             break;
 132         case 'r':
 133             metodo = &rk4;
 134             break;
 135         case 'n':
 136             metodo = &nystrom;
 137             break;
 138         default:
 139             cerr << "Debe especificar un método válido:" << endl;
 140             cerr << "\te: Euler" << endl;
 141             cerr << "\tr: Runge-Kutta 4" << endl;
 142             cerr << "\tn: Nystrom" << endl;
 143             return EXIT_FAILURE;
 144     }
 145
 146     // Se inicializan los datos.
 147     Numero N = DEFAULT_N;
 148     Datos D;
 149     D.to = DEFAULT_to;
 150     D.tf = DEFAULT_tf;
 151     D.k  = ( D.tf - D.to ) / N;
 152     D.xo = DEFAULT_xo;
 153     D.yo = DEFAULT_yo;
 154     D.D  = DEFAULT_D;
 155     D.n  = DEFAULT_n;
 156     D.g  = DEFAULT_g;
 157     D.f  = DEFAULT_f;
 158     D.L  = DEFAULT_L;
 159     D.Le = DEFAULT_Le;
 160     D.G  = DEFAULT_G;
 161     D.fx = DEFAULT_fx;
 162     D.fy = DEFAULT_fy;
 163     D.fi = DEFAULT_fi;
 164
 165     // Si se pasaron datos como argumento, se los va agregando.
 166     switch ( argc ) {
 167         case 15: D.Le = Numero( atof( argv[14] ) );
 168         case 14: D.L  = Numero( atof( argv[13] ) );
 169         case 13: D.f  = Numero( atof( argv[12] ) );
 170         case 12: D.g  = Numero( atof( argv[11] ) );
 171         case 11: D.n  = Numero( atof( argv[10] ) );
 172         case 10: D.D  = Numero( atof( argv[9] ) );
 173         case 9:  D.yo = Numero( atof( argv[8] ) );
 174         case 8:  D.xo = Numero( atof( argv[7] ) );
 175         case 7:  D.to = Numero( atof( argv[6] ) );
 176         case 6:  N    = Numero( atof( argv[5] ) );
 177         case 5:  D.tf = Numero( atof( argv[4] ) );
 178                  // Se recalcula el paso (k) si se cambio to, N o tf.
 179                  D.k  = ( D.tf - D.to ) / N;
 180         case 4:  D.G  = Numero( atof( argv[3] ) );
 181         case 3:  switch ( char( argv[2][0] ) ) { // Tipo de fricción
 182                     case 'n':
 183                         D.fi = &friccionNula;
 184                         break;
 185                     case 'l':
 186                         D.fi = &friccionLaminar;
 187                         break;
 188                     case 't':
 189                         D.fi = &friccionTurbulenta;
 190                         break;
 191                     case 'd':
 192                         D.fi = &friccionTurbulentaConDepositos;
 193                         break;
 194                     default:
 195                         cerr << "Debe especificar un tipo de fricción válido:" << endl;
 196                         cerr << "\tn: Nula" << endl;
 197                         cerr << "\tl: Laminar" << endl;
 198                         cerr << "\tt: Turbulenta" << endl;
 199                         cerr << "\td: Turbulenta (con depósitos)" << endl;
 200                         return EXIT_FAILURE;
 201                  }
 202     }
 203
 204     // Imprime el paso utilizado.
 205     cerr << "Paso k = " << D.k << endl;
 206
 207     // Ejecuta el método correspondiente con los datos correspondientes.
 208     metodo( D );
 209
 210     return EXIT_SUCCESS;
 211
 212 }
 213
 214 Numero friccionNula( Datos& D, Numero x, Numero y ) {
 215
 216     return 0.0;
 217
 218 }
 219
 220 Numero friccionLaminar( Datos& D, Numero x, Numero y ) {
 221
 222     return 32 * D.n / ( D.D * D.D );
 223
 224 }
 225
 226 Numero friccionTurbulenta( Datos& D, Numero x, Numero y ) {
 227
 228     return D.f * fabs( y ) / ( 2 * D.D );
 229
 230 }
 231
 232 Numero friccionTurbulentaConDepositos( Datos& D, Numero x, Numero y ) {
 233
 234     return D.f * fabs( y ) * D.Le / ( 2 * D.D * D.L);
 235
 236 }
 237
 238 Numero funcionX( Datos& D, Numero t, Numero x, Numero y ) {
 239
 240     return y;
 241
 242 }
 243
 244 Numero funcionY( Datos& D, Numero t, Numero x, Numero y ) {
 245
 246     return - D.fi( D, x, y ) * y - D.g * D.G * x / D.L;
 247
 248 }
 249
 250 void euler( Datos& D ) {
 251
 252     Numero  xo = D.xo,
 253             yo = D.yo,
 254             x  = 0.0,
 255             y  = 0.0,
 256             t  = D.to;
 257
 258     while ( t < D.tf ) {
 259
 260         // Calculo los datos para este punto.
 261         x = xo + D.k * D.fx( D, t, xo, yo );
 262         y = yo + D.k * D.fy( D, t, xo, yo );
 263
 264         // Imprimo resultados.
 265         cout << t << " " << x << " " << y << endl;
 266
 267         // Reemplazo valores iniciales.
 268         xo = x;
 269         yo = y;
 270         t += D.k;
 271
 272     }
 273
 274 }
 275
 276 void rk4( Datos& D ) {
 277
 278     Numero  x   = D.xo,
 279             y   = D.yo,
 280             t   = D.to,
 281             qx1 = 0.0,
 282             qx2 = 0.0,
 283             qx3 = 0.0,
 284             qx4 = 0.0,
 285             qy1 = 0.0,
 286             qy2 = 0.0,
 287             qy3 = 0.0,
 288             qy4 = 0.0,
 289             unSexto = 1.0 / 6.0;
 290
 291     // Imprimo datos iniciales.
 292     cout << t << " " << x << " " << y << endl;
 293
 294     while ( t < D.tf ) {
 295
 296         // Calculo los datos para este punto.
 297         qx1 = D.k * D.fx( D, t, x, y );
 298         qy1 = D.k * D.fy( D, t, x, y );
 299
 300         qx2 = D.k * D.fx( D, t + D.k / 2.0, x + qx1 / 2.0, y + qy1 / 2.0 );
 301         qy2 = D.k * D.fy( D, t + D.k / 2.0, x + qx1 / 2.0, y + qy1 / 2.0 );
 302
 303         qx3 = D.k * D.fx( D, t + D.k / 2.0, x + qx2 / 2.0, y + qy2 / 2.0 );
 304         qy3 = D.k * D.fy( D, t + D.k / 2.0, x + qx2 / 2.0, y + qy2 / 2.0 );
 305
 306         qx4 = D.k * D.fx( D, t + D.k, x + qx3, y + qy3 );
 307         qy4 = D.k * D.fy( D, t + D.k, x + qx3, y + qy3 );
 308
 309         x += unSexto * ( qx1 + 2 * qx2 + 2 * qx3 + qx4 );
 310         y += unSexto * ( qy1 + 2 * qy2 + 2 * qy3 + qy4 );
 311         t += D.k;
 312
 313         // Imprimo resultados.
 314         cout << t << " " << x << " " << y << endl;
 315
 316     }
 317
 318 }
 319
 320 void nystrom( Datos& D ) {
 321
 322     Numero  gGk2_L = D.g * D.G * D.k * D.k / D.L,
 323             xo  = D.xo,
 324             x1  = ( 1 - gGk2_L * 0.5 ) * D.xo,
 325             x2  = 0.0,
 326             y   = 0.0,
 327             fi  = 0.0,
 328             t   = D.to;
 329
 330     // Imprimo valores iniciales.
 331     cout << t << " " << xo << " " << y << endl;
 332     y = ( x1 - xo ) / D.k;
 333     cout << t << " " << x1 << " " << y << endl;
 334
 335     while ( t < D.tf ) {
 336
 337         // Calculo los datos para este punto.
 338         fi  = D.fi( D, x1, y );
 339         x2 = ( ( fi - 1 ) * xo + ( 2 - gGk2_L ) * x1 ) / ( fi + 1 );
 340
 341         // Prepara para próxima iteración
 342         t += D.k;
 343         xo = x1;
 344         x1 = x2;
 345         y  = ( x1 - xo ) / D.k;
 346
 347         // Imprimo resultados.
 348         cout << t << " " << x2 << " " << y << endl;
 349
 350     }
 351
 352 }