Programa Principal y Rutinas Auxiliares (comprimido)

Programa Principal.

/*
  Programa ejemplo del método de colocación asimétrica, para la
  construcción de un interpolante de la función de prueba u(x,y)=
  (5/4)+cos(5.4y)/6+6(3x-1)^2, utilizando una función radial
  capa delgada r^4 log(r).
  
  Para mayor referencia del algoritmo ver AlgInterpTPS.pdf. 

  Entrada: Número de puntos aleatorios.
  Salida: Valores de la interpolación en los puntos aleatorios. 

   Autor: Daniel A. Cervantes Cabrera. Esta dirección de correo electrónico está protegida contra spambots. Usted necesita tener Javascript activado para poder verla.
   Revisión: Pedro Gonzalez Casanova. Esta dirección de correo electrónico está protegida contra spambots. Usted necesita tener Javascript activado para poder verla.

*/

#include <stdlib.h>
#include <iostream>
#include "Interpolante.h"
#include "Matrix.h"
#include "GrammMatrix.h"
#include <time.h>
#include <fstream>
#include <math.h>
#include <string>

using namespace Interpolante;

float s(float x, float y)
{
  if (0 <= x &&  x <= 1 && 0 <= y && y <= 1)
    return (5.0/4.0+cos(5.4*y))    / (6 + 6* pow((3*x-1),2));  
  else
    return 0;  
}

float interpolante(float x, float y,
     std::vector<float> xv,std::vector<float> yv,
     std::vector<float> lambda,
     float(*tp)(float,float,float,float))
{
  int N = xv.size();
  float interp = 0.0; 

  for(int i=0; i<N; i++)
    interp += lambda[i]*(*tp)(x,y,xv[i],yv[i]);
    
  interp += lambda[N]*x*x+lambda[N+1]*y*y+
            lambda[N+2]*x*y+lambda[N+3]*x+
            lambda[N+4]*y+lambda[N+5];

  return interp;  
}

int main(int argc, char * argv[]){

  Matrix M;  
  Matrix P; 
  
  std::vector<float> u;
  std::vector<float> mx;
  std::vector<float> my;  
  
  
  if (argc <= 2){
    std::cout << "usar ./main na na.txt nm" << std::endl
       << "na = numero de nodos aleatorios" << std::endl
       << "na.txt = nombre del archivo de salida "<< std::endl
       << "nm = numero de elementos en la malla"// << std::endl 
       << "(por defecto nm = na) "<< std::endl;
       
      exit(1);
  }
 
  //Se lee el numero de puntos aleatorios.
  int N = atoi(argv[1]);  
  string narch = argv[2];
  int Nm; // Numero de elementos en la malla de evaluacion
 

  if(argc == 4)
    Nm = atoi(argv[3]);
  else
    Nm = N;

  //Se establece el tama~o en los vectores de nodos
  x.resize(N);
  y.resize(N);

  // Se establece el tama~o de la matrices.
  M.resize(N,N);
  P.resize(N,6); 

  //Se estable el tama~o de los vectores lambda y u.
  lambda.resize(N+6);
  u.resize(N+6);
  
  //Se generan los nodos aleatorios entre 0 y 1.  
  srand(time(0));  
  for (int i=0; i< N; i++){
    x[i] = 1- float(rand()) / RAND_MAX;
    y[i] = 1- float(rand()) / RAND_MAX;  
  }
  
  
  //Se crea la matriz M   
  for(int i=0; i<N; i++)
    for(int j=0; j<N; j++)
      M(i,j) = TPS(x[i],y[i],x[j],y[j]);
 

  //Se crea la matriz P   
  for(int i=0; i<N; i++){      
    P(i,0) = x[i]*x[i];
    P(i,1) = y[i]*y[i];
    P(i,2) = x[i]*y[i];
    P(i,3) = x[i];
    P(i,4) = y[i];
    P(i,5) = 1.0;      
  }

  //Se crea la matriz de Gramm
  GrammMatrix G(M,P);   

  //  std::cout << G ;
    
  
  //Se calcula el vector de soluciones exactas
  for (int i=0; i< N; i++)
    u[i] = s(x[i],y[i]);
  
  for(int i=0; i<6; i++)
    u[N+i] = 0.0;    
  
 

  //Solucion del sistema  
   lambda = G/u;

   
   float xmi,ymi;
   
   ofstream sal(narch.c_str());
   
   for(int i=0; i <= Nm; i++)
     for(int j=0; j <= Nm; j++){
       ymi = float(i)/Nm;
       xmi = float(j)/Nm;
       sal << xmi <<" "<< ymi <<"  "
   << interpolante(xmi,ymi,x,y,lambda,TPS) 
   << std::endl;       
     }        

   sal.close();

}