Кубический сплайн

void CubDay(double h, double *y, double dcof[4]) // * Построение кубического многочлена по 4-м точкам *
{  double a, z0, z1, z2, h2;
 
    a = dcof[0] = y[1];
    z0 = y[0] - a;
    z1 = y[2] - a;
    z2 = y[3] - a;
    dcof[1] = ( - z2 + 6*z1 - 2*z0 ) / (6*h) ;
    h2 = h*h;
    dcof[2] = (  z1 + z0) / (2*h2) ;
    dcof[3] = (  z2 - 3*z1 - z0) / (6*h*h2) ;
}

#ifndef SPLINE_H
#define SPLINE_H
#include <cmath>
 
 
//Модуль дробного числа
 double AbsD(const double& x)
{
    double k = 1;
    if (x < 0) k = -1;
    return x*k;
}
 
//Численное интегрирование методом Симпсона
//Вертает числовое значение функции заданной таблично y[]
//x[] - массив значений аргумента с равным шагом h=x[i+1]-x[i]
//y[] - массив значений функции
//n   - число узлов матрицы (Обязательно чётное число)
static double NumIntegralMethSimpson(const double *x, const double *y, const int& n)
{
    double yout = 0;
    double a = x[0];
    double b = x[n-1];
 
    double y0 = 0;
    double y1 = 0;
    double y2 = 0;
    double y3 = 0;
 
    y0 = y[0];
    y3 = y[n-1];
 
    for (int i = 1; i < n-1; i = i+2)
    {
        y1 = y1 + y[i];
    }
 
 
    for (int i = 2; i < n-2; i = i+2)
    {
        y2 = y2 + y[i];
    }
 
    y1 = y1*4;
    y2 = y2 + y2;
 
    yout = (b-a)*(y0 + y1 + y2 + y3)/(6*n);
 
 
    return yout;
}
 
//Заполнение массивов производных
void inic_cubicSpline(const double *xi, const double *yi, double *m, const int& ni)
{
    //Заполняем работчие массивы
 
    double d = 0;
    double e = 0;
    double h = 0;
    double f = 0;
    double p = 0;
 
    //Создаём работчие массивы
    double *l=new double[ni];
    double *r=new double[ni];
    double *s=new double[ni];
 
    d = xi[1]-xi[0];
    e = yi[1]-yi[0];
    e = e/d;
 
    for (int k = 1; k <= ni-2; k++)
    {
        h = d;
        d = xi[k+1]-xi[k];
        f = e;
        e = yi[k+1]-yi[k];
        e = e/d;
        l[k] = d/(d+h);
        r[k] = 1-l[k];
        s[k] = 6*(e-f)/(h+d);
    }
 
    for (int k = 1; k <= ni-2; k++)
    {
        p = 1/(r[k]*l[k-1]+2);
        l[k] = -l[k]*p;
        s[k] = (s[k]-r[k]*s[k-1])*p;
    }
 
    m[ni-1] = 0;
    l[ni-2] = s[ni-2];
    m[ni-2] = l[ni-2];
 
 
    for (int k=ni-3;k>=0;k--)
    {
        l[k] = l[k]*l[k+1]+s[k];
        m[k] = l[k];
    }
}
 
 
//Расчёт значения функции для аргументта x
static double Calc_cubicSpline(const double *xi, const double *yi,
                               const double *m, const int& ni, const double& x)
{
    double yout = 0;
 
    double d = 0;
    double h = 0;
    double p = 0;
    double rr = 0;
 
    int i = 0;
    int j = 0;
 
    if (x <= xi[0])
    {
        d = xi[1]-xi[0];
        yout = -d*m[1]/6+(yi[1]-yi[0])/d;
        yout = yout*(x-xi[0])+yi[0];
 
    }
    else
    {
        i = 0;
 
        while(x > xi[i]) i++;
        j = i - 1;
 
        d = xi[i]-xi[j];
        h = x-xi[j];
        rr = xi[i]-x;
        p = d*d/6;
        yout = (m[j]*rr*rr*rr+m[i]*h*h*h)/6/d;
        yout = yout + ((yi[j]-m[j]*p)*rr+(yi[i]-m[i]*p)*h)/d;
    }
 
    return yout;
}
 
