数值分析作业-三次样条插值

python三次样条插值函数一、什么是插值函数插值函数是一种数学方法，用于通过给定数据点之间的间隔来估计未知数据点的值。

在Python中，我们可以使用三次样条插值函数来进行这样的估计。

二、三次样条插值三次样条插值是一种数值分析方法，用于在给定数据点之间构造一个平滑的多项式函数。

这个函数被称为样条函数，由许多小的多项式片段组成。

在每个数据点之间，这些多项式片段满足一定的条件，使得整个函数是连续且光滑的。

2.1 样条函数的性质三次样条插值函数具有以下性质： - 在每个数据点处，函数值等于给定的数据点的函数值。

- 在每个数据点处，函数的一阶导数值等于给定数据点的一阶导数值。

- 在每个数据点处，函数的二阶导数值等于给定数据点的二阶导数值。

- 在数据点之间，函数是一个三次多项式。

2.2 插值函数的构造要构造三次样条插值函数，我们需要以下步骤： 1. 首先，给定一些数据点，这些数据点包含要插值的函数的值。

2. 然后，计算每个数据点之间的插值多项式的系数。

3. 接下来，定义一个样条函数，它由这些插值多项式组成。

4. 最后，使用这个样条函数来估计未知数据点的值。

三、三次样条插值函数的Python实现在Python中，我们可以使用SciPy库中的interp1d函数来实现三次样条插值。

interp1d函数接受一维数组作为输入，并返回一个能够进行插值的函数对象。

3.1 安装SciPy库要使用interp1d函数，首先需要安装SciPy库。

可以使用以下命令来安装SciPy：pip install scipy3.2 使用interp1d函数进行插值以下是使用interp1d函数进行三次样条插值的示例代码：import numpy as npfrom scipy.interpolate import interp1d# 定义一些数据点x = np.array([1, 2, 3, 4, 5])y = np.array([2, 3, 5, 8, 9])# 使用interp1d函数进行插值f = interp1d(x, y, kind='cubic')# 估计新的数据点的值x_new = np.array([1.5, 2.5, 3.5, 4.5])y_new = f(x_new)print(y_new)以上代码中，我们首先定义了一些数据点，然后使用interp1d函数创建了一个插值函数对象f。

《数值分析课程设计-三次样条插值》报告掌握三次样条插值函数的构造方法，体会三次样条插值函数对被逼近函数的近似。

三次样条插值函数边界条件由实际问题对三次样条插值在端点的状态要求给出。

以第1 边界条件为例，用节点处二阶导数表示三次样条插值函数，用追赶法求解相关方程组。

通过Matlab 编制三次样条函数的通用程序，可直接显示各区间段三次样条函数体表达式，计算出已给点插值并显示各区间分段曲线图。

引言分段低次样条插值虽然计算简单、稳定性好、收敛性有保证且易在电子计算机上实现，但只能保证各小段曲线在连接处的连续性，不能保证整件曲线的光滑性。

利用样条插值，既可保持分段低次插值多项式，又可提高插值函数光滑性。

故给出分段三次样条插值的构造过程算法步骤，利用Matlab软件编写三次样条插值函数通用程序，并通过数值算例证明程序的正确性。

三次样条函数的定义及特征定义：设[a,b] 上有插值节点，a＝x1＜x2<…xn＝b，对应函数值为y1,y2,⋯yn。

若函数S(x) 满足S(xj) ＝ yj （ j ＝ 1,2, ⋯,n ）， S(x) 在[xj,xj+1] （ j ＝1,2,⋯,n-1）上都是不高于三的多项式（为了与其对应j 从1 开始，在Matlab 中元素脚标从1 开始）。

当S(x) 在 [a,b] 具有二阶连续导数。

则称S(x) 为三次样条插值函数。

要求S(x) 只需在每个子区间[xj,xj+1] 上确定 1 个三次多项式，设为：Sj(x)＝ajx3＋bjx2＋cjx＋dj, (j＝1,2,⋯,n-1) (1)其中aj,bj,cj,dj 待定，并要使它满足：S(xj)＝yj, S(xj-0)=S(xj＋0), (j=2,⋯,n-1) (2)S'(xj-0)=S'(xj＋0), S"(xj-0)=S"(xj＋0), (j=2,⋯,n-1) (3)式(2)、(3)共给出n＋3(n-2)＝4n-6 个条件，需要待定4(n-1) 个系数，因此要唯一确定三次插值函数，还要附加2个边界条件。

