/*no se usa en el código pero una implementación a futuro podría ser implementar este tipo de interpolación
#include <cmath>
#include <iostream>
#include<iomanip>
//#include "algebralineal.h"
//#include "interpolacion.h"
//#include "gnuplot_i_calnum.h"
using namespace std;
using namespace CALNUM;
inline double fun(double x, int ft)
{
switch(ft)
{
case 0:
return 1./(1.+25*x*x);
break;
case 1:
return abs(x);
break;
case 2:
return sin(x);
case 3:
return 2*atan(x) +exp(5-x) + 5*cos(x);
}
}
template <class T> int
splines(Vector<T> &x, Vector<T> &a, Vector<T> &b, Vector<T> &c, Vector<T> &d)
{
cout<<"splines naturales"<<endl;
int n=x.size(); //n numero de puntos
int nint= n-1; //numero de intervalos
int neq =nint -1; //numero de ecuaciones
//Creacion vector diagonal superior
Vector<T> h(nint);
for( int i=0; i<nint; i++) h[i]=x[i+1]-x[i];
//Creacion Vector diagonal
Vector<T> u(nint), v(nint); //u(0), v(0) no estan definidos
for( int i=1; i<nint; i++) u[i]=2*(h[i-1]+h[i]);
for( int i=0; i<neq; i++)
v[i+1]=3*( (a[i+2]-a[i+1])/h[i+1] - (a[i+1]-a[i])/h[i] );
//Eliminacion de Gauss
for( int i=1; i<neq; i++)
{
u[i+1]=u[i+1]-h[i]*h[i]/u[i];
v[i+1]=v[i+1]-h[i]*v[i]/u[i];
}
//Sustitucion hacia atras
c[nint-1]=v[nint-1]/u[nint-1];
d[nint-1]= -c[nint-1]/(3*h[nint-1]);
for( int k=nint-2; k>0; k--)
{
c[k]=(v[k]-h[k]*c[k+1])/u[k];
d[k]=(c[k+1]-c[k])/(3*h[k]);
}
//Valores especiales en x0
c[0]= T(0);
d[0]=c[1]/(3*h[0]);
b[0]=(a[1]-a[0])/h[0]-c[1]*h[0]/3;
// Calculo de b[k]
for( int k=0; k<nint-1; k++) b[k+1]=b[k]+(c[k]+c[k+1])*h[k];
//Impresion de coeficientes
cout<< " a "<<" b "<<" c "<< " d "<<endl;
for(int i=0; i<nint; i++) cout<<setprecision(4)<<setw(8)<<a[i]<<" "<<setw(8)<<b[i]<<" "<<setw(8)<<c[i]<<" "<<setw(8)<<d[i]<<endl;
cout<<endl;
return 0;
}
template <class T> int
splines_sujetos(Vector<T> &x, Vector<T> &a, Vector<T> &b,
Vector<T> &c, Vector<T> &d, T& fp0, T& fp1)
{
cout<<"splines sujetos"<<endl;
int n=x.size(); //n numero de puntos
int nint= n-1; //numero de intervalos
// neq =nint -1; //numero de ecuaciones
//Creacion vector diagonal superior
Vector<T> h(nint);
for( int i=0; i<nint; i++) h[i]=x[i+1]-x[i];
//Creacion Vector diagonal y de terminos independientes
Vector<T> u(n), v(n);
// Valores especiales de u(0), v(0)
u[0]=2*h[0];
v[0]= 3*( (a[1]-a[0])/h[0] -fp0);
for( int i=1; i<nint; i++) u[i]=2*(h[i-1]+h[i]);
for( int i=0; i<nint-1; i++)
v[i+1]=3*( (a[i+2]-a[i+1])/h[i+1] - (a[i+1]-a[i])/h[i] );
//Valores especiales de u(n-1), v(n-1)
v[nint]=3*(fp1- (a[nint]-a[nint-1])/h[nint-1] );
u[nint]=2*h[nint-1];
//Eliminacion de Gauss
for( int i=0; i<nint; i++)
{
u[i+1]=u[i+1]-h[i]*h[i]/u[i];
v[i+1]=v[i+1]-h[i]*v[i]/u[i];
}
//Calculo de c[k] y d[k]
//Sustitucion hacia atras
c[nint]=v[nint]/u[nint];
for( int k=nint-1; k>=0; k--)
{
c[k]=(v[k]-h[k]*c[k+1])/u[k];
d[k]=(c[k+1]-c[k])/(3*h[k]);
}
//Calculo de b[k]
//Valores especiales en x0
b[0]=fp0;
for( int k=0; k<nint-1; k++) b[k+1]=b[k]+(c[k]+c[k+1])*h[k];
//Impresion de coeficientes
cout<< " a "<<" b "<<" c "<< " d "<<endl;
for(int i=0; i<nint; i++) cout<<setprecision(4)<<setw(8)<<a[i]<<" "<<
setw(8)<<b[i]<<" "<<setw(8)<<c[i]<<" "<<setw(8)<<d[i]<<endl;
return 0;
}
template <class T> T
eval_splines(T xx, const Vector<T> &x, const Vector<T> &a,
const Vector<T> &b, const Vector<T> &c, const Vector<T> &d)
{
//buscar indice k
int n=x.size();
int k=0;
int kinf=0;
int ksup=n-1;
while(ksup-kinf>1)
{
k=(ksup+kinf)/2;
if (x[k]>xx) ksup=k;
else kinf=k;
}
k=kinf;
// cout<<"xx= "<<xx<<" k= "<<k;
xx=xx-x[k];
return a[k]+ xx*( b[k]+ xx*( c[k] +d[k]*xx ) );
}
int main(){
int nd;
cout<< "Entrar numero de nodos"<<endl;
cin>>nd;
//Construccion de vectores de nodos
Vector<double> x(nd);
cout<<"Entrar valores de los nodos"<<endl;
for (int i=0; i<nd; i++) cin>>x[i];
double x0=x[0], x1=x[nd-1];
//Eleccion de la funcion a interpolar
int ft=-1;
while((ft!=0)&&(ft!=1)&&(ft!=2)&&(ft!=3))
{
cout<<"Entrar funcion: 0- Runge 1- Bernstein 2- sin(x)"
<<" 3- user "<<endl;
cin>>ft;
}
//Calculo de las ordenadas de interpolacion y de los coeficientes del
//polinomio
Vector < double > y(nd);
for (int i=0; i< nd; i++) y[i]=fun(x[i],ft);
//Derivadas en los extremos
double fp0,fp1;
//Calculo de las derivadas en los extremos
//Se puede fijar el error de la derivada por extrapolacion
double hhh=0.00001;
fp0= (fun(x0+hhh,ft)-fun(x0-hhh,ft))/(2*hhh);
fp1=(fun(x1+hhh,ft)-fun(x1-hhh,ft))/(2*hhh);
int nint= nd - 1;
Vector<double> b(nint), c(nd), d(nint);
//Impresion de los nodos
cout<<" x y"<<endl;
for(int i=0;i<nd;i++)
cout<<setprecision(4)<<setw(8)<<x[i]<<" "<<setw(8)<<y[i]<<endl;
cout<<endl<<endl;
//Splines naturales
//Calculo de los coeficientes
Vector < double > a=y;
splines<double>(x,y,b,c,d);
//evaluacion del polinomio y funcion en los puntos de dibujo (200)
int n=200;
double hh=(x1-x0)/n;
n++;
Vector <double> xx(n), yy(n) ,ff(n);
cout<<" x " << " f(x) "<< " splin(x)" <<endl;
for (int i = 0; i <n; i++)
{
xx[i]=x0+ i*hh;
yy[i] = eval_splines<double>(xx[i], x, y, b,c,d);
ff[i] = fun(xx[i],ft);
if(!(i-(i/5)*5))
cout<<setw(6)<<xx[i]<<" "<<setprecision(8)<<
setw(10)<<ff[i]<<" "<<setw(8)<<yy[i]<<endl;
}
cout<<endl;
//Dibujo del polinomio y la funcion
Gnuplot g1 = Gnuplot();
// Seleccionamos escala logaritmica para ciertas funciones
// if( (ft==0)||(ft==1)) g1.cmd("set logscale y");
// if( (ft==0)) g1.cmd("set logscale y");
g1.set_style("lines");
g1.plot_xy(xx,ff," ");
g1.set_style("points");
g1.plot_xy(xx,yy," splines naturales ");
//Splines sujetos
//Calculo de los coeficientes
splines_sujetos<double>(x,y,b,c,d, fp0,fp1);
cout<<endl<<" x " << " f(x) "<< " splin(x)" <<endl;
for (int i = 0; i <n; i++)
{
xx[i]=x0+ i*hh;
yy[i] = eval_splines<double>(xx[i], x, y, b,c,d);
ff[i] = fun(xx[i],ft);
if(!(i-(i/5)*5))
cout<<setw(6)<<xx[i]<<" "<<setprecision(8)<<
setw(10)<<ff[i]<<" "<<setw(8)<<yy[i]<<endl;
}
//Dibujo del polinomio y la funcion
Gnuplot g2 = Gnuplot();
// Seleccionamos escala logaritmica para ciertas funciones
// if( (ft==0)||(ft==1)) g1.cmd("set logscale y");
// if( (ft==0)) g1.cmd("set logscale y");
g2.set_style("lines");
g2.plot_xy(xx,ff," ");
g2.set_style("points");
g2.plot_xy(xx,yy,"splines sujetos ");
//Eleccion del tiempo de permanencia de la grafica en pantalla
int SLEEP_LGTH=15;
sleep(SLEEP_LGTH);
return 0;
}
*/