三次样条插值法《数值分析》上机实验作业

昆明理工大学研究生《数值分析》上机实验作业姓名：学号：专业：一、课题名称课题七三次样条插值法1、问题提出设已知数据如下：0.20.40.60.81.00.97986520.91777100.80803480.63860920.3843735求f(x)的三次样条插值函数S(x)。2、要求（1）满足自然边界条件0)0.1()2.0(ss；（2）满足第一类边界条件20271.0)2.0(s，55741.1)0.1(s。（3）打印输出用追赶法解出的弯曲向量(0M,1M,2M,3M,4M)和)1.02.0(iS（i=0,1,2,3,4,5,6,7,8）的值。并画出)(xSy的图形。二、班级、姓名、学号三、目的和意义由于航空、造船等工程设计的需要而发展起来所谓样条插值方法，既保留了分段低次插值多项式的各种优点，又提高了插值函数的光滑性，而且具有较好的稳定性。今天，样条插值方法已成为数值逼近的一个极其重要的分支，在许多领域里得到越来越广泛地应用。其中，尤以三次样条插值函数应用最为广泛，如在高速飞机的机翼形体和船体放样等方面的应用，同时在计算机作图方面更是大有作为。它能够解决一些既有二阶光滑度，又有二阶连续导数的方程，具有良好的收敛性和稳定性。1.通过本次实验进一步了解三次样条插值函数，并通过求解三弯矩方程组得出曲线函数组；2.通过MATLAB编程实现求三次样条插值函数的算法，分别考虑不同的边界条件，同时用追赶法解出弯曲向量和)1.02.0(iS（i=0,1,2,3,4,5,6,7,8）的值。四、计算公式首先我们利用)(xS的二阶导数值),,2,1,0()(jnjMxSj表达)(xS,因为在区间]x,[x1jj上)()(jxSxS是不高于三次的多项式，其二阶导数)(xS必是线性函数，所以可表示为：]x,[xx,hxxMhxxM(x)S1jjjj1jj1jj对)(xS积分两次并利用1j1jjjy)(xS,y)(xS，可定出积分常数，于是得三次样条表达式。jj2j1j1jj1j2jjjj3j1jj31jjhxx)h6M(yhxx)h6M(y6h)x(xM6h)x(xM(x)S1n,0,1,2,，j这里),,1,0(njMj是未知的。为了确定),,1,0(njMj，对)(xS求导得)M(M6hhyy2h)x(xM2hx)(xM(x)Sj1jjjj1jj2j1jj21jj由此可得jj1jjjj1jjM3hM6hhyy0)(xS类似的可求出)(xS在区间]x,[xj1-j上的表达式，进而得j1j1j1j1j1jjjM3hM6hhyy0)(xS利用0)(xS0)(xSjj可得j1jjj1jjdMλ2MMμ，1,n,1,j（三弯矩方程）其中j1j1jjhhhμ，j1jjjhhhλ，1jj1jj1jj1j1jjjx,x,x6fhhx,xfx,xf6d.1,n,1,j其中有（1n）个未知数nMMM,...,10，而方程只有（n-1）个，当满足第一种边界条件时，可得另两个方程010010,62fxxfhMM，nnnnnnxxffhMM,62111如果令nnnnnnxxffhgfxxfhd,6,1,,6,111010000，将上述方程综合后的一下矩阵形式：n1n10n1n10n1n1n110ddddMMMM2μλ2μλ2μλ2可以证明此方程组满足追赶法的条件，我们用追赶法可得M的)1(n值，将其带入公式即得)(xs。对第二种边界条件，直接的端点方程nnfMfM,00并且令nnnfdfd2,2,0000，则又得三弯矩方程同理即可求得解。五、结构程序设计1.满足自然边界条件时自定义函数：followup.m%追赶法求m%A为线性方程组的系数矩阵%b为常数向量functionm=followup(A,b)n=rank(A);fori=1:nifA(i,i)==0disp('error:对角元素中有数据为0');return;endendd=ones(n,1);a=ones(n-1,1);c=ones(n-1);fori=1:n-1a(i,1)=A(i+1,i);c(i,1)=A(i,i+1);d(i,1)=A(i,i);endd(n,1)=A(n,n);fori=2:nd(i,1)=d(i,1)-(a(i-1,1)/d(i-1,1))*c(i-1,1);b(i,1)=b(i,1)-(a(i-1,1)/d(i-1,1))*b(i-1,1);endm(n,1)=b(n,1)/d(n,1);fori=(n-1):-1:1m(i,1)=(b(i,1)-c(i,1)*m(i+1,1))/d(i,1);end自定义函数：thrsample2.m%a为要求的插值点%f为区间内的插值函数%f0为输入点处的插值%m为追赶法解出的弯矩向量functionthrsample2(a)x=[0.2:0.2:1.0];y=[0.97986520.91777100.80803480.63860920.3843735];s02=0;s10=0;x0=a;n=length(x);fori=1:nif(x(i)=x0)&&(x(i+1)=x0)index=i;break;endendA=diag(2*ones(1,n));A(1,2)=1;A(n,n-1)=1;u=zeros(n-2,1);lamda=zeros(n-1,1);c=zeros(n,1);fori=2:n-1u(i-1)=(x(i)-x(i-1))/(x(i+1)-x(i-1));lamda(i)=(x(i+1)-x(i))/(x(i+1)-x(i-1));c(i)=3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i));A(i,i+1)=u(i-1);A(i,i-1)=lamda(i);endc(1)=3*(y(2)-y(1))/(x(2)-x(1))-(x(2)-x(1))*s02/2;c(n)=3*(y(n)-y(n-1))/(x(n)-x(n-1))-(x(n)-x(n-1))*s10/2;m=followup(A,c)h=x(index+1)-x(index);symst;f=y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h+m(index)*(t-x(index))*(x(index+1)-t)^2/h/h-m(index+1)*(x(index+1)-t)*(t-x(index))^2/h/hf0=subs(f,'t',x0)运行结果（方程S、弯矩M和插值函数f的值）为：S0=0.9798652000S1=0.9525726703S2=0.9177710000S3=0.8695049390S4=0.8080348000S5=0.7334085485S6=0.6386092000S7=0.5187304296S8=0.3843735000f=125*((8825841099186633*t)/4503599627370496-8825841099186633/45035996273704960)*(t-2/5)^2-125*((8266546267222895*t)/4503599627370496-8266546267222895/9007199254740992)*(t-1/5)^2-25*((7396583675403915*t)/18014398509481984-1479316735080783/9007199254740992)*(t-1/5)^2-25*(t-2/5)^2*((2290591133669*t)/8796093022208-2290591133669/43980465111040)2.满足第一类边界条件时自定义函数：thrsample1.m%a为要求的插值点%f为区间内的插值函数%f0为输入点处的插值%m为追赶法解出的弯矩向量functionthrsample1(a)x=[0.2:0.2:1.0];y=[0.97986520.91777100.80803480.63860920.3843735];s02=0.20271;s10=1.55741;x0=a;n=length(x);fori=1:nif(x(i)=x0)&&(x(i+1)=x0)index=i;break;endendA=diag(2*ones(1,n));u=zeros(n-2,1);lamda=zeros(n-1,1);c=zeros(n,1);fori=2:n-1u(i-1)=(x(i)-x(i-1))/(x(i+1)-x(i-1));lamda(i)=(x(i+1)-x(i))/(x(i+1)-x(i-1));c(i)=3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i));A(i,i+1)=u(i-1);A(i,i-1)=lamda(i);endc(1)=2*s02;c(n)=2*s10;m=followup(A,c)h=x(index+1)-x(index);symst;f=y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h+m(index)*(t-x(index))*(x(index+1)-t)^2/h/h-m(index+1)*(x(index+1)-t)*(t-x(index))^2/h/hf0=subs(f,'t',x0)运行结果（方程S、弯矩M和插值函数f的值）为：S0=0.9798652000S1=0.9685577748S2=0.9177710000S3=0.8590474261S4=0.8080348000S5=0.7592534958S6=0.6386092000S7=0.4258081535S8=0.3843735000f=125*((8825841099186633*t)/4503599627370496-8825841099186633/45035996273704960)*(t-2/5)^2-125*((8266546267222895*t)/4503599627370496-8266546267222895/9007199254740992)*(t-1/5)^2+25*((20271*t)/100000-20271/500000)*(t-2/5)^2-25*((1321529499150801*t)/2251799813685248-1321529499150801/5629499534213120)*(t-1/5)^23.画出y=S(x)的图形（1）满足自然边界条件时程序为：x=0.2:0.2:1;y=[0.97986520.91777100.80803480.63860930.3843735];xi=0.2:0.01:1.0;yi=interp1(x,y,xi,'variational');plot(x,y,'o',xi,yi,'k-');输出图形为：图1自然边界条件下的y=S(x)的图形（2