Код: Выделить всё
double* xv(double a,double b,double step,int np)
{
double *x,t;
x=new double [np];
t=a;
for(int i=0;i<np;i++)
{
x[i]=t;
t+=step;
}
return x;
}
double *yv(double *x,int np,double (*f)(double))
{
double *y;
y=new double [np];
for(int i=0;i<np;i++)
{
y[i]=f(x[i]);
}
return y;
}
double* gaus(int n,double **a,double *b) //методом гауса
{
double *x,t,s;
x=new double [n];
int i,j,k;
for(k=0;k<n-1;k++)
{
for(i=k+1;i<n;i++)
{
t=a[i][k]/a[k][k];
b[i]=b[i]-t*b[k];
for(j=0;j<n;j++)
a[i][j]=a[i][j]-t*a[k][j];
}
}
for(k=n-1;k>=0;k--)
{
s=0;
for(j=k+1;j<n;j++)
s+=a[k][j]*x[j];
x[k]=(b[k]-s)/a[k][k];
}
return x;
}
double proiz(int i,int m,double *x,double *y)
{
double res;
if(i==0)
res=(y[1]-y[0])/(x[1]-x[0]);
else if(i==m-1)
res=(y[m-1]-y[m-2])/(x[m-1]-x[m-2]);
else
res=(y[i+1]-y[i-1])/(x[i+1]-x[i-1]);
return res;
}
double *mit(int m,double *x,double *y)
{
int i,j;
double *mi,**a,*b,h;
a=new double *[m];
for(i=0;i<m;i++) a[i]=new double [m];
b=new double [m];
for(i=0;i<m;i++)
{
for(j=0;j<m;j++)
a[i][j]=0;
if(i==0 || i==m-1) a[i][i]=1;
else
{
a[i][i-1]=1;
a[i][i]=4;
a[i][i+1]=1;
h=x[i]-x[i-1];
if(i != 0 && i != m-1)
b[i]=3*(y[i+1]-y[i-1])/h; //система из n+1 линейных уравнений относительно неизвестных m(i)
}
if(i==0) b[i]=proiz(0,m,x,y);
else if(i==m-1) b[i]=proiz(m-1,m,x,y);
}
mi=gaus(m,a,b);
if(a)
{
for(i=0;i<m;i++) delete []a[i];
delete []a;
}
if(b) delete []b;
return mi;
}
double spline_it(double xp,int m,double *x,double *y,double *mi)
{
int i;
double sp, drob1,drob2,drob3,drob4,xstep;
for(i=1;i<m;i++)
{
if(x[i-1]<=xp && xp<=x[i])
{
xstep=x[i]-x[i-1];// Интерполяционный кубический сплайн
drob1=(xp-x[i])*(xp-x[i])*(2*(xp-x[i-1])+xstep)/pow(xstep,3);
drob2=(xp-x[i-1])*(xp-x[i-1])*(2*(x[i]-xp)+xstep)/pow(xstep,3);
drob3=(xp-x[i])*(xp-x[i])*(xp-x[i-1])/pow(xstep,2);
drob4=(xp-x[i-1])*(xp-x[i-1])*(xp-x[i])/pow(xstep,2);
sp=y[i-1]*drob1+y[i]*drob2+mi[i-1]*drob3+mi[i]*drob4;
}
}
return sp;
}
double spline(double xp,int m,double *x,double *y)
{
double *mi,sp;
mi= mit(m,x,y);
sp=spline_it(xp,m,x,y,mi);
if(mi) delete []mi;
return sp;
}