计算方法实验7：实现三次样条插值算法

千里澄江

已于 2024-05-12 12:27:41 修改

阅读量901

点赞数 32

分类专栏： C/C++ 计算方法实验文章标签：算法线性代数矩阵 c++

于 2024-05-12 12:24:52 首次发布

本文链接：https://blog.csdn.net/weixin_73004416/article/details/138751326

版权

C/C++ 同时被 2 个专栏收录

30 篇文章 1 订阅

订阅专栏

计算方法实验

10 篇文章 0 订阅

订阅专栏

任务

point.txt文件中包含了21个压铁的位置信息

利用大M法计算出木条在压铁控制下的曲线，边界条件取自然边界条件；
将第10个压铁的位置移动至(0,10)，计算出新的曲线，观察每个区间内的三次函数是否改变。

算法

$\mu_{i}M_{i-1}+2M_{i}+\lambda_{i}M_{i+1}=d_{i},i=1,2,\cdots,n-1$

其中

$\mu_{i}=1-\lambda_{i}$

$\frac{h_{i}}{h_{i}+h_{i}}$

$d_{i}=\frac{6}{h_{i}+h_{i-1}}\left({\frac{y_{i+1}-y_{i}}{h_{i}}-\frac{y_{i}-y_{i-1}}{h_{i-1}}}\right)=6f\left(x_{i-1},x_{i},x_{i+1}\right)$

给定 $M_{0}, M_{n}$ 的值，此时 $n - 1$ 个方程组有 $n - 1$ 个未知量 $\left\{M_i,i=1,\right.2,\cdots,n-1\}.当M_{0}=0,M_{n}=0$ 时，称为自然边界条件.
$\begin{bmatrix}2&\lambda_{1}\\\mu_{2}&2&\lambda_{2}\\&\ddots&\ddots&\ddots\\&&\mu_{n-2}&2&\lambda_{n-2}\\&&&\mu_{n-1}&2\end{bmatrix}\begin{bmatrix}M_{1}\\M_{2}\\\\M_{n-2}\\\\M_{n-1}\end{bmatrix}=\begin{bmatrix}d_{1}-\mu_{1}M_{0}\\d_{2}\\\vdots\\d_{n-2}\\d_{n-1}-\lambda_{n-1}M_{n}\end{bmatrix}$
$\begin{aligned}S\left(x\right)=&\frac{\left(x_{i+1}-x\right)^{3}M_{i}+\left(x-x_{i}\right)^{3}M_{i+1}}{6h_{i}}+\frac{\left(x_{i+1}-x\right)y_{i}+\left(x-x_{i}\right)y_{i+1}}{h_{i}}\\&-\frac{h_{i}}{6}[(x_{i+1}-x)M_{i}+(x-x_{i})M_{i+1}],\quad x\in[x_{i},x_{i+1}]\end{aligned}$

代码

#include<bits/stdc++.h>
using namespace std;

//三次样条插值
vector<long double> Spline_Interpolation(const vector<long double>& x, const vector<long double>& y);
// 列主元消去法求解线性方程组
vector<long double> Column_Elimination(vector<vector<long double>> A, vector<long double> b);