三次样条插值的方法和思路摘要：1.三次样条插值的基本概念2.三次样条插值的数学原理3.三次样条插值的实现步骤4.三次样条插值的优缺点5.三次样条插值在实际应用中的案例正文：在日常的科学研究和工程应用中，我们经常会遇到需要对一组数据进行插值的问题。

插值方法有很多，其中三次样条插值是一种常见且有效的方法。

本文将从基本概念、数学原理、实现步骤、优缺点以及实际应用案例等方面，全面介绍三次样条插值的方法和思路。

一、三次样条插值的基本概念三次样条插值（Cubic Spline Interpolation）是一种基于分段多项式的插值方法。

它通过在各个节点上构建一条三次多项式曲线，使得这条曲线在节点之间满足插值条件，从而达到拟合数据的目的。

二、三次样条插值的数学原理三次样条插值的数学原理可以分为两个部分：一是分段三次多项式的构建，二是插值条件的满足。

1.分段三次多项式的构建假设有一组数据点序列为（x0，y0），（x1，y1），（x2，y2），（x3，y3），我们可以将这些数据点连接起来，构建一条分段三次多项式曲线。

分段三次多项式在每个子区间上都是一个三次多项式，它们之间通过节点值进行连接。

2.插值条件的满足为了使分段三次多项式在节点之间满足插值条件，我们需要在每个子区间上满足以下四个条件：（1）端点条件：三次多项式在区间的端点上分别等于节点值；（2）二阶导数条件：三次多项式在区间内的二阶导数等于节点间的斜率；（3）三阶导数条件：三次多项式在区间内的三阶导数等于节点间的曲率；（4）内部点条件：三次多项式在区间内部满足插值函数的连续性。

通过求解这四个条件，我们可以得到分段三次多项式的系数，从而实现插值。

三、三次样条插值的实现步骤1.确定插值节点：根据数据点的位置，选取合适的节点；2.构建分段三次多项式：根据节点值和插值条件，求解分段三次多项式的系数；3.计算插值结果：将待插值点的横坐标代入分段三次多项式，得到插值结果。

