三次样条插值作业题

1、例1 1设f(x)为定义在0,30,3上的函数，有下列函数值表:X Xi i0 01 12 23 3y yi i0 00.50.52 21.51.5且f(x0)=0.2,f(X3)=-1,试求区间0,30,3上满足上述条件的三次样条插值函数s(x)本算法求解出的三次样条插值函数将写成三弯矩方程的形式:Mjhj2Mj12MiM(yj-a)(yj1-).，, ,6 6, ,6 6将由MatlabMatlab6hj6hjhjhjCoefs_2Coefs_3Coefs_2Coefs_3以及Coefs_4Coefs_4的值求出以下为MatlabMatlab代码：%=%本段代码解决作业题的例i%=clea

2、rallclc%自变量x与因变量y,两个边界条件的取值IndVar=0,1,2,3;DepVar=0,0.5,2,1.5;LeftBoun=0.2;RightBoun=-1;%区间长度向量，其各元素为自变量各段的长度h=zeros(1,length(IndVar)-1);fori=1:length(IndVar)-1h(i)=IndVar(i+1)-IndVar(i);end%为向量科赋值s(x)=Mj6hj3Mj13(jT)所(x)M白2(yj6hjMj后)(yji-)-(xj1-x)-6(x-xj)hj其中，方程中的系数代码中的变量Coefs_1Coefs_1、mu=zeros(1,len

3、gth(h);fori=1:length(mu)-1mu(i)=h(i)/(h(i)+h(i+1);endmu(i+1)=1;%为向量入赋值lambda=zeros(1,length(h);lambda(1)=1;fori=2:length(lambda)lambda(i)=h(i)/(h(i-1)+h(i);end%为向量d赋值d=zeros(1,length(h)+1);d(1)=6*(DepVar(2)-DepVar(1)/(IndVar(2)-IndVar(1)-LeftBoun)/h(1);fori=2:length(h)a=(DepVar(i)-DepVar(i-1)/(IndVa

4、r(i)-IndVar(i-1);b=(DepVar(i+1)-DepVar(i)/(IndVar(i+1)-IndVar(i);c=(b-a)/(IndVar(i+1)-IndVar(i-1);d(i)=6*c;endd(i+1)=6*(RightBoun-(DepVar(i+1)-DepVar(i)/(IndVar(i+1)-IndVar(i)%为矩阵A赋值%将主对角线上的元素全部置为2A=zeros(length(d),length(d);fori=1:length(d)A(i,i)=2;end%将向量入的各元素赋给主对角线右侧第一条对角线fori=1:length(d)-1A(i,i+

5、1)=lambda(i);end%将向量d的各元素赋给主对角线左侧第一条对角线fori=1:length(d)-1A(i+1,i)=mu(i);end%求解向量MM=Ad;%求解每一段曲线的函数表达式fori=1:length(h)Coefs_1=M(i)/(6*h(i);Part_1=conv(Coefs_1,.conv(-1,IndVar(i+1),.)/h(i);conv(-1,IndVar(i+1),-1,IndVar(i+1);S_1=polyval(Part_1,IndVar(i):0.01:IndVar(i+1);Coefs_2=M(i+1)/(6*h(i);Part_2=con

6、v(Coefs_2,.conv(1,-IndVar(i),.conv(1,-IndVar(i),1,-IndVar(i);S_2=polyval(Part_2,IndVar(i):0.01:IndVar(i+1);Coefs_3=(DepVar(i)-M(i)*h(i)A2/6)/h(i);Part_3=conv(Coefs_3,-1,IndVar(i+1);S_3=polyval(Part_3,IndVar(i):0.01:IndVar(i+1);Coefs_4=(DepVar(i+1)-M(i+1)*h(i42/6)/h(i);Part_4=conv(Coefs_4,1,-IndVar(i

7、);S_4=polyval(Part_4,IndVar(i):0.01:IndVar(i+1);S=S_1+S_2+S_3+S_4;plot(IndVar(i):0.01:IndVar(i+1),S,LineWidth,1.25)%在样条插值曲线的相应位置标注该段曲线的函数表达式text(i-1,polyval(Part_1,3),.itS,num2str(i),(x)=,num2str(Coefs_1),(,num2str(IndVar(i+1),-x)A3+,.num2str(Coefs_2),(x-,num2str(IndVar(i),F3+,num2str(Coefs_3),.(,nu

8、m2str(IndVar(i+1),-x)+,num2str(Coefs_4),(x-,num2str(IndVar(i),),.FontName,TimesNewRoman,FontSize,14)holdonend%过x=1和x=2两个横轴点作垂线%line(1,1,2.5,-0.5,LineStyle,-);line(2,2,2.5,-0.5,LineStyle,-);%为x轴和y轴添加标注xlabel(itx,FontName,TimesNewRoman,.FontSize,14,FontWeight,bold);ylabel(its(x),FontName,TimesNewRoman

9、,.Rotation,0,FontSize,14,FontWeight,bold);最终，三次样条插值函数s(x)|ls(x)|l达式为：0.06(1X)3十0.42x3+0.06(1x)+0.08x,xwb,l,s(x)=0.42(2x)3-0.62(x-1)3+0.08(2x)+2.62(x1),xw1,2】，0.62(3x)3+0.06(x-2)3+2.62(3-x)+1.44(x-2),xw12,31曲线的图像如图所示：例2 2已知函数值表:X Xi i1 12 24 45 5y yi i1 13 34 42 2试求在区间1,51,5上满足上述函数表所给出的插值条件的三次自然样条插值函

10、数s(x)本算法求解出的三次样条插值函数将写成三弯矩方程的形式:Mjh2Mj1h2MiM(yj-a)(yj1-).，, ,6 6, ,6 6将由MatlabMatlab6hj6hjhjhjCoefs_2Coefs_3Coefs_2Coefs_3以及Coefs_4Coefs_4的值求出以下为MatlabMatlab代码：%=%本段代码解决作业题的例2%=clearallclc%自变量X与因变量y的取值IndVar=1,2,4,5;DepVar=1,3,4,2;%区间长度向量，其各元素为自变量各段的长度h=zeros(1,length(IndVar)-1);fori=1:length(IndVar

11、)-1h(i)=IndVar(i+1)-IndVar(i);end%为向量科赋值mu=zeros(1,length(h);fori=1:length(mu)-1s(x)=Mj6hj3Mj13(jT)所(X)M讨(yj6hjMj后)(山1一)-(Xj1-X)-6(X-Xj)hj其中，方程中的系数代码中的变量Coefs_1Coefs_1、mu(i)=h(i)/(h(i)+h(i+1);endmu(i+1)=0;%为向量入赋值lambda=zeros(1,length(h);lambda(1)=0;fori=2:length(lambda)lambda(i)=h(i)/(h(i-1)+h(i);en

12、d%为向量d赋值d=zeros(1,length(h)+1);d(1)=0;fori=2:length(h)a=(DepVar(i)-DepVar(i-1)/(IndVar(i)-IndVar(i-1);b=(DepVar(i+1)-DepVar(i)/(IndVar(i+1)-IndVar(i);c=(b-a)/(IndVar(i+1)-IndVar(i-1);d(i)=6*c;endd(i+1)=0;%为矩阵A赋值%将主对角线上的元素全部置为2A=zeros(length(d),length(d);fori=1:length(d)A(i,i)=2;end%将向量入的各元素赋给主对角线右侧第

13、一条对角线fori=1:length(d)-1A(i,i+1)=lambda(i);end%将向量d的各元素赋给主对角线左侧第一条对角线fori=1:length(d)-1A(i+1,i)=mu(i);end%求解向量MM=Ad;%求解每一段曲线的函数表达式fori=1:length(h)Coefs_1=M(i)/(6*h(i);Part_1=conv(Coefs_1,.conv(-1,IndVar(i+1),.conv(-1,IndVar(i+1),-1,IndVar(i+1);S_1=polyval(Part_1,IndVar(i):0.01:IndVar(i+1);Coefs_2=M(i

14、+1)/(6*h(i);Part_2=conv(Coefs_2,.conv(1,-IndVar(i),.conv(1,-IndVar(i),1,-IndVar(i);S_2=polyval(Part_2,IndVar(i):0.01:IndVar(i+1);Coefs_3=(DepVar(i)-M(i)*h(i)A2/6)/h(i);Part_3=conv(Coefs_3,-1,IndVar(i+1);S_3=polyval(Part_3,IndVar(i):0.01:IndVar(i+1);Coefs_4=(DepVar(i+1)-M(i+1)*h(i42/6)/h(i);Part_4=co

15、nv(Coefs_4,1,-IndVar(i);S_4=polyval(Part_4,IndVar(i):0.01:IndVar(i+1);S=S_1+S_2+S_3+S_4;plot(IndVar(i):0.01:IndVar(i+1),S,LineWidth,1.25)%在样条插值曲线的相应位置标注该段曲线的函数表达式text(i,polyval(Part_1,5),.itS,num2str(i),(x)=,num2str(Coefs_1),(,num2str(IndVar(i+1),-x)A3+,.num2str(Coefs_2),(x-,num2str(IndVar(i),F3+,nu

16、m2str(Coefs_3),.(,num2str(IndVar(i+1),-x)+,num2str(Coefs_4),(x-,num2str(IndVar(i),),.FontName,TimesNewRoman,FontSize,14)holdonend%过x=2和x=4两个横轴点作垂线%line(2,2,4.5,0.5,LineStyle,-);line(4,4,4.5,0.5,LineStyle,-);%为x轴和y轴添加标注xlabel(itx,FontName,TimesNewRoman,.FontSize,14,FontWeight,bold);ylabel(its(x),Font

17、Name,TimesNewRoman,.Rotation,0,FontSize,14,FontWeight,bold);最终，三次自然样条插值函数s(xs(x层达式为：-0.125(x-1)3+(2-x)+3.125(x-1),xw1,2】，s(x)=0.0625(4x)3-0.1875(x-2)3+1.75(4-x)+2.75(x-2),xw12,4,0.375(5x)3+4.375(5x)+2(x4),xw4,51曲线的图像如图所示：例3 3课后习题与思考题第7 7题X Xi i0.250.250.300.300.390.390.450.450.530.53y yi i0.50000.50

18、000.54770.54770.62450.62450.67080.67080.72800.7280试求在区间0.25,0.530.25,0.53比满足上述函数表所给出的插值条件的三次自然样条插值函数s(x)s(x)求解出的三次样条插值函数将写成三弯矩方程的形式:MjhjMi3Mli3-一)Ns(x)(Xj1-X)(x-Xj)6一(Xj.1-X)6%6%)-(x-xj)本题采用和例2 2基本相同的MatlabMatlab代码，只改变初始条件。最终，三次自然样条插值函数s(xs(x层达式为:s(x)=-6,2652(x-0,25)3+10(0.3-x)+10.9697(x-0.25),xb,25

19、,0.30，33-3.4806(0.39-x)-1.5993(x-0.3)+6,1137(0.39-x)+6.9518(x-0.3),33-2.399(0.45-x)-2.859(x-0.39)+10.417(0.45x)+11.1903(x0.39),-2.1442(0.53-x)3十8.3987(0.53x)+9,1(x0.45),x=0.45,0.53】.0.30,0,391,0.39,0.45】,曲线的图像如图所示:例4 4课后习题与思考题第6 6题求y=的二次插值式P2(x),使:F2(100)=10,P2(121)=11,P2(144)=12(144);并计算Y Y市的近似值并估计

20、误差。本题采用拉格朗日二次插值法进行计算:以下为MatlabMatlab代码：%=%本段代码解决课本第2章习题与思考题第6题%=clearallclc%自变量x与因变量y的取值IndVar=100,121,144;DepVar=10,11,12;%构造拉格朗日插值函数Coefs_1=DepVar(1)/.(IndVar(1)-IndVar(2)*(IndVar(1)-IndVar(3);Part_1=conv(Coefs_1,.conv(1,-IndVar(2),1,-IndVar(3);f_1=polyval(Part_1,IndVar(1):0.01:IndVar(3);Coefs_2=D

21、epVar(2)/.(IndVar(2)-IndVar(1)*(IndVar(2)-IndVar(3);Part_2=conv(Coefs_2,.conv(1,-IndVar(1),1,-IndVar(3);f_2=polyval(Part_2,IndVar(1):0.01:IndVar(3);Coefs_3=DepVar(3)/.(IndVar(3)-IndVar(1)*(IndVar(3)-IndVar(2);Part_3=conv(Coefs_3,.conv(1,-IndVar(1),1,-IndVar(2);f_3=polyval(Part_3,IndVar(1):0.01:IndVar(3);f=f_1+f_2+f_3;plot(IndVar(1):0.01:IndVar(3),f,LineWidth,1.25);%在样条插值曲线的相应位置标注该段曲线的函数表达式text(110,polyval(Part_1,110)+polyval(Part_2,110)+polyval(Part_3,110),.itP_2(x)=