計算方法實驗7：實現三次樣條插值算法

任務

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

本文來自互聯網用戶投稿，該文觀點僅代表作者本人，不代表本站立場。本站僅提供信息存儲空間服務，不擁有所有權，不承擔相關法律責任。
如若轉載，請注明出處：http://www.pswp.cn/web/11427.shtml
繁體地址，請注明出處：http://hk.pswp.cn/web/11427.shtml
英文地址，請注明出處：http://en.pswp.cn/web/11427.shtml

如若內容造成侵權/違法違規/事實不符，請聯系多彩編程網進行投訴反饋email:809451989@qq.com，一經查實，立即刪除！