东南大学_数值分析_第四章_三次样条插值.docVIP

第四章 多项式插值与函数最佳逼近 ——曲线拟合之3次样条插值 *****（学号） *****（姓名） 上机题目要求见教材P195,37题。 一、算法原理 题目要求编写第一边界条件的3次样条插值函数的通用程序，同时根据汽车门曲线值点构造三次紧压样条曲线函数。其基本原理如下 定义 设有N+1个点，其中。如果存在N个三次多项式，系数为满足如下性质： 则成为三次样条函数。 现证明其存在： 由于是分段三次多项式，其二阶导数是在区间内是分段线性的。根据线性拉格朗日插值可以表示为： 用代入上式，得 将上式积分两次，会引入两个积分常数，可得到如下形式： 将代入上式，并利用可得两个方程： 求解，并将所得的结果带入方程得 求式（7）的导数，并化简得 由上述方程可得如下方程 重组上述方程，得三角线性方程组，表示为 该式具有严格对角优势。算出系数后可由如下公式计算的样条系数。 二、流程图 用已知的N+1个点构造三次紧压样条曲线的问题。其通用程序流程图1所示。 图1关于求解三阶样条曲线算法流程图 具体步骤如下： 计算，， 联系端点约束条件求解样条函数的系数。 三、计算代码 核心代码 for k=2:N-1 temp=A(k-1)/B(k-1); B(k)=B(k)-temp*C(k-1); U(k)=U(k)-temp*U(k-1); end for k=N-2:-1:1 M(k+1)=(U(k)-C(k)*M(k+2))/B(k); for k=0:N-1 S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1)); S(k+1,2)=M(k+1)/2; S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6; S(k+1,4)=Y(k+1); end 完整代码 function S=Three_fit(X,Y,dx0,dxn) %Input - X is a vector that contains a list of abscissas % - Y is a vector that contains a list of ordinates % - dx0=S(x0) first derivative boundary condition % - dxn=S(xn) first derivative boundary condition %Output - S rows of S are the coefficients , in descending order ,for % the cubic interpolants N=length(X)-1; H=diff(X); D=diff(Y)./H; A=H(2:N-1); B=2*(H(1:N-1)+H(2:N)); C=H(2:N); U=6*diff(D); %clamped spline endpoint constraints B(1)=B(1)-H(1)/2; U(1)=U(1)-3*(D(1)-dx0); B(N-1)=B(N-1)-H(N)/2; U(N-1)=U(N-1)-3*(dxn-D(N)); %Gauss elimination to solve mk for k=2:N-1 temp=A(k-1)/B(k-1); B(k)=B(k)-temp*C(k-1); U(k)=U(k)-temp*U(k-1); end M(N)=U(N-1)/B(N-1); for k=N-2:-1:1 M(k+1)=(U(k)-C(k)*M(k+2))/B(k); end M(1)=3*(D(1)-dx0)/H(1)-M(2)/2; M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2; for k=0:N-1 S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1)); S(k+1,2)=M(k+1)/2; S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6; S(k+1,4)=Y(k+1); end 四、计算结果及分析 根据题目数据调用函数如下 X=0:10; Y=[2.51,3.30,4.04,4.70,5.22,5.54,5.78,5.40,5.57,5.70,5.80]; dx0=0.8;dxn=0.2; S=Three_fit(X,Y,dx0,dxn)