/* 三次样条插值计算算法*/#include "math.h "#include "stdio.h "#include "stdlib.h "/*N：已知节点数N+1R：欲求插值点数R+1x，y为给定函数f(x)的节点值{x(i)} (x(i) <x(i+1)) ,以及相应的函数值{f(i)} 0 <=i <=NP0=f(x0)的二阶导数；Pn=f(xn)的二阶导数u:存插值点{u(i)} 0 <=i <=R求得的结果s(ui)放入s[R+1] 0 <=i <=R返回0表示成功，1表示失败*/int SPL(int N,int R,double x[],double y[],double P0,double Pn,double u[],double s[]){/*声明局部变量*/double *h; /*存放步长:{hi} 0 <=i <=N-1 */double *a; /*存放系数矩阵{ai} 1 <=i <=N ；分量0没有利用*/ double *c; /*先存放系数矩阵{ci} 后存放{Bi} 0 <=i <=N-1 */double *g; /*先存放方程组右端项{gi} 后存放求解中间结果{yi} 0 <=i <=N */double *af; /*存放系数矩阵{a(f)i} 1 <=i <=N ；*/double *ba; /*存放中间结果0 <=i <=N-1*/double *m; /*存放方程组的解{m(i)} 0 <=i <=N ；*/int i,k;double p1,p2,p3,p4;/*分配空间*/if(!(h=(double*)malloc(N*sizeof(double)))) exit(1);if(!(a=(double*)malloc((N+1)*sizeof(double)))) exit(1);if(!(c=(double*)malloc(N*sizeof(double)))) exit(1);if(!(g=(double*)malloc((N+1)*sizeof(double)))) exit(1);if(!(af=(double*)malloc((N+1)*sizeof(double)))) exit(1);if(!(ba=(double*)malloc((N)*sizeof(double)))) exit(1);if(!(m=(double*)malloc((N+1)*sizeof(double)))) exit(1);/*第一步：计算方程组的系数*/for(k=0;k <N;k++)h[k]=x[k+1]-x[k];for(k=1;k <N;k++)a[k]=h[k]/(h[k]+h[k-1]);for(k=1;k <N;k++)c[k]=1-a[k];for(k=1;k <N;k++)g[k]=3*(c[k]*(y[k+1]-y[k])/h[k]+a[k]*(y[k]-y[k-1])/h[k-1]); c[0]=a[N]=1;g[0]=3*(y[1]-y[0])/h[0]-P0*h[0]/2;g[N]=3*(y[N]-y[N-1])/h[N-1]+Pn*h[N-1]/2;/*第二步：用追赶法解方程组求{m(i)} */ba[0]=c[0]/2;g[0]=g[0]/2;for(i=1;i <N;i++){af[i]=2-a[i]*ba[i-1];g[i]=(g[i]-a[i]*g[i-1])/af[i];ba[i]=c[i]/af[i];}af[N]=2-a[N]*ba[N-1];g[N]=(g[N]-a[N]*g[N-1])/af[N];m[N]=g[N]; /*P110 公式：6.32*/ for(i=N-1;i> =0;i--)m[i]=g[i]-ba[i]*m[i+1];/*第三步：求值*/for(i=0;i <=R;i++){/*判断u(i)属于哪一个子区间，即确定k */if(u[i] <x[0] || u[i]> x[N]){/*释放空间*/free(h);free(a);free(c);free(g);free(af);free(ba);free(m);return 1;}k=0;while(u[i]> x[k+1])k++;//p1=(h[k]+2*(u[i]-x[k])*pow((u[i]-x[k+1]),2)*y[k])/pow(h[k],3); //p2=(h[k]-2*(u[i]-x[k+1])*pow((u[i]-x[k]),2)*y[k+1])/pow(h[k],3);p1=(h[k]+2*(u[i]-x[k]))*pow((u[i]-x[k+1]),2)*y[k]/pow(h[k],3);p2=(h[k]-2*(u[i]-x[k+1]))*pow((u[i]-x[k]),2)*y[k+1]/pow(h[k],3); p3=(u[i]-x[k])*pow((u[i]-x[k+1]),2)*m[k]/pow(h[k],2);p4=(u[i]-x[k+1])*pow((u[i]-x[k]),2)*m[k+1]/pow(h[k],2);s[i]=p1+p2+p3+p4;}/*释放空间*/free(h);free(a);free(c);free(g);free(af);free(ba);free(m);return 0;}void main(){int N,R;double *x,*y,*u,*s;double P0,Pn;int i;/*验证算法：*/N=7;R=6;/*分配空间*/if(!(x=(double*)malloc((N+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}if(!(y=(double*)malloc((N+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}if(!(u=(double*)malloc((R+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}if(!(s=(double*)malloc((R+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}x[0]=0.5;x[1]=0.7;x[2]=0.9;x[3]=1.1;x[4]=1.3;x[5]=1.5;x[6]=1.7;x[7]=1.9;y[0]=0.4794;y[1]=0.6442;y[2]=0.7833;y[3]=0.8912;y[4]=0.9636;y[5]=0.9975;y[6]=0.9917;y[7]=0.9 463;u[0]=0.6;u[1]=0.8;u[2]=1.0;u[3]=1.2;u[4]=1.4;u[5]=1.6;u[6]=1.8;P0=-0.4794;Pn=-0.9463;if(!SPL( N, R, x, y, P0, Pn, u, s)){/*打印结果*/printf( "\nx= ");for(i=0;i <=N;i++)printf( "%8.1f ",x[i]);printf( "\ny= ");for(i=0;i <=N;i++)printf( "%8.4f ",y[i]);printf( "\n\nu= ");for(i=0;i <=R;i++)printf( "%9.2f ",u[i]);printf( "\ns= ");for(i=0;i <=R;i++)printf( "%9.5f ",s[i]);printf( "\nsin= ");for(i=0;i <=R;i++)printf( "%9.5f ",sin(u[i]));}/*释放空间*/free(x);free(y);free(u);free(s);}/* 测试数据来自课本55页例5 《数值分析》清华大学出版社第四版*/ //输入327.7 4.128 4.329 4.130 3.013.0 -4.0//输出输出三次样条插值函数：1: [27.7 , 28]13.07*(x - 28)^3 + 0.22*(x - 27.7)^3+ 14.84*(28 - x) + 14.31*(x - 27.7)2: [28 , 29]0.066*(29 - x)^3 + 0.1383*(x - 28)^3+ 4.234*(29 - x) + 3.962*(x - 28)3: [29 , 30]0.1383*(30 - x)^3 - 1.519*(x - 29)^3+ 3.962*(30 - x) + 4.519*(x - 29)//三次样条插值函数#include<iostream>#include<iomanip>using namespace std;const int MAX = 50;float x[MAX], y[MAX], h[MAX];float c[MAX], a[MAX], fxym[MAX];float f(int x1, int x2, int x3){float a = (y[x3] - y[x2]) / (x[x3] - x[x2]);float b = (y[x2] - y[x1]) / (x[x2] - x[x1]);return (a - b)/(x[x3] - x[x1]);} //求差分void cal_m(int n){ //用追赶法求解出弯矩向量M……float B[MAX];B[0] = c[0] / 2;for(int i = 1; i < n; i++)B[i] = c[i] / (2 - a[i]*B[i-1]);fxym[0] = fxym[0] / 2;for(i = 1; i <= n; i++)fxym[i] = (fxym[i] - a[i]*fxym[i-1]) / (2 - a[i]*B[i-1]);for(i = n-1; i >= 0; i--)fxym[i] = fxym[i] - B[i]*fxym[i+1];}void printout(int n);int main(){int n,i; char ch;do{cout<<"Please put in the number of the dots:";cin>>n;for(i = 0; i <= n; i++){cout<<"Please put in X"<<i<<':';cin>>x[i]; //cout<<endl;cout<<"Please put in Y"<<i<<':';cin>>y[i]; //cout<<endl;}for(i = 0; i < n; i++) //求步长h[i] = x[i+1] - x[i];cout<<"Please 输入边界条件\n 1: 已知两端的一阶导数\n 2:两端的二阶导数已知\n 默认:自然边界条件\n";int t;float f0, f1;cin>>t;switch(t){case 1:cout<<"Please put in Y0\' Y"<<n<<"\'\n";cin>>f0>>f1;c[0] = 1; a[n] = 1;fxym[0] = 6*((y[1] - y[0]) / (x[1] - x[0]) - f0) / h[0];fxym[n] = 6*(f1 - (y[n] - y[n-1]) / (x[n] - x[n-1])) / h[n-1];break;case 2:cout<<"Please put in Y0\" Y"<<n<<"\"\n";cin>>f0>>f1;c[0] = a[n] = 0;fxym[0] = 2*f0; fxym[n] = 2*f1;break;default:cout<<"不可用\n";//待定};//switchfor(i = 1; i < n; i++)fxym[i] = 6 * f(i-1, i, i+1);for(i = 1; i < n; i++){a[i] = h[i-1] / (h[i] + h[i-1]);c[i] = 1 - a[i];}a[n] = h[n-1] / (h[n-1] + h[n]);cal_m(n);cout<<"\n输出三次样条插值函数：\n";printout(n);cout<<"Do you to have anther try ? y/n :";cin>>ch;}while(ch == 'y' || ch == 'Y');return 0;}void printout(int n){cout<<setprecision(6);for(int i = 0; i < n; i++){cout<<i+1<<": ["<<x[i]<<" , "<<x[i+1]<<"]\n"<<"\t";/*cout<<fxym[i]/(6*h[i])<<" * ("<<x[i+1]<<" - x)^3 + "<<<<" * (x - "<<x[i]<<")^3 + "<<(y[i] - fxym[i]*h[i]*h[i]/6)/h[i]<<" * ("<<x[i+1]<<" - x) + "<<(y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i]<<"(x - "<<x[i]<<")\n";cout<<endl;*/float t = fxym[i]/(6*h[i]);if(t > 0)cout<<t<<"*("<<x[i+1]<<" - x)^3";else cout<<-t<<"*(x - "<<x[i+1]<<")^3";t = fxym[i+1]/(6*h[i]);if(t > 0)cout<<" + "<<t<<"*(x - "<<x[i]<<")^3";else cout<<" - "<<-t<<"*(x - "<<x[i]<<")^3";cout<<"\n\t";t = (y[i] - fxym[i]*h[i]*h[i]/6)/h[i];if(t > 0)cout<<"+ "<<t<<"*("<<x[i+1]<<" - x)";else cout<<"- "<<-t<<"*("<<x[i+1]<<" - x)";t = (y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i];if(t > 0)cout<<" + "<<t<<"*(x - "<<x[i]<<")";else cout<<" - "<<-t<<"*(x - "<<x[i]<<")";cout<<endl<<endl;}cout<<endl;}。

