任务
point.txt文件中包含了21个压铁的位置信息
- 利用大M法计算出木条在压铁控制下的曲线,边界条件取自然边界条件;
- 将第10个压铁的位置移动至(0,10),计算出新的曲线,观察每个区间内的三次函数是否改变。
算法
μ i M i − 1 + 2 M i + λ i M i + 1 = d i , i = 1 , 2 , ⋯ , n − 1 \mu_{i}M_{i-1}+2M_{i}+\lambda_{i}M_{i+1}=d_{i},i=1,2,\cdots,n-1 μiMi−1+2Mi+λiMi+1=di,i=1,2,⋯,n−1
其中
μ i = 1 − λ i \mu_{i}=1-\lambda_{i} μi=1−λi
h i h i + h i \frac{h_{i}}{h_{i}+h_{i}} hi+hihi
d i = 6 h i + h i − 1 ( y i + 1 − y i h i − y i − y i − 1 h i − 1 ) = 6 f ( x i − 1 , x i , x i + 1 ) 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) di=hi+hi−16(hiyi+1−yi−hi−1yi−yi−1)=6f(xi−1,xi,xi+1)
给 定
M
0
,
M
n
M_{0}, M_{n}
M0,Mn的值,此时
n
−
1
n-1
n−1个方程组有
n
−
1
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时,称为自然边界条件.
[
2
λ
1
μ
2
2
λ
2
⋱
⋱
⋱
μ
n
−
2
2
λ
n
−
2
μ
n
−
1
2
]
[
M
1
M
2
M
n
−
2
M
n
−
1
]
=
[
d
1
−
μ
1
M
0
d
2
⋮
d
n
−
2
d
n
−
1
−
λ
n
−
1
M
n
]
\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}
2μ2λ12⋱λ2⋱μn−2⋱2μn−1λn−22
M1M2Mn−2Mn−1
=
d1−μ1M0d2⋮dn−2dn−1−λn−1Mn
S
(
x
)
=
(
x
i
+
1
−
x
)
3
M
i
+
(
x
−
x
i
)
3
M
i
+
1
6
h
i
+
(
x
i
+
1
−
x
)
y
i
+
(
x
−
x
i
)
y
i
+
1
h
i
−
h
i
6
[
(
x
i
+
1
−
x
)
M
i
+
(
x
−
x
i
)
M
i
+
1
]
,
x
∈
[
x
i
,
x
i
+
1
]
\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}
S(x)=6hi(xi+1−x)3Mi+(x−xi)3Mi+1+hi(xi+1−x)yi+(x−xi)yi+1−6hi[(xi+1−x)Mi+(x−xi)Mi+1],x∈[xi,xi+1]
代码
#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;
}