//Сплайн итнерполяция
//Функция реализующая кубическиий сплайн.
//Вертает значение функции для соответствующего x в массиве y
//и максимальное значение y по модулю
//xi[]- входной массив значений xi контрольных точек
//yi[]- входной массив значений yi контрольных точек
//ni  - количество контрольных точек
//x[] - масссив значений аргумента для вычисления функции
//y[] - выходной массив значений функций
//n   - количество точек для расчёта
static double cubicSpline(const double *xi, const double *yi, const int &ni,
                          const double *x, double *y, const int& n)
{
    double yout = 0;
    double xinp = 0;
 
    double ymaxout = 0;
 
    //Создаём работчие массивы
    double *m=new double[ni];
 
    int k = 0;
 
    //Заполняем раьотчие массивы
    inic_cubicSpline(xi,yi,m,ni);
 
    //Расчитываем значения аргументов для x[]
    for(k = 0;k <= n; k++)
    {
       xinp = x[k];
       yout = Calc_cubicSpline(xi,yi,m,ni,xinp);
       y[k] = yout;
 
        if (ymaxout < AbsD(yout)) ymaxout =  AbsD(yout);
    }
    return ymaxout;
}
 
 
//Функция реализующая кубическиий сплайн.
//Вертает значение функции для соответствующего x
//и максимальное значение y по модулю
//xi[]- входной массив значений xi контрольных точек
//yi[]- входной массив значений yi контрольных точек
//ni  - количество контрольных точек
//x   - аргумент для которого необходимо расчитать значение функции
static double cubicSplineX(const double *xi, const double *yi, const int& ni, const double& x)
{
    double yout = 0;
    double xinp = x;;
 
    //Заполняем работчие массивы
    double *m=new double[ni];
 
    inic_cubicSpline(xi,yi,m,ni);
 
    //Расчитываем значения аргументов для x
    xinp = x;    
    yout = Calc_cubicSpline(xi,yi,m,ni,xinp);
 
    return yout;
}
 
//Функция реализующая определение производной кубического сплайна.
//Вертает значение производной функции для соответствующего x
//и максимальное значение y по модулю
//xi[]- входной массив значений xi контрольных точек
//yi[]- входной массив значений yi контрольных точек
//ni  - количество контрольных точек
//nx  - dx = x/nx
static double DiffCubicSplineX(const double *xi, const double *yi, const int& ni,
                               const double x, const int& nx)
{
    double yout1 = 0;
    double yout2 = 0;
    double xinp  = x;
    double dx    = xinp/(double)nx;
 
    //Заполняем работчие массивы
    double *m=new double[ni];
 
    inic_cubicSpline(xi,yi,m,ni);
 
    //Расчитываем значения аргументов для x
    xinp = x;
    yout1 = Calc_cubicSpline(xi,yi,m,ni,xinp);
    //Расчитываем значения аргументов для x
    xinp = x+dx;
    yout2 = Calc_cubicSpline(xi,yi,m,ni,xinp);
 
    return (yout2-yout1)/dx;
}
 
 
//Функция реализующая кубическиий сплайн функции 2-х переменных.
//Вертает значение функции для соответствующих x и y
//и максимальное значение y по модулю
//xi[]  - входной массив значений xi контрольных точек
//yi[]  - входной массив значений yi контрольных точек
//yi[][]- входной массив значений yi контрольных точек
//xn    - количество контрольных точек по оси x
//yn    - количество контрольных точек по оси y
//x,y   - аргументы для которых необходимо расчитать значение функции
static double cubicSplineZ_XY(const double *xi, const double *yi,
                              const double **zi, const int& xn, int yn,
                              const double& x, const double& y)
{
    int i = 0; //Ось X
    int j = 0; //Ось Y
 
    double yout = 0;
 
    double *y_x=new double[yn];
    double *x_in=new double[xn];
 
    //Заполняем столбец для требуемого x
    for (i = 0; i < yn; i++)
    {
        for (j = 0; j <xn; j++)
        {
            x_in[j] = zi[j][i];
        }
 
        y_x[i] = cubicSplineX(xi,x_in,xn,x);
    }
 
    yout = cubicSplineX(yi,y_x,yn,y);
 
    return yout;
}
 
#endif // SPLINE_H

#include <math.h>
#include <iostream>
 
const int n = 10;
const double a = 0;
const double b = 1;
const double h = 0.1; //h = (b-a)/n
const double param = 0.3;
const double arr_p_inv[] = {0.0000001, 0.000276, 0.000021, 0.00014, 0.00000001, 0.0001083, 0.00018, 0.9, 0.002425, 0.00013, 0.000000001};
const double arr_rand[] = {0.0001221, -0.002012, 0.000231, -0.0001234, -0.0051021, -0.0055103, 0.0001528, -0.00521, -0.0051342, -0.008132, 0.000001};
 
