数值分析 三次样条插值及实现

 分析：

第一问，给出的是第一类边界条件

第二问，给出的是第二类边界条件

我们按照想要的步骤，分别求第一类与第二类边界条件下的三次样条插值函数即可

为了不重复计算，且易于扩展，我们用C++编程，循环实现即可。

 （这肯定不能手算的，手算必手酸）

#include <cstdio>
#include <iostream>

using namespace std;

const int N = 5;

int main()
{
    // 录入数据
    double x[N]={0.25,0.30,0.39,0.45,0.53};
    double y[N]={0.5,0.5477,0.6245,0.6708,0.7280};
    // 输出数据
    printf("X(i),i=0...%d
",N-1);
    for(int i=0;i<N;++i)
    {
        printf("%f ",x[i]);
    }
    printf("

Y(i),i=0...%d
",N-1);
    for(int i=0;i<N;++i)
    {
        printf("%f ",y[i]);
    }
    printf("

");

    // 计算h[i]
    double h[N];
    for(int i=0;i<N-1;++i)
    {
        h[i]=x[i+1]-x[i];
    }
    // 输出数据
    printf("h(i),i=0...%d
",N-2);
    for(int i=0;i<N-1;++i)
    {
        printf("%f ",h[i]);
    }
    putchar('
');
    putchar('
');

    // 计算u[i],r[i]
    double u[N],r[N];
    for(int i=1;i<N;++i)
    {
        u[i]=h[i-1]/(h[i-1] + h[i]);
        r[i]=h[i]/(h[i-1] + h[i]);
    }
    // 输出数据
    printf("u(i),i=1...%d
",N-1);
    for(int i=1;i<N;++i)
    {
        printf("%f ",u[i]);
    }
    putchar('
');
    putchar('
');
    printf("r(i),i=1...%d
",N-1);
    for(int i=1;i<N;++i)
    {
        printf("%f ",r[i]);
    }
    putchar('
');
    putchar('
');

    // 计算f[x(i-1),x(i)]
    // 用f[i]表示f[x(i-1),x(i)]
    double f[N+1];
    for(int i=1;i<N;++i)
    {
        f[i]=(y[i] - y[i-1])/(x[i] - x[i-1]);
    }
    // 输出数据
    printf("f[x(i),x(i-1)],i=1...%d
",N-1);
    for(int i=1;i<N;++i)
    {
        printf("%f ",f[i]);
    }
    putchar('
');
    putchar('
');

    // 记录两个导数
    f[0]=1;
    f[N]=0.6868;

    // 计算d(i)
    double d[N+1];
    d[0]=6*(f[1]-f[0])/h[0];
    for(int i=1;i<N-1;++i)
    {
        d[i]=6*(f[i+1]-f[i])/(h[i-1] + h[i]);
    }
    d[N]=6*(f[N]-f[N-1])/h[N-2];

    printf("d[i],i=0...%d
",N);
    for(int i=0;i<=N;++i)
    {
        printf("%f ",d[i]);
    }

    // AX=B
    // X=[A^(-1)]*B
    return 0;
}

求出上述数据后，代入矩阵求解即可得到所有的M( i ) 的值，对应矩阵的求逆与矩阵乘法的操作，用python或matlab都可以实现，这里不再重复实现，只需将上述计算所得的数据，写入相应的文件，然后在python中读取调用numpy库函数或用matlab实现都可。

所有数据都已经求出后，带入，三次样条函数即为所求

第二问给出了第二类边界条件

M₀与M_n即为给出的第二类边界条件的值，这里都是0，计算更加的简单，等式右边的直接与之前一问所求的d相同，无需重复计算

同样带入如下矩阵，

利用矩阵的逆运算，即可求出剩余的M_1,...,n-1.

同样，所有数据已知，带入，三次样条函数即为所求。