三次样条插值函数（边界二）分析及代码zorkey的博客-

三次样条插值函数（边界二）分析及代码
根据三次样条插值公式与题目中所给出的X Y值，分别计算所需要的数值。

1. 步长：hi=xi+1-xi (i=0，1，…，n-1)，

2. μ值：μ[k] = h[k – 1] / (h[k – 1] + h[k])

3. λ值 ：λ[k] = 1 – μ[k];

4. 计算一阶差商，二阶差商。

5. 根据边界条件二与所求的的二阶差商计算出C，C=6*二阶差商，（c[i] = (6 / (h[i – 1] + h[i]) * ((y[i + 1] – y[i]) / h[i] – (y[i] – y[i – 1]) / h[i – 1]))。

6. 再利用λ和µ以及C用追赶法计算出矩阵的解M。

7. 对于所求出的 μ, λ,列成矩阵A[][],为初始矩阵。

8. 利用追赶法，对初始矩阵A[][]进行LU分解，并进行计算 。

9. s″(x0)=f0″， s″(xn)=fn″，利用边界二进行求解。

10. 若第2种端点条件取为：M0=Mn=0(s″(x0)=s″(xn)=0)

11. 最后，利用求解出来的M[i],计算是s（x）；
    求解过程
    1． 设置step()函数，利用for循环，计算出步长h, μ值，λ值；
    2． 设置connect（）函数，计算C；
    3． 设置det（）函数，初始化矩阵A；
    4． 设置run()函数，用追赶法计算求得M;
    5． 设置get()函数，根据s（x）公式，分别计算出每个段的插值。

程序代码

在这里插入代码片 #include<iostream> #include<math.h> using namespace std; const int max = 19; double x[max], y[max];//已知 自变量 因变量 double h[max], b[max], m[max], n[max];//步长 常量2  μ λ double c[max], Y[max], M[max];//解向量  追赶中间的Y 求得的解 double A[max][max] = { 0 }, l[max][max] = { 0 }, u[max][max];//追赶法矩阵 double  s; //s(x)的结果 void step(double x[]) { for (int i = 0; i < max - 1; i++)   b[i] = 2; for (int i = 0; i < max - 1; i++) {   h[i] = x[i + 1] - x[i]; //cout << h[i] << " "; } for (int k = 1; k < max - 1; k++) {   m[k] = h[k - 1] / (h[k - 1] + h[k]);//μ值    n[k] = 1 - m[k];//λ值 //cout << m[k] <<  " "; } } void conect(double y[]) {  c[0] = 0;  c[max - 1] = 0; for (int i = 1; i < max - 1; i++) {   c[i] = (6 / (h[i - 1] + h[i]) * ((y[i + 1] - y[i]) / h[i] - (y[i] - y[i - 1]) / h[i - 1])); //cout << c[i] << " "; } } /*void connect(double y[], double x[]){               //求一阶导二阶导  for (int i = 0; i <= max - 2; i++)   fir[i] = (y[i + 1] - y[i]) / h[i];   for (int i = 0; i <= max - 3; i++)   sec[i] = (fir[i + 1] - fir[i]) / (x[i + 2] - x[i]);  for (int i = 1; i <= max - 2; i++)  {   c[i] = 6 * sec[i - 1];   cout << c[i] << " ";  } }*/ void det() { //初始化 for (int i = 0; i < max; i++) for (int j = 0; j < max; j++) if (i == j)//对角线为2 {     A[i][i] = b[i]; } else if (j - i == 1) {     A[i][j] = n[i + 1]; } else if (i - j == 1) {     A[i][j] = m[i + 2]; } else //其余为0     A[i][j] = 0; } //追赶法计算线性方程组 void run() { for (int i = 0; i < max - 2; i++) { for (int j = 0; j < max - 2; j++) { if (i == j) //计算l，u矩阵中间对角线部分 {     l[i][j] = 1; if (i == 0) {      u[i][j] = A[i][j] = 2; } else {      u[i][j] = A[i][j] - l[i][j - 1] * A[i - 1][j]; } } else if (j == i - 1) //中间对角线下面的对角线 {     l[i][j] = A[i][j] / u[i - 1][j];     u[i][j] = 0; } else if (j == i + 1) //中间对角线上面对角线 {     u[i][j] = A[i][j];     l[i][j] = 0; }    Else    { //其余部分都是0     u[i][j] = 0;     l[i][j] = 0; } } }  Y[0] = c[1]; for (int i = 1; i < max - 2; i++) {   Y[i] = c[i + 1] - l[i][i - 1] * Y[i - 1]; } for (int i = 0; i < max; i++)   cout << "Y:" << Y[i] << endl;  M[0] = M[max - 1] = 0;//边界二条件  M[max - 2] = Y[max - 3] / u[max - 3][max - 3]; for (int i = max - 3; i > 0; i--) {   M[i] = (Y[i - 1] - u[i - 1][i] * M[i + 1]) / u[i - 1][i - 1]; } for (int i = 0; i < max; i++) {   cout << "M[i]" << M[i] << endl; } } double get(double h[max], double M[max], double x[max], double y[max], double l) { double s = 0; for (int i = 0; i < max; i++) { if (l > x[i] && l < x[i + 1]) {    s = (M[i] * (pow((x[i + 1] - l), 3))) / (6 * h[i]) + (M[i + 1] * (pow((l - x[i]), 3))) / (6 * h[i]) + (y[i] - (M[i] * h[i] * h[i] / 6)) * (x[i + 1] - l) / h[i] + (y[i + 1] - (M[i + 1] * h[i] * h[i] / 6)) * (l - x[i]) / h[i]; } } return s; } int main() { double x[19] = { 0.520,3.1,8.0,17.95,28.65,39.62,50.65,78,104.6,156.6,208.6,260.7,312.5,364.4,416.3,468,494,507,520 };//自变量 double y[19] = { 5.288,9.4,13.84,20.20,24.90,28.44,31.10,35,36.9,36.6,34.6,31.0,26.34,20.9,14.8,7.8,3.7,1.5,0.2 };//因变量 step(x); conect(y); det(); run();  cout << endl;  cout << "所求的边界二结果为" << endl << endl;  cout << "x=2.3，S(x)=" << get(h, M, x, y, 2.3) << endl;  cout << "x=130，S(x)=" << get(h, M, x, y, 130) << endl;  cout << "x=350，S(x)=" << get(h, M, x, y, 350) << endl;  cout << "x=515，S(x)=" << get(h, M, x, y, 515) << endl; //sx(x); } MATLAB程序 function y=yangtiao() x0=[0.52,3.1,8.0,17.95,28.65,39.62,50.65,78,104.6,156.6,208.6,260.7,312.5,364.4,416.3,468,494,507,520]; y0= [5.288,9.4,13.84,20.20,24.90,28.44,31.10,35,36.9,36.6,34.6,31.0,26.34,20.9,14.8,7.8,3.7,1.5,0.2]; %pp=spline(x0,y0) %pp.coefs y=spline(x0,y0) plot(x0,y0)