任务

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

1. 利用大M法计算出木条在压铁控制下的曲线，边界条件取自然边界条件；
2. 将第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$个未知量 { M i , i = 1 , 2 , ⋯ , n − 1 } . 当 M 0 = 0 , M n = 0 \left\{M_i,i=1,\right.2,\cdots,n-1\}.当M_{0}=0,M_{n}=0 {Mi​,i=1,2,⋯,n−1}.当M0​=0,Mn​=0时，称为自然边界条件.
$$\left[\begin{array}{ccccc}2& {\lambda }_1& & & \\ {\mu }_2& 2& {\lambda }_2& & \\ & \ddots & \ddots & \ddots & \\ & & {\mu }_{n-2}& 2& {\lambda }_{n-2}\\ & & & {\mu }_{n-1}& 2\end{array}\right]\left[\begin{array}{c}{M}_1\\ M_2\\ \\ M_{n-2}\\ \\ M_{n-1}\end{array}\right]=\left[\begin{array}{c}{d}_1-{\mu }_1{M}_0\\ d_2\\ ⋮\\ d_{n-2}\\ d_{n-1}-{\lambda }_{n-1}{M}_n\end{array}\right]$$
$$\begin{array}{rl}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}],x\in [x_i,x_{i+1}]\end{array}$$

代码

#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;
}