using namespace std;
 
double f (double x)
{
    double res = cos(7*x)*exp(pow(x,3)) + pow(x,2);
    return res;
}
 
double dfdx (double x)
{
    double res = (-7*sin(7*x) + 3*pow(x,2)*cos(7*x))*exp(pow(x,3)) + 2*x;
    return res;
}
 
double d2fdx2 (double x)
{
    double sin_x = sin(7*x);
    double cos_x = cos(7*x);
    double res = (-49*cos_x - 42*pow(x,2)*sin_x + 6*x*cos_x + 9*pow(x,4)*cos_x)*exp(pow(x,3)) + 2;
    return res;
}
 
int def_elem_interval (double x)
{
    for (int i = 1; i <= n; i++)
        if ((x >= a + (i-1)*h)&&(x <= a + i*h))
            return i;
    return -1;
}
 
double* get_moments_cub ()
{
    double* moments = (double*) malloc ((n+1)*sizeof(double));
    //a,b,c,k,g
    double A = h/6;
    double B = 2*h/3;
    double C = h/6;
    double K = 0;
    double* F_knots = (double*) malloc ((n+1)*sizeof(double));
    for (int i = 0; i < n+1; i++)
        F_knots[i] = f(a + i*h);
    double* G = (double*) malloc ((n+1)*sizeof(double));
    G[0] = 0;
    for (int i = 1; i < n; i++)
        G[i] = (F_knots[i+1] - 2*F_knots[i] + F_knots[i-1])/h;
    G[n] = 0;
    //direct
    double* alpha = (double*) malloc ((n+1)*sizeof(double));
    double* beta = (double*) malloc ((n+1)*sizeof(double));
    alpha[0] = 0;
    beta[0] = 0;
    alpha[1] = K;
    beta[1] = G[0];
    for (int i = 1; i <= n-1; i++) {
        alpha[i+1] = - C/(A*alpha[i] + B);
        beta[i+1] = (G[i] - A*beta[i])/(A*alpha[i] + B);
    }
    moments[n] = (G[n] + K*beta[n])/(1 - K*alpha[n]);
    for (int i = n-1; i >= 0; i--)
        moments[i] = alpha[i+1]*moments[i+1] + beta[i+1];
    return moments;
}
 
double* disc_get_moments_cub ()
{
    double* moments = (double*) malloc ((n+1)*sizeof(double));
    //a,b,c,k,g
    double A = h/6;
    double B = 2*h/3;
    double C = h/6;
    double K = 0;
    double* F_knots = (double*) malloc ((n+1)*sizeof(double));
    for (int i = 0; i < n+1; i++)
        F_knots[i] = f(a + i*h) + arr_rand[i];
    double* G = (double*) malloc ((n+1)*sizeof(double));
    G[0] = 0;
    for (int i = 1; i < n; i++)
        G[i] = (F_knots[i+1] - 2*F_knots[i] + F_knots[i-1])/h;
    G[n] = 0;
    //direct
    double* alpha = (double*) malloc ((n+1)*sizeof(double));
    double* beta = (double*) malloc ((n+1)*sizeof(double));
    alpha[0] = 0;
    beta[0] = 0;
    alpha[1] = K;
    beta[1] = G[0];
    for (int i = 1; i <= n-1; i++) {
        alpha[i+1] = - C/(A*alpha[i] + B);
        beta[i+1] = (G[i] - A*beta[i])/(A*alpha[i] + B);
    }
    moments[n] = (G[n] + K*beta[n])/(1 - K*alpha[n]);
    for (int i = n-1; i >= 0; i--)
        moments[i] = alpha[i+1]*moments[i+1] + beta[i+1];
    return moments;
}
 
double cubic_spl (double x)
{
    double* moments = get_moments_cub();
    int i = def_elem_interval(x);
    double m_1 = moments[i-1];
    double m_2 = moments[i];
    double x_1 = a + (i-1)*h;
    double x_2 = a + i*h;
    double f_1 = f(x_1);
    double f_2 = f(x_2);
    double res = m_1*pow(x_2 - x, 3)/(6*h) + m_2*pow(x - x_1, 3)/(6*h) + ((f_2-f_1)/h - (m_2-m_1)*h/6)*(x-x_2) + f_2 - m_2*pow(h, 2)/6;
    return res;
}
 
