NumericalMethods/lab9/lab9.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef struct{
    double a, b, c, d;
    double x;
}CubicSpline;    

double *EquallySpaced(double a, double b, int node){   
    int i;
    double step;
    step = (b-a)/(node-1);
    double* x = (double*)malloc(node*sizeof(double));
    for (i=0; i<node; i++){
        x[i] = a + i*step;
    } 
    return x;
}

void FindCubicSplines(CubicSpline* splines, int n, double *xk, double *yk, double dfa, double dfb)
{
    double* h = (double*)malloc((n-1)*sizeof(double));
    
    double* b = (double*)malloc(n * sizeof(double));

    double* c = (double*)malloc(n * sizeof(double));

    double* d = (double*)malloc(n * sizeof(double));

    double* f = (double*)malloc(n * sizeof(double));

    double* M = (double*)malloc(n * sizeof(double));

    int i;
    b[0] = 0;
    d[n-1] = 0;
    for (i = 0; i < n - 1; i++)
    {
        h[i] = xk[i+1] - xk[i];
    }
    c[0] = -h[0] / 3;
    d[0] = -h[0] / 6;
    b[n-1] = h[n-2] / 6;
    c[n-1] = h[n-2] / 3;
    f[0] = dfa - (yk[1]-yk[0]) / h[0];
    f[n-1] = dfb - (yk[n-1]-yk[n-2]) / h[n-2];

    for(i = 1; i < n - 1; i++)
    {
        b[i] = h[i-1] / (h[i-1] + h[i]);
        c[i] = 2;
        d[i] = h[i] / (h[i-1] + h[i]);
        f[i] = (6 / (h[i-1] + h[i])) * ((yk[i+1] - yk[i]) / h[i] - (yk[i] - yk[i-1]) / h[i-1]);
    } 

    //Найти Mi
    for(i = 1; i < n; i++)
    {
        c[i] = c[i] - b[i] * d[i-1] / c[i-1];
        f[i] = f[i] - b[i] * f[i-1]/c[i-1];
        b[i] = 0;        
    }

    M[n-1] = f[n-1] / c[n-1];
    for(i = n - 2; i >= 0; i--)
    {
        M[i] = (f[i] - d[i] * M[i+1]) / c[i];
    }
    
    //Найти коэффициенты
    for(i = 0; i < n - 1; i++){
        splines[i].a = M[i] / (6*h[i]);
        splines[i].b = M[i+1] / (6*h[i]);
        splines[i].c = (yk[i+1] - yk[i]) / h[i] - h[i] * (M[i+1]-M[i]) / 6;
        splines[i].d = yk[i] - M[i] * h[i] * h[i] / 6;
        splines[i].x = xk[i];
    } 
    splines[n-1].x = xk[n-1];
    free(h);
    free(b);
    free(c);
    free(d);
    free(f);
    free(M);
}

double cubicSplineInterpolate(CubicSpline* splines, int n, double x) {
    int i;
    double dx1, dx2;
    for (i = 0; i < n - 1; i++) {

        if (x >= splines[i].x && x <= splines[i+1].x) {
            dx1 = x - splines[i].x;
            dx2 = splines[i+1].x - x;
            return splines[i].a*dx2*dx2*dx2 + splines[i].b*dx1*dx1*dx1 + splines[i].c*dx1 + splines[i].d;
        }

    }
    return 0; 
}
double Func1(double x) {
    return 1.0*(exp(x) - exp(-x))/(exp(x)+exp(-x));
}

double Func2(double x) {
    double sign_x;
    if (x>0) sign_x =1;
    else if (x==0) sign_x =0;
    else sign_x = -1;

    return 3 * sign_x * pow(x*1.0, 4) + 4 * pow(x*1.0, 3) - 12 * pow(x*1.0, 2) + 4;
}

int main() { 
    FILE *f = fopen("data.txt", "w");
    if (f == NULL) {
        printf("Error opening file!\n");
        return 1;
    }

    int nodes[] = {5,7,10}; 
    int i,j,n;
    double dfa, dfb, a, b, x, y;

    dfa = 48;
    dfb = 31.5;
    a = -1; b = 1.5;

    for(i=0; i<3;i++){
        n = nodes[i];
        double *xk;
        CubicSpline* splines = (CubicSpline*)malloc(n * sizeof(CubicSpline));        
        double* yk = (double*)malloc(n*sizeof(double));
        
        fprintf(f,"%d ", n);
        xk = EquallySpaced(a, b, n);
        for(j =0; j< n; j++) yk[j] = Func2(xk[j]);
        FindCubicSplines(splines, n, xk, yk, dfa, dfb);
        

        for(j=0; j<1000;j++){
            x = a+j*(b-a)/999;
            y = cubicSplineInterpolate(splines, n, x);
            fprintf(f, "%.13f ", y);
        }
        fprintf(f,"\n");
        free(splines);
        free(xk);
        free(yk);
    }
    fclose(f);
    return 0;
}