一、引言在计算机编程和数据处理领域，插值是一种常见的数值分析方法，用于在已知数据点之间估算未知点的数值。

而三次样条插值是插值方法中的一种重要技术，它可以在使用较少插值节点的情况下，实现更为平滑和精确的插值结果。

本文将着重探讨三次样条插值的原理和C++代码实现，并给出详细的注释和解释。

二、三次样条插值的原理三次样条插值是一种分段插值方法，它将整个插值区间分割为若干个小区间，每个小区间内采用三次多项式进行插值。

这样做的好处是可以在每个小区间内实现更为细致和精确的插值，从而提高插值的准确性和平滑性。

而三次样条插值的核心在于确定每个小区间内的三次多项式的系数，一般采用自然边界条件进行求解。

在具体实现中，我们需要先对给定的插值节点进行排序，并求解出每个小区间内的三次多项式系数。

最终将这些系数整合起来，就可以得到整个插值区间的三次样条插值函数。

三、C++代码实现及注释接下来，我们将给出使用C++语言实现三次样条插值的代码，并对每个关键步骤进行详细注释和解释。

```cpp// include necessary libraries#include <iostream>#include <vector>using namespace std;// define the function for cubic spline interpolationvector<double> cubicSplineInterpolation(vector<double> x, vector<double> y) {// initialize necessary variables and containersint n = x.size();vector<double> h(n-1), alpha(n), l(n), mu(n), z(n), c(n), b(n), d(n);vector<double> interpolatedValues;// step 1: calculate the differences between x valuesfor (int i = 0; i < n-1; i++) {h[i] = x[i+1] - x[i];}// step 2: calculate alpha valuesfor (int i = 1; i < n-1; i++) {alpha[i] = (3/h[i]) * (y[i+1] - y[i]) - (3/h[i-1]) * (y[i] - y[i-1]); }// step 3: calculate l, mu, and z valuesl[0] = 1;mu[0] = 0;z[0] = 0;for (int i = 1; i < n-1; i++) {l[i] = 2*(x[i+1] - x[i-1]) - h[i-1]*mu[i-1];mu[i] = h[i]/l[i];z[i] = (alpha[i] - h[i-1]*z[i-1])/l[i];}l[n-1] = 1;z[n-1] = 0;c[n-1] = 0;// step 4: calculate coefficients for the cubic polynomials for (int j = n-2; j >= 0; j--) {c[j] = z[j] - mu[j]*c[j+1];b[j] = (y[j+1] - y[j])/h[j] - h[j]*(c[j+1] + 2*c[j])/3;d[j] = (c[j+1] - c[j])/(3*h[j]);}// step 5: interpolate values using the cubic polynomials for (int i = 0; i < n-1; i++) {double xi = x[i];while (xi < x[i+1]) {double dx = xi - x[i];double interpolatedValue = y[i] + b[i]*dx + c[i]*dx*dx + d[i]*dx*dx*dx;interpolatedValues.push_back(interpolatedValue);xi += 0.1; // adjust the step size for finer interpolation }}return interpolatedValues;}// main function for testing the cubic spline interpolation int main() {vector<double> x = {1, 2, 3, 4, 5};vector<double> y = {3, 6, 8, 10, 15};vector<double> interpolatedValues = cubicSplineInterpolation(x, y);for (int i = 0; i < interpolatedValues.size(); i++) {cout << "Interpolated value " << i << " : " << interpolatedValues[i] << endl;}return 0;}```四、总结与展望通过本文的学习，我们了解了三次样条插值的原理和C++代码实现。