double disc_cubic_spl (double x)
{
    double* moments = disc_get_moments_cub();
    int i = def_elem_interval(x);
    double m_1 = moments[i-1];
    double m_2 = moments[i];
    double x_1 = a + (i-1)*h;
    double x_2 = a + i*h;
    double f_1 = f(x_1) + arr_rand[i-1];
    double f_2 = f(x_2) + arr_rand[i];
    double res = m_1*pow(x_2 - x, 3)/(6*h) + m_2*pow(x - x_1, 3)/(6*h) + ((f_2-f_1)/h - (m_2-m_1)*h/6)*(x-x_2) + f_2 - m_2*pow(h, 2)/6;
    return res;
}
 
void print_cubic_spl()
{
    cout << "        Cubic spline" << endl;
    cout << "  x        f(x)       s(x)       |(f-s)(x)|" << endl;
    double max1 = 0;
    for (int i = 0; i < n; i++) {
        double x = a + i*h + param*h;
        double f_ = f(x);
        double s = cubic_spl(x);
        cout << x << "    " << f_ << "    " << s << "    " << fabs(f_-s) << endl;
        if (fabs(f_-s) > max1)
            max1 = fabs(f_-s);
    }
    cout << "Max deviation equals " << max1 << endl;
    cout << "        Cubic spline in the case when f has a discrepancy" << endl;
    cout << "  x        f(x)       s(x)       |(f-s)(x)|" << endl;
    double max2 = 0;
    for (int i = 0; i < n; i++) {
        double x = a + i*h + param*h;
        double f_ = f(x);
        double s = disc_cubic_spl(x);
        cout << x << "    " << f_ << "    " << s << "    " << fabs(f_-s) << endl;
        if (fabs(f_-s) > max2)
            max2 = fabs(f_-s);
    }
    cout << "Max deviation equals " << max2 << endl;
}
 
int getmax_abs_i(int k, double** arr)
{
    int max_i = k;
    double max_a = fabs(arr[k][k]);     
    for (int i = k; i < n-1; i++)
        for (int j = k; j < n-1; j++)
            if (fabs(arr[i][j]) > max_a) {
                max_a = fabs(arr[i][j]);
                max_i = i;
            }
    return max_i;
}
 
int getmax_abs_j(int k, double** arr)
{
    int max_j = k;
    double max_a = fabs(arr[k][k]);        
    for (int i = k; i < n-1; i++)
        for (int j = k; j < n-1; j++)
            if (fabs(arr[i][j]) > max_a) {
                max_a = fabs(arr[i][j]);
                max_j = j;
            }
    return max_j;
}
 
double* Gauss_method (double** arr_a, double* vect)
{
    int* order;
    order = (int*)malloc((n-1)*sizeof(int));
    for (int i = 0; i < n-1; i++)
        order[i] = i;
    double** arr_a_k_minus_1;
    arr_a_k_minus_1 = (double**)malloc((n-1)*sizeof(double*));
    for (int i = 0; i < n-1; i++)
        arr_a_k_minus_1[i] = (double*)malloc(n*sizeof(double));
    for (int i = 0; i < n-1; i++)
        for (int j = 0; j <= n-1; j++)
            arr_a_k_minus_1[i][j] = (j==(n-1))?vect[i]:arr_a[i][j];
    for (int k = 0; k < n-1; k++) {
        int ik = getmax_abs_i(k, arr_a_k_minus_1);
        int jk = getmax_abs_j(k, arr_a_k_minus_1);
        if (ik != k) {
            for (int j = k; j <= n-1; j++) {
                double temp = arr_a_k_minus_1[k][j];
                arr_a_k_minus_1[k][j] = arr_a_k_minus_1[ik][j];
                arr_a_k_minus_1[ik][j] = temp;
            }
        }
        if (jk != k) {
            for (int i = 0; i < n-1; i++) {
                double temp = arr_a_k_minus_1[i][k];
                arr_a_k_minus_1[i][k] = arr_a_k_minus_1[i][jk];
                arr_a_k_minus_1[i][jk] = temp;
            }
            int temp = order[k];
            order[k] = order[jk];
            order[jk] = temp;
        }
        double diag = 1/arr_a_k_minus_1[k][k];
        for (int j = k; j <= n; j++) {
            arr_a_k_minus_1[k][j] *= diag;
        }
        for (int i = 0; i < n-1; i++)
            if (i != k) {
                double temp = arr_a_k_minus_1[i][k];
                for (int j = k; j <= n-1; j++)
                    arr_a_k_minus_1[i][j] -= temp*arr_a_k_minus_1[k][j];
            }
    }
    double* x;
    x = (double*) malloc ((n-1)*sizeof(double));
    for (int i = 0; i < n-1; i++) { 
        x[order[i]] = arr_a_k_minus_1[i][n-1];
    }
    return x;
}
 
