數值分析 三次樣條插值及實現

 分析：

第一問，給出的是第一類邊界條件

第二問，給出的是第二類邊界條件

我們按照想要的步驟，分別求第一類與第二類邊界條件下的三次樣條插值函數即可

為了不重復計算，且易於擴展，我們用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",N-1);
    for(int i=0;i<N;++i)
    {
        printf("%f ",x[i]);
    }
    printf("\n\nY(i),i=0...%d\n",N-1);
    for(int i=0;i<N;++i)
    {
        printf("%f ",y[i]);
    }
    printf("\n\n");

    // 計算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",N-2);
    for(int i=0;i<N-1;++i)
    {
        printf("%f ",h[i]);
    }
    putchar('\n');
    putchar('\n');

    // 計算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",N-1);
    for(int i=1;i<N;++i)
    {
        printf("%f ",u[i]);
    }
    putchar('\n');
    putchar('\n');
    printf("r(i),i=1...%d\n",N-1);
    for(int i=1;i<N;++i)
    {
        printf("%f ",r[i]);
    }
    putchar('\n');
    putchar('\n');

    // 計算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",N-1);
    for(int i=1;i<N;++i)
    {
        printf("%f ",f[i]);
    }
    putchar('\n');
    putchar('\n');

    // 記錄兩個導數
    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",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.

同樣，所有數據已知，帶入，三次樣條函數即為所求。