数值分析三次样条插值函数【问题】对函数f x =ex, x∈[0,1]构造等距节点的三次样条插值函数，对以下两种类型的样条函数1. 三次自然样条2. 满足S′ 0 =1,S′ 1 =e的样条并计算如下误差：max{ f x1 −S x1 ,i=1,…,N} i−i−i这里xi−1为每个小区间的中点。

对N=10,20,40比较以上两组节点的结果。

讨论你的结果。

【三次样条插值】在每一个区间[t1,t2],…,[tn−1,tn]上，S都是不同的三次多项式，我们把在[ti−1,ti]上表示S的多项式记为Si，从而，S0 x x∈[t0,t1]∈[t1,t2] S x = S1 x x…Sn−1 x x∈[tn−1,tn]通过在节点处函数值、一阶导数和二阶导数的连续性可以得到：Si−1 ti = yi= Si ti 1≤i≤ n−1Si−1′ ti = Si′ tix→ti+limS′′ x =zi=limS′′(x) x→ti−再给定z0和zn 的值就构成了4n个条件，而三次样条插值函数共4n个系数，故可以通过这4n个条件求解三次样条函数的系数，从而求得该三次样条插值函数。

特别的，当z0=zn=0 时称为自然三次样条。

文本预览：一、自然三次样条插值【自然三次样条插值算法】1.由上面的分析可知，求解三次样条函数实际上就是求解一个矩阵：u 1h 1h1u2h2h2u3…v1 z1 v2 z2 z3=v3 … z…hn−2 n−2 vn−2 z vn−1 un−1 n−1ih3…hn−3un−2hn−26…其中hi=ti+1−ti，ui=2(hi+hi−1)，ui=h(yi+1−yi)，vi=bi−bi−1 所以自然三层次样条插值的算法就是在得到端点的函数值，一次导数值和二次导数值，然后根据上述求解矩阵得到v，代入自然三次样条的表达式即可。