double** build_matr()
{
    double** matr;
    matr = (double**) malloc((n-1)*sizeof(double*));
    for (int i = 0; i < n-1; i++)
        matr[i] = (double*) malloc((n-1)*sizeof(double));
    for (int i = 0; i < n-1; i++)
        for (int j = 0; j < n-1; j++) {
            matr[i][j] = 0;
            if (i == j)
                matr[i][j] += 2*h/3;
            if ((i-j == 1) || (j-i == 1))
                matr[i][j] += h/6;
            double A = 0, B = 0, C = 0;
            if (i == j-2)
                A = 1/h;
            if (i == j-1) {
                A = -2/h;
                B = 1/h;
            }
            if (i == j) {
                A = 1/h;
                B = -2/h;
                C = 1/h;
            }
            if (i == j+1) {
                B = 1/h;
                C = -2/h;
            }
            if (i == j+2)
                C = 1/h;
            matr[i][j] += (arr_p_inv[j]*A - 2*arr_p_inv[j+1]*B + arr_p_inv[j+2]*C)/h;
        }
    return matr;
}
 
double* build_hf ()
{
    double* f_ = (double*) malloc ((n+1)*sizeof(double));
    for (int i = 0; i <= n; i++)
        f_[i] = f(a + i*h) + arr_rand[i];
    double* hf = (double*) malloc ((n-1)*sizeof(double));
    for (int i = 0; i < n-1; i++) 
        hf[i] = (f_[i] - 2*f_[i+1] + f_[i+2])/h;
    return hf;
}
 
double* get_moments_smooth()
{
    double* moments = (double*) malloc ((n-1)*sizeof(double));
    double** matr = build_matr();
    double* hf = build_hf();
    moments = Gauss_method (matr, hf);
    return moments;
}
 
double* get_min_G(double* moments)
{
    double* y = (double*) malloc ((n+1)*sizeof(double));
    for (int i = 0; i <= n; i++) {
        y[i] = f(a + i*h) + pow(h, 10) + arr_p_inv[i]*moments[i]/h + arr_rand[i];
        if (i > 0)
            y[i] -= 2*arr_p_inv[i]*moments[i-1]/h;
        if (i > 1)
            y[i] += arr_p_inv[i]*moments[i-2]/h;
    }
    return y;
}
 
double smoothing_spl (double x, double* moments, double* y)
{
    int i = def_elem_interval(x);
    double m_1 = moments[i-1];
    double m_2 = moments[i];
    double x_1 = a + (i-1)*h;
    double x_2 = a + i*h;
    double y_1 = y[i-1];
    double y_2 = y[i];
    double res = m_1*pow(x_2 - x, 3)/(6*h) + m_2*pow(x - x_1, 3)/(6*h) + (y_1 - pow(h, 2)*m_1/6)*(x_2 - x)/h + (y_2 - pow(h, 2)*m_2/6)*(x - x_1)/h;
    return res;
}
 
double smoothing_spl (double x)
{
    double* moments = get_moments_smooth();
    double* y = get_min_G(moments);
    double res = smoothing_spl (x, moments, y);
    return res;
}
 
void print_smooth_spl()
{
    cout << "        Smoothing spline" << endl;
    cout << "  x        f(x)       s(x)       |(f-s)(x)|" << endl;
    double max = 0;
    for (int i = 0; i < n; i++) {
        double x = a + i*h + param*h;
        double f_ = f(x);
        double s = smoothing_spl(x);
        cout << x << "    " << f_ << "    " << s << "    " << fabs(f_-s) << endl;
        if (fabs(f_-s) > max)
            max = fabs(f_-s);
    }
    cout << "Max deviation equals " << max << endl;
}
 
int _main()
{
    print_cubic_spl();
    print_smooth_spl();
    return 0;
}