77891.cpp

   1 // vim: set tabstop=4 softtabstop=4 shiftwidth=4 expandtab:
   2 //
   3 // Trabajo Práctico I de Análisis Numérico I
   4 // Este programa resuelve un sistema de ecuaciones lineales ralo
   5 // resultante de la discretización mediante diferencias finitas
   6 // de ecuaciones diferenciales.
   7 // Copyright (C) 2002 Leandro Lucarella <leandro@lucarella.com.ar>
   8 //
   9 // Este programa es Software Libre; usted puede redistribuirlo
  10 // y/o modificarlo bajo los términos de la "GNU General Public
  11 // License" como lo publica la "FSF Free Software Foundation",
  12 // o (a su elección) de cualquier versión posterior.
  13 //
  14 // Este programa es distribuido con la esperanza de que le será
  15 // útil, pero SIN NINGUNA GARANTIA; incluso sin la garantía
  16 // implícita por el MERCADEO o EJERCICIO DE ALGUN PROPOSITO en
  17 // particular. Vea la "GNU General Public License" para más
  18 // detalles.
  19 //
  20 // Usted debe haber recibido una copia de la "GNU General Public
  21 // License" junto con este programa, si no, escriba a la "FSF
  22 // Free Software Foundation, Inc.", 59 Temple Place - Suite 330,
  23 // Boston, MA  02111-1307, USA.
  24 //
  25 // $URL$
  26 // $Date$
  27 // $Rev$
  28 // $Author$
  29 //
  30
  31 #include <stdio.h>
  32 #include <stdlib.h>
  33 #include <math.h>
  34
  35 // Tipos de datos.
  36 typedef unsigned int Indice;
  37 typedef float Numero;
  38 typedef struct {
  39     Indice  N;
  40     Numero  w;
  41     Numero  rtol;
  42     Numero  TS;
  43     Numero  TO;
  44     Numero  TN;
  45     Numero  TE;
  46     Numero* X;
  47 } Datos;
  48
  49 // Constantes.
  50 const Numero DEFAULT_W    = 1.0,
  51              DEFAULT_RTOL = 0.01,
  52              DEFAULT_T    = 1.0;
  53
  54 // Calcula la iteración para un nodo.
  55 Numero procesarNodo( Datos&, Indice );
  56
  57 // Calcula la iteración para un nodo dependiendo del tipo de nodo.
  58 Numero procesarNodo( Datos&, Indice, bool, bool, bool, bool, Numero = 0 );
  59
  60 // Calcula la norma 2 de un vector.
  61 Numero norma2( Numero*, Indice );
  62
  63 // Imprime un vector por la salida estándar.
  64 void imprimirVector( Indice, Numero* );
  65
  66 // Imprime un vector por la salida estándar.
  67 void imprimirVector( Numero*, Indice );
  68
  69 // Resuelve el sistema de ecuaciones lineales por el método S.O.R.
  70 // (o Gasuss-Seidel si w=1).
  71 void resolver( Datos& );
  72
  73 int main( int argc, char* argv[] ) {
  74
  75     // Se fija que tenga los argumentos necesarios para correr.
  76     if ( argc < 2 ) {
  77         printf( "Faltan argumentos. Modo de uso:\n" );
  78         printf( "\t%s N w RTOL TS TO TN TE\n", argv[0] );
  79         printf( "Desde w en adelante son opcionales. Los valores por defecto son:\n" );
  80         printf( "\tw                 = %f\n", DEFAULT_W );
  81         printf( "\tRTOL              = %f\n", DEFAULT_RTOL );
  82         printf( "\tTS = TO = TN = TE = %f\n", DEFAULT_T );
  83         return EXIT_FAILURE;
  84     }
  85
  86     // Se inicializan los datos.
  87     Datos D;
  88     D.N      = Indice( atol( argv[1] ) );
  89     if ( D.N < 3 ) {
  90         printf( "N Debe ser un entero mayor o igual a 3\n" );
  91         return EXIT_FAILURE;
  92     }
  93     D.w    = DEFAULT_W;
  94     D.rtol = DEFAULT_RTOL;
  95     D.TS   =
  96     D.TO   =
  97     D.TN   =
  98     D.TE   = DEFAULT_T;
  99     // Se aloca uno más de lo necesario para usar el rango [1,n]
 100     D.X    = new Numero[Indice(pow(D.N-1,2))+1];
 101
 102     // Si se pasaron datos como argumento, se los va agregando.
 103     switch ( argc ) {
 104         case 8: D.TE    = Numero( atof( argv[7] ) );
 105         case 7: D.TN    = Numero( atof( argv[6] ) );
 106         case 6: D.TO    = Numero( atof( argv[5] ) );
 107         case 5: D.TS    = Numero( atof( argv[4] ) );
 108         case 4: D.rtol = Numero( atof( argv[3] ) );
 109         case 3: D.w     = Numero( atof( argv[2] ) );
 110     }
 111
 112     resolver( D );
 113
 114     delete D.X;
 115
 116     return EXIT_SUCCESS;
 117
 118 }
 119
 120 Numero procesarNodo( Datos& D, Indice i ) {
 121
 122     Indice N = D.N;
 123
 124     // Si es el primero, es SO.
 125     if ( i == 1 )
 126         return procesarNodo( D, i,
 127                 false, false, true, true,
 128                 D.TS + D.TO );
 129
 130     // Si es N - 1 es NO.
 131     if ( i == ( N - 1 ) )
 132         return procesarNodo( D, i,
 133                 false, true, false, true,
 134                 D.TN + D.TO );
 135
 136     // Si es N² - 3N + 3 es SE.
 137     if ( i == ( N*N - 3*N + 3 ) )
 138         return procesarNodo( D, i,
 139                 true, false, true, false,
 140                 D.TS + D.TE );
 141
 142     // Si es N² - 2N + 1 es NE.
 143     if ( i == ( N*N - 2*N + 1 ) )
 144         return procesarNodo( D, i,
 145                 true, true, false, false,
 146                 D.TN + D.TE );
 147
 148     // Si pertenece al intervalo [2;N-1] es O.
 149     if ( ( 1 < i ) && ( i < ( N - 1 ) ) )
 150         return procesarNodo( D, i,
 151                 false, true, true, true,
 152                 D.TO );
 153
 154     // Si pertenece al intervalo [N² - 3N + 3;N² - 2N + 1] es E.
 155     if ( ( ( N*N - 3*N + 3 ) < i ) && ( i < ( N*N - 2*N + 1 ) ) )
 156         return procesarNodo( D, i,
 157                 true, true, true, false,
 158                 D.TE );
 159
 160     // Si i - 1 es múltiplo de N - 1 (y no es NE, SE, NO, SO) es S.
 161     if ( ( ( i - 1 ) % ( N - 1 ) ) == 0 )
 162         return procesarNodo( D, i,
 163                 true, false, true, true,
 164                 D.TS );
 165
 166     // Si i es múltiplo de N - 1 (y no es NE, SE, NO, SO) es N.
 167     if ( ( i % ( N - 1 ) ) == 0 )
 168         return procesarNodo( D, i,
 169                 true, true, false, true,
 170                 D.TN );
 171
 172     return procesarNodo( D, i, true, true, true, true );
 173
 174 }
 175
 176 Numero procesarNodo( Datos& D, Indice i,
 177                      bool i_n1, bool i_1, bool i1, bool in_1,
 178                      Numero b ) {
 179
 180     Numero x  = 0,
 181            Xo = D.X[i];
 182
 183     // Voy sumando los elementos que corresponden.
 184     if ( i_n1 )
 185         x += D.X[i-D.N+1];
 186     if ( i_1 )
 187         x += D.X[i-1];
 188     if ( i1 )
 189         x += D.X[i+1];
 190     if ( in_1 )
 191         x += D.X[i+D.N-1];
 192
 193     // Calculo el elemento i.
 194     D.X[i] = ( 1 - D.w ) * D.X[i] + ( b + x ) * D.w / 4;
 195
 196     // Devuelve parte de la sumatoria para calcular la norma 2.
 197     return pow( D.X[i] - Xo, 2 );
 198
 199 }
 200
 201 Numero norma2( Numero* X, Indice n ) {
 202
 203     Numero sum = 0;
 204
 205     for ( Indice i = 1; i <= n; i++ )
 206         sum += pow( X[i], 2 );
 207
 208     return sqrt( sum );
 209
 210 }
 211
 212 void imprimirVector( Indice n, Numero* X ) {
 213
 214     for ( Indice i = 1; i <= n; i++ )
 215         printf( "%.7e%s", X[i], ( i == n ) ? "\n" : " " );
 216
 217 }
 218
 219 void resolver( Datos& D ) {
 220
 221     Indice ite  = 0,
 222            n    = Indice( pow( D.N - 1, 2) );
 223     Numero r    = D.rtol + 1,
 224            Ek   = 0,
 225            Ek_1 = 0;
 226
 227     // Imprime datos iniciales.
 228     // Formato de salida:
 229     // iteración R S <vector X>
 230     // (si no se puede calcular S, se muestra 0)
 231     printf( "%d %.2e 0.00000000 ", ite, r );
 232     imprimirVector( n, D.X );
 233
 234     while ( r > D.rtol ) {
 235
 236         for ( Indice i = 1; i <= n; i++ )
 237             Ek += procesarNodo( D, i );
 238         Ek = sqrt( Ek );
 239         r = Ek / norma2( D.X, n );
 240         printf( "%d %.2e %.8f ", ++ite, r, Ek_1 ? ( Ek / Ek_1 ) : 0 );
 241         imprimirVector( n, D.X );
 242         Ek_1 = Ek;
 243         Ek   = 0;
 244
 245     }
 246
 247 }