2.根据题目中所给出的误差估计，计算在区间中点处的最大误差。

【实验】通过Mathematica编写程序得到如下结果：N=101. 计算得到zi的值为：由此可以得到各个区间的自然三次样条插值函数。

实验报告：牛顿差值多项式&三次样条... . (1)问题：在区间［-1,1］上分别取n=10、20用两组等距节点对龙格函数f (x)---作多项式插25 x 2值及三次样条插值对每个n值，分别画出插值函数矽(x)的图形。

实验目的：通过编程实现牛顿插值方法和三次样条方法，加深对多项式插值的理解。

应用所编程序解决实际算例。

实验要求：1.认真分析问题，深刻理解相关理论知识并能熟练应用；2.编写相关程序并进行实验；3.调试程序，得到最终结果；4.分析解释实验结果；5.按照要求完成实验报告。

实验原理：详见《数值分析第5版》第二章相关容。

实验容：(1)牛顿插值多项式1.1 当 n=10 时：在Matlab下编写代码完成计算和画图。

结果如下：代码：clear allclcx1=-1:0.2:1;y1=1./(1+25.*x1.八2);n=length(x1);f=y1(：);for j=2:nfor i=n:-1:jf(i) = (f(i)-f(i-1))/(x1(i)-x1(i-j+1));endendsyms F x p;F(1)=1;p(1)=y1(1);for i=2:nF(i)=F(i-1)*(x-x1(i-1));p(i)=f(i)*F(i);endsyms PP=sum(p);P10=vpa(expand(P),5);x0=-1:0.001:1;y0=subs(P,x,x0);y2=subs(1/(1+25火x八2),x,x0);plot(x0,y0,x0,y2)grid onxlabel('x')ylabel('y')P10即我们所求的牛顿插值多项式，其结果为：P10（x ）=-220.94*x A10+494.91*x A8-9.5065e-14*x A7-381.43*x A6-8.504e-14*x A5+123.36*x A4+2.0202e-14*x A3-16.855*x A2-6.6594e-16*x+1.0并且这里也能得到该牛顿插值多项式的在［-1，1］上的图形，并和原函数进行对比（见Fig.1）。

数值计算方法作业实验4.3 三次样条差值函数实验目的：掌握三次样条插值函数的三弯矩方法。

