Método implicito, TPS + Pol de grado 6
*
Source:    cd_imp_tps.cpp
                  Conveccion - Difusion en 2-d 
                  Metodo implicito 
                  Nucleo de placa delgada + Pol de grado 6 
Compile:
                   c++ -o cditps cd_imp_tps.cpp -lm 

Example:
                   ./cditps 30 1 0.1 10 1   0.001 0.1 
Output:
-----------------------------------------------------
Convection-Difussion - Implicit - TPS
n           = 30 , total(nxn) =900
a           = 1.000000
b           = 0.100000
convectivo  = 10.000000
difusivo    = 1.000000
stept       = 0.001000
tmax        = 0.100000

t=0.100000, e_rms       = 0.000676
t=0.100000, e_inf       = 0.002434
-----------------------------------------------------

Los siguientes archivos de datos que
se pueden ver con GNUPLOT
          
            splot  'real'
            splot  'numerica'
            splot  'realnumerica'
            splot  'realt'
            splot  'numericat'
            splot  'realnumericat'

real   -> corresponde a la solucion analitica de la EDP en t=0
realt  -> corresponde a la solucion analitica de la EDP en t=tmax
numerica   -> corresponde a la solucion numerica de la EDP en t=0
numericat   -> corresponde a la solucion numerica de la EDP en t=tmax
realnumerica ->     analitica - numerica  en t=0
realnumericat ->     analitica - numerica  en t=tmax

Autor: Jose Antonio Muñoz Gómez
Supervisión: Pedro González Casanova


*/
#include
#include
#include

#define INTERIOR       0
#define FRONTERA_REAL  1

#include "../libs/jag_mem.cpp"
#include "../libs/jag_phi.cpp"
#include "../libs/jag_vector.cpp"
#include "../libs/jag_cond.cpp"
#include "../libs/qSort.cpp"
#include "../libs/jag_file.cpp"

#include "cd_edp_ut.cpp"


#include "../libs/varios1.cpp"



//--------------------------------------------------------------------
void reconstruye(double *x,double *y, double *lambda,double c,int N, double *sol_numerica)
{
int     i,j;
double xc,yc;
double s;

 for(i=0;i      \n\n");
   exit(1);     
 }
 
    n = atoi(argv[1]);
a_cd = atof(argv[2]);  
b_cd = atof(argv[3]);  
Convectivo = atof(argv[4]);  
Difusivo       = atof(argv[5]);  
stept = atof(argv[6]);  
tmax  = atof(argv[7]);  

 tini = 0.0;
 Pol  = 6;
 printf("-----------------------------------------------------\n");
 printf("n           = %d , total(nxn) =%d\n",n,n*n);
 printf("a           = %f\n",a_cd);   
 printf("b           = %f\n",b_cd);   
 printf("convectivo  = %f\n",Convectivo);   
 printf("difusivo    = %f\n",Difusivo);   
 printf("stept       = %f\n",stept);   
 printf("tmax        = %f\n",tmax);    


//Creo el grid cartesiano
 xmin =  0;
 ymin =  0;
 xmax =  1;
 ymax =  1;

  malla  = make_mesh(xmin,xmax,ymin,ymax,n);
  
  N = malla.N;
  x = malla.x;
  y = malla.y;
tipo= malla.tipo;  
 ni = malla.ni; 
 
 
//Reservo la memoria dinamica
  A           = make_dmatrix(N+Pol,N+Pol);
  M           = make_dmatrix(N+Pol,N+Pol);
  P           = make_dmatrix(N+Pol,N+Pol);
  f           = make_dvector(N+Pol);
  lambda      = make_dvector(N+Pol);
  pivot_index = make_ivector(N+Pol);
  
  sol_analitica          = make_dvector(N);
  sol_numerica           = make_dvector(N);
  sol_analitica_numerica = make_dvector(N);   
   
  
//Creo la matriz de interpolacion:  condicion inicial

//Relleno de zeros
  fill(f,N+Pol,0.0);


 for(int i=0;i<(N+Pol);i++)
   fill(A[i],N+Pol,0.0);


//Creo la matriz de interpolacion  
for(int i=0;i