int main()
{
    ifstream file("point.txt");
    if (!file) {
        cerr << "Unable to open file point.txt";
        return 1;
    }

    vector<long double> x, y, lambda, d, h1, miu, h2;
    long double a, b;
    while (file >> a >> b)
    {
        x.push_back(a);
        y.push_back(b);
        lambda.push_back(0);
        d.push_back(0);
        h1.push_back(0);
        h2.push_back(0);
        miu.push_back(0);
    }
    
    cout << "原始数据插值结果：" << endl;
    vector<long double> M1 = Spline_Interpolation(x, y);
    vector<vector<long double>> S1(M1.size()-1, vector<long double>(4, 0));
    for(int i = 0; i < M1.size(); i++)
    {
        cout << "M1[" << i << "]:" << M1[i] << endl;
        h1[i] = x[i + 1] - x[i];
    }    
    cout << endl;
    for(int i = 0; i < M1.size() - 1; i++)
    {
        cout << "S1[" << i << "]:";
        S1[i][0] = (M1[i + 1] - M1[i]) / (6 * h1[i]);
        cout << S1[i][0] << "x^3";
        S1[i][1] = (x[i + 1] * M1[i] - x[i] * M1[i + 1]) / (2 * h1[i]);
        if(S1[i][1] >= 0)
            cout << "+" ;
        cout << S1[i][1] << "x^2";
        S1[i][2] = (3 * M1[i+1] * x[i] * x[i] - 3 * M1[i] * x[i + 1] * x[i + 1] - 6 * y[i] + 6 * y[i+1] + h1[i] * h1[i] * M1[i] - h1[i] * h1[i] * M1[i + 1]) / (6 * h1[i]);
        if(S1[i][2] >= 0)
            cout << "+" ;
        cout << S1[i][2] << "x";
        S1[i][3] = (x[i + 1] * x[i + 1] * x[i + 1] * M1[i] - x[i] * x[i] * x[i] * M1[i + 1] + 6 * x[i+1] * y[i] - 6 * x[i] * y[i+1] - h1[i] * h1[i] * M1[i] * x[i+1] + h1[i] * h1[i] * M1[i+1] * x[i]) / (6 * h1[i]);
        if(S1[i][3] >= 0)
            cout << "+" ;
        cout << S1[i][3];
        cout << endl;
    }

    //修改第十个压铁的坐标为(0,10)
    cout << endl << "修改第十个压铁的坐标为(0,10)后：" << endl;
    y[9] = 10;
    vector<long double> M2 = Spline_Interpolation(x, y);
    vector<vector<long double>> S2(M2.size()-1, vector<long double>(4, 0));
    for(int i = 0; i < M2.size(); i++)
    {
        cout << "M2[" << i << "]:" << M2[i] << endl;
        h2[i] = x[i + 1] - x[i];
    }
    cout << endl;
    for(int i = 0; i < M2.size() - 1; i++)
    {
        cout << "S2[" << i << "]:";
        S2[i][0] = (M2[i + 1] - M2[i]) / (6 * h1[i]);
        cout << S2[i][0] << "x^3";
        S2[i][1] = (x[i + 1] * M1[i] - x[i] * M2[i + 1]) / (2 * h1[i]);
        if(S2[i][1] >= 0)
            cout << "+" ;
        cout << S2[i][1] << "x^2";
        S2[i][2] = (3 * M2[i+1] * x[i] * x[i] - 3 * M2[i] * x[i + 1] * x[i + 1] - 6 * y[i] + 6 * y[i+1] + h2[i] * h2[i] * M2[i] - h2[i] * h2[i] * M2[i + 1]) / (6 * h2[i]);
        if(S2[i][2] >= 0)
            cout << "+" ;
        cout << S2[i][2] << "x";
        S2[i][3] = (x[i + 1] * x[i + 1] * x[i + 1] * M2[i] - x[i] * x[i] * x[i] * M2[i + 1] + 6 * x[i+1] * y[i] - 6 * x[i] * y[i+1] - h2[i] * h2[i] * M2[i] * x[i+1] + h2[i] * h2[i] * M2[i+1] * x[i]) / (6 * h2[i]);
        if(S2[i][3] >= 0)
            cout << "+" ;
        cout << S2[i][3];
        cout << endl;
    }
    return 0;
}

//三次样条插值
vector<long double> Spline_Interpolation(const vector<long double>& x, const vector<long double>& y) 
{
    int len = x.size();
    int n = len - 1;
    vector<long double> h(n), lambda(n), miu(n), d(n);
    for(int i = 0; i < n; i++)
        h[i] = x[i + 1] - x[i];
    for(int i = 1; i < n; i++)
    {
        lambda[i] = h[i] / (h[i] + h[i - 1]);
        miu[i] = 1 - lambda[i];
        d[i] = 6 / (h[i] + h[i - 1]) * ((y[i + 1] - y[i]) / h[i] - (y[i] - y[i - 1]) / h[i - 1]);
    }
    vector<vector<long double>> A(n - 1, vector<long double>(n - 1, 0));
    for(int i = 0; i < n - 1; i++)
    {
        A[i][i] = 2;
        if(i != 0)
            A[i][i - 1] = miu[i + 1];
        if(i != n - 2)
            A[i][i + 1] = lambda[i + 1];
    }
    vector<long double> B(n - 1, 0);
    for(int i = 0; i < n - 1; i++)
        B[i] = d[i + 1];
    vector<long double> M = Column_Elimination(A, B);
    return M;
}

// 列主元消去法求解线性方程组
vector<long double> Column_Elimination(vector<vector<long double>> A, vector<long double> b)
{
    int n = A.size();
    vector<long double> x(n + 2, 0);
    vector<vector<long double>> matrix(n, vector<long double>(n + 1, 0));

    for(int i = 0; i < n; i++)
    {
        matrix[i][n] = b[i];
        for(int j = 0; j < n; j++)
            matrix[i][j] = A[i][j];
    }
        
    for(int col = 0; col < n; col++)//找到列主元
    {
        long double maxnum = abs(matrix[col][col]);
        int maxrow = col;
        for(int row = col + 1; row < n; row++)
        {
            if(abs(matrix[row][col]) > maxnum)
            {
                maxnum = abs(matrix[row][col]);
                maxrow = row;
            }
        }
        swap(matrix[col], matrix[maxrow]);
        for(int row = col + 1; row < n; row++)
        {
            long double res = matrix[row][col] / matrix[col][col];
            for(int loc = col; loc <= n; loc++)
                matrix[row][loc] -= matrix[col][loc] * res; 
        }
    }
    for(int row = n - 1; row >= 0; row--)//回代
    {
        for(int col = row + 1; col < n; col++)
        {
            matrix[row][n] -= matrix[col][n] * matrix[row][col] / matrix[col][col];
            matrix[row][col] = 0;
        }
        matrix[row][n] /= matrix[row][row];
        matrix[row][row] = 1;
    }
    for(int i = 0; i < n; i++)
        x[i+1] = matrix[i][n];
    return x;
}