实验函数:dt ex f xt ⎰∞--=2221)(π实验内容：（1） 编程实现求三次样条插值函数的算法，分别考虑不同的边界条件； （2） 计算各插值节点的弯矩值；（3） 在同一坐标系中绘制函数f(x)，插值多项式，三次样条插值多项式的曲线比较插值结果。

实验4.5 三次样条差值函数的收敛性实验目的：多项式插值不一定是收敛的，即插值的节点多，效果不一定好。

对三次样条插值函数如何呢？理论上证明三次样条插值函数的收敛性是比较困难的，通过本实验可以证明这一理论结果。

实验内容：按照一定的规则分别选择等距或非等距的插值节点，并不断增加插值节点的个数。

实验要求：（1） 随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情况，分析所得结果并与拉格朗日插值多项式比较；（2） 三次样条插值函数的思想最早产生于工业部门。

作为工业应用的例子，考虑如下例子：某汽车制造商根据三次样条插值函数设计车门曲线，其中一段数据如下：kx012345678910 ky0.00.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29ky'0.80.2算法描述：拉格朗日插值：其中是拉格朗日基函数，其表达式为：()∏≠=--=nijj jiji xxxxxl)()(牛顿插值：))...()(](,...,,[....))(](,,[)0](,[)()(11211211----++--+-+=nnnxxxxxxxxxxfxxxxxxxfxxxxfxfxN其中⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎨⎧--=--=--=-)/(]),...,[],...,[(]...,[..],[],[],,[)()(],[11211xxxxxfxxxfxxxfxxxxfxxfxxxfxxxfxfxxfnnnnikjikjkjijijiji三样条插值：所谓三次样条插值多项式Sn(x)是一种分段函数，它在节点Xi(a<X0<X1……<Xn<b)分成的每个小区间[xi-1,xi]上是三次多项式，其在此区间上的表达式如下：],[),6()6(]6)([6)(6)()(111113131iiiiiiiiiiiiiiiiiiiiixxxhyMhMhhyxMMhhyyhxxMihxxMxS-------∈-+-+---+-+-=式中Mi=)(ixS''.因此，只要确定了Mi 的值，就确定了整个表达式，Mi 的计算方法如下：令⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=---+=+=+=+--++++++],,[6)(6111111111i i i i i i i i i i i i i i i i i i i ix x x f h y y h y y h h d h h h h h h λμ则Mi 满足如下n-1个方程：1,...2,1,211-==+++-n i d M M M i i i i i i λμ 常用的边界条件有如下几类：（1） 给定区间两端点的斜率m 0,m n ,即n n n m y x S m y x S ='='='=')(,)(000 （2） 给定区间两端点的二阶导数M0，Mn,即n n n M y x S M y x S =''=''=''='')(,)(000 （3） 假设y=f(x)是以b-a 为周期的周期函数，则要求三次样条插值函数S （x ）也为周期函数，对S （x ）加上周期条件2,1,0),0()0()(0)(=-=+p x S x S n p p对于第一类边界条件有⎪⎪⎩⎪⎪⎨⎧--=+--=+--)(62)(6211001110n n n n n n i h y y mn h M M m h y y h M M对于第二类边界条件有⎩⎨⎧=+=+-n n n n d M M d M M 221100μλ其中n n n n nnn M u x x f m h d M m x x f h d )1(2]),[(6)1(2)],[(6100001010-+-=-+-=-μλλ那么解就可以为⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----n n n n n n n d d d d d M M M M M 1210121011...2...............2............................1..2.1......0..2μλμλμλ对于第三类边界条件，)0()0(,,000-=+==n n n x S x S M M y y ，由此推得0010012d M M M n =-++μλ，其中]),1[],[(6,,101010110n n nn n n x x f x x f h h d h h h h h h --+=+=+=μλ，那么解就可以为： ⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛-------1221012101221100...2.............2..............................2..,,.......,..22n n n n n n n d d d d d M M M M M n μλλμλμμλ 程序代码： 1拉格朗日插值函数Lang.mfunction f=lang(X,Y,xi) %X 为已知数据的横坐标 %Y 为已知数据的纵坐标 %xi 插值点处的横坐标%f 求得的拉格朗日插值多项式的值 n=length(X); f=0; for i=1:n l=1; for j=1:i-1l=l.*(xi-X(j))/(X(i)-X(j)); end ; for j=i+1:nl=l.*(xi-X(j))/(X(i)-X(j)); end ;%拉格朗日基函数 f=f+l*Y(i); endfprintf('%d\n',f) return2 牛顿插值函数newton.mfunction f=newton(X,Y,xi) %X 为已知数据的横坐标 %Y 为已知数据的纵坐标%xi插值点处的横坐标%f求得的拉格朗日插值多项式的值n=length(X);newt=[X',Y'];%计算差商表for j=2:nfor i=n:-1:1if i>=jY(i)=(Y(i)-Y(i-1))/(X(i)-X(i-j+1));else Y(i)=0;endendnewt=[newt,Y'];end%计算牛顿插值f=newt(1,2);for i=2:nz=1;for k=1:i-1z=(xi-X(k))*z;endf=f+newt(i-1,i)*z;endfprintf('%d\n',f)return3三次样条插值第一类边界条件Threch.mfunction S=Threch1(X,Y,dy0,dyn,xi)% X为已知数据的横坐标%Y为已知数据的纵坐标%xi插值点处的横坐标%S求得的三次样条插值函数的值%dy0左端点处的一阶导数% dyn右端点处的一阶导数n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:n%求函数的一阶差商h(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i))/h(i);endfor i=2:n%求函数的二阶差商f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i)=6*f2(i);endd(1)=6*(f1(1)-dy0)/h(1);d(n+1)=6*(dyn-f1(n-1))/h(n-1);%¸赋初值A=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1));C(i)=1-B(i);endA(1,2)=1;A(n+1,n)=1;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=A\d;syms x;for i=1:nSx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))...+M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3);digits(4);Sx(i)=vpa(Sx(i));%三样条插值函数表达式endfor i=1:ndisp('S(x)=');fprintf('%s (%d,%d)\n',char(Sx(i)),X(i),X(i+1));endfor i=1:nif xi>=X(i)&&xi<=X(i+1)S=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2+(M(i+1)-M(i))/(6 *h(i))*(xi-X(i))^3;endenddisp('xi S');fprintf('%d,%d\n',xi,S);return4 三次样条插值第二类边界条件Threch2.mfunction [Sx]=Threch2(X,Y,d2y0,d2yn,xi)X为已知数据的横坐标%Y为已知数据的纵坐标%xi插值点处的横坐标%S求得的三次样条插值函数的值%d2y0左端点处的二阶导数% d2yn右端点处的二阶导数n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:n%求一阶差商h(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i))/h(i);endfor i=2:n%求二阶差商f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i)=6*f2(i);endd(1)=2*d2y0;d(n+1)=2*d2yn;%赋初值A=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1));C(i)=1-B(i);endA(1,2)=0;A(n+1,n)=0;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=A\d;syms x;for i=1:nSx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))... +M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3);digits(4);Sx(i)=vpa(Sx(i));endfor i=1:ndisp('S(x)=');fprintf('%s (%d,%d)\n',char(Sx(i)),X(i),X(i+1));endfor i=1:nif xi>=X(i)&&xi<=X(i+1)S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2+(M(i+1)-M(i)) /(6*h(i))*(xi-X(i))^3;endenddisp('xi S');fprintf('%d,%d\n',xi,S);return5插值节点处的插值结果main3.mclearclcX=[0.0,0.1,0.2,0.3,0.4];Y=[0.5000,0.5398,0.5793,0.6179,0.7554];xi=0.13;%xi=0.36;disp('xi=0.13');%disp('xi=0.36');disp('拉格朗日插值结果');lang(X,Y,xi);disp('牛顿插值结果');newton(X,Y,xi);disp('三次样条第一类边界条件插值结果');Threch1(X,Y,0.40,0.36,xi);%0.4,0.36分别为两端点处的一阶导数disp('三次样条第二类边界条件插值结果');Threch2(X,Y,0,-0.136,xi);%0,-0.136分别为两端点处的二阶导数6将多种插值函数即原函数图像画在同一张图上main2.mclearclcX=[0.0,0.1,0.2,0.3,0.4];Y=[0.5000,0.5398,0.5793,0.6179,0.7554];a=linspace(0,0.4,21);NUM=21;L=zeros(1,NUM);N=zeros(1,NUM);S=zeros(1,NUM);B=zeros(1,NUM);for i=1:NUMxi=a(i);L(i)=lang(X,Y,xi);% 拉格朗日插值N(i)=newton(X,Y,xi);%牛顿插值B(i)=normcdf(xi,0,1);%原函数S(i)=Threch1(X,Y,0.4,0.36,xi);%三次样条函数第一类边界条件endplot(a,B,'--r');hold on;plot(a,L,'b');hold on;plot(a,N,'r');hold on;plot(a,S,'r+');hold on;legend('原函数','拉格朗日插值','牛顿插值','三次样条插值',2);hold off7增加插值节点观察误差变化main4.mclear;clc;N=5;%4.5第一问Ini=zeros(1,1001);a=linspace(-1,1,1001);Ini=1./(1+25*a.^2);for i=1:3 %节点数量变化次数N=2*N;t=linspace(-1,1,N+1);%插值节点ft=1./(1+25*t.^2);%插值节点函数值val=linspace(-1,1,101);for j=1:101L(j)=lang(t,ft,val(j));S(j)=Threch1(t,ft,0.074,-0.074,val(j));%三样条第一类边界条件插值endplot(a,Ini,'k')%原函数图象hold onplot(val,L,'r')%拉格朗日插值函数图像hold onplot(val,S,'b')%三次样条插值函数图像str=sprintf('插值节点为%d时的插值效果',N);title(str);legend('原函数','拉格朗日插值','三次样条插值');%显示图例hold offfigureend8车门曲线main5.mclearclcX=[0,1,2,3,4,5,6,7,8,9,10];Y=[0.0,0.79,1.53,2.19,2.71,3.03,3.27,2.89,3.06,3.19,3.29]; dy0=0.8;dyn=0.2;n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:nh(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i))/h(i);endfor i=2:nf2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i)=6*f2(i);endd(1)=6*(f1(1)-dy0)/h(1);d(n+1)=6*(dyn-f1(n-1))/h(n-1); A=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1));C(i)=1-B(i);endA(1,2)=1;A(n+1,n)=1;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=A\d;x=zeros(1,n);S=zeros(1,n);for i=1:nx(i)=X(i)+0.5;S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x(i)-X(i))+M(i)/2*(x(i)-X(i))^2+(M(i+1)-M (i))/(6*h(i))*(x(i)-X(i))^3;endplot(X,Y,'k'); hold on;plot(x,S,'o');title('三次样条插值效果图');legend('已知插值节点','三次样条插值');hold off实验结果：4.31计算插值节点处的函数值xi=0.13时Xi=0.36时2将多种插值函数即原函数图像画在同一张图上4.5.1增加插值节点观察误差变化从上面三张图可以看出增加插值节点并不能改善差之效果4.5.2 车门曲线（注：专业文档是经验性极强的领域，无法思考和涵盖全面，素材和资料部分来自网络，供参考。