三次样条插值作业题

1、例 1设 f ( x)为定义在 0,3 上的函数，有下列函数值表：x i0123yi0525.01且 f ( x0 ) 0.2， f (x3)1 ，试求区间 0,3上满足上述条件的三次样条插值函数 s( x)本算法求解出的三次样条插值函数将写成三弯矩方程的形式：M jM j 1( y jM j hj2 )( yj 1M j 1h2j )( xj 1x)3( x x j )36( x j 1x)6( x xj )s( x)6hjhjhj6h jM jM j 1( y jM j h2j )( yj 1M j 1hj2)66其中，方程中的系数，将由6hj6hjhjhjMatlab 代码中的变量 Co

2、efs_1、Coefs_2 、Coefs_3以及 Coefs_4的值求出。以下为 Matlab代码：%=% 本段代码解决作业题的例1%=clearallclc% 自变量 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);for i = 1 : length(IndVar) - 1h(i) = IndVar(i + 1) - IndVar(i);end%

3、为向量赋值mu = zeros(1, length(h);for i = 1 : length(mu) - 1mu(i) = h(i) / (h(i) + h(i + 1);endmu(i + 1) = 1;% 为向量赋值lambda = zeros(1, length(h);lambda(1) = 1;for i = 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) ) / ( In

4、dVar(2) - IndVar(1) ) - LeftBoun) / h(1);for i = 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) = 6 *(RightBoun - ( DepVar(i + 1)

5、 - DepVar(i) ) / ( IndVar(i + 1) -IndVar(i) ) / h(i);% 为矩阵 A 赋值% 将主对角线上的元素全部置为 2 A = zeros( length(d), length(d) );for i = 1 : length(d) A(i, i) = 2;end% 将向量的各元素赋给主对角线右侧第一条对角线for i = 1 : length(d) - 1A(i, i + 1) = lambda(i);end% 将向量 d 的各元素赋给主对角线左侧第一条对角线for i = 1 : length(d) - 1A(i + 1, i) = mu(i);en

6、d% 求解向量 M M =A d;% 求解每一段曲线的函数表达式for i = 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 + 1)/(6 * h(i);Part_2 = conv(Coefs_2,.conv(1, -IndVar(i),.

7、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)2 / 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(i)2 / 6) / h(i);Part_4

8、 = conv(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 - 1, polyval(Part_1, 3),.itS, num2str(i),(x)= , num2str(Coefs_1),( , num2str( IndVar(i + 1) ),-x

9、)3+, .num2str(Coefs_2),(x- , num2str( IndVar(i) ),)3+ , num2str(Coefs_3),.(, num2str( IndVar(i + 1) ),-x)+ , num2str(Coefs_4),(x- , num2str( IndVar(i) ),), .FontName, Times New Roman, FontSize, 14)holdonend% 过 x=1 和 x=2 两个横轴点作垂线%line(1, 1, 2.5, -0.5,LineStyle, - );line(2, 2, 2.5, -0.5,LineStyle, - )

10、;% 为 x 轴和 y 轴添加标注xlabel(itx, FontName, Times New Roman, .FontSize, 14,FontWeight, bold );ylabel(its(x), FontName, Times New Roman, .Rotation, 0, FontSize, 14,FontWeight, bold );最终，三次样条插值函数s(x)表达式为：0.06(1x)30.42x30.06(1 x)0.08x , x0,1 ,s( x)0.42(2x) 30.62(x1)30.08(2x)2.62( x1) , x1,2 ,0.62(3x)30.06(x

11、2) 32.62(3x)1.44( x2) ,x 2 ,3 .曲线的图像如图所示：例 2 已知函数值表：x i1245y i1342试求在区间 1,5 上满足上述函数表所给出的插值条件的三次自然样条插值函数 s(x)本算法求解出的三次样条插值函数将写成三弯矩方程的形式：M jM j 1( y jM j hj2 )( yj 1M j 1h2j )( xj 1x)3( x x j )36( x j 1x)6( x xj )s( x)6hjhjhj6h jM jM j 1( y jM j h2j )( yj 1M j 1hj2)66其中，方程中的系数，将由6hj6hjhjhjMatlab 代码中的变

12、量 Coefs_1、Coefs_2 、Coefs_3以及 Coefs_4的值求出。以下为 Matlab代码：%=% 本段代码解决作业题的例2%=clearallclc% 自变量 x 与因变量 y 的取值IndVar = 1, 2, 4, 5;DepVar = 1, 3, 4, 2;% 区间长度向量，其各元素为自变量各段的长度h = zeros(1, length(IndVar) - 1);for i = 1 : length(IndVar) - 1h(i) = IndVar(i + 1) - IndVar(i);end% 为向量赋值mu = zeros(1, length(h);for i =

13、 1 : length(mu) - 1mu(i) = h(i) / (h(i) + h(i + 1);endmu(i + 1) = 0;% 为向量赋值lambda = zeros(1, length(h);lambda(1) = 0;for i = 2 : length(lambda)lambda(i) = h(i) / (h(i - 1) + h(i);end% 为向量 d 赋值d = zeros(1, length(h) + 1);d(1) = 0;for i = 2 : length(h)a = ( DepVar(i) - DepVar(i - 1) ) / ( IndVar(i) -

14、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 赋值% 将主对角线上的元素全部置为 2 A = zeros( length(d), length(d) );for i = 1 : length(d) A(i, i) = 2;end% 将向量的各元素赋给主对角线右侧第一条对角线for i = 1 : length(

15、d) - 1A(i, i + 1) = lambda(i);end% 将向量 d 的各元素赋给主对角线左侧第一条对角线for i = 1 : length(d) - 1A(i + 1, i) = mu(i);end% 求解向量 MM =A d;% 求解每一段曲线的函数表达式for i = 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,

16、IndVar(i) : 0.01 : IndVar(i + 1);Coefs_2 = M(i + 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)2 / 6) / h(i);Part_3 = conv(Coefs_3, -1, IndVar(i + 1);S_3 = polyval

17、 (Part_3, IndVar(i) : 0.01 : IndVar(i + 1);Coefs_4 = (DepVar(i + 1) - M(i + 1) * h(i)2 / 6) / h(i);Part_4 = conv(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)% 在样条插值曲线的相应位置标注该段曲线的函数表

18、达式text(i, polyval(Part_1, 5),.itS, num2str(i),(x)= , num2str(Coefs_1),( , num2str( IndVar(i + 1) ),-x)3+, .num2str(Coefs_2),(x- , num2str( IndVar(i) ),)3+ , num2str(Coefs_3),.(, num2str( IndVar(i + 1) ),-x)+ , num2str(Coefs_4),(x- , num2str( IndVar(i) ),), .FontName, Times New Roman, FontSize, 14)ho

19、ldonend% 过x=2和x=4两个横轴点作垂线%line(2, 2, 4.5, 0.5, line(4, 4, 4.5, 0.5,LineStyle LineStyle, - );, - );% 为 x 轴和 y 轴添加标注xlabel(itx, FontName, Times New Roman, .FontSize, 14,FontWeight, bold );ylabel( its(x) Rotation, FontName , 0, FontSize, Times New Roman , 14, FontWeight, ., bold);最终，三次自然样条插值函数s(x)表达式为：

20、0.125(x1)3(2x)3.125( x1) , x1,2 ,s( x)0.0625( 4x)30.1875( x 2) 31.75(4x)2.75( x 2) , x 2, 4 ,0.375(5x)34.375(5x) 2( x4) ,x4,5 .曲线的图像如图所示：例 3课后习题与思考题第7 题i50953x23345.00000y i07408004725782055667.00000试求在区间 0.25,0.53上满足上述函数表所给出的插值条件的三次自然样条插值函数 s(x)求解出的三次样条插值函数将写成三弯矩方程的形式：M jM j 1( y jM j h2j )( y j 1M

21、 j 1hj2 )s( x)( x j 1x)3(x x j )36( x j 1x)6(x x j )6hj6h jhjhj本题采用和例 2 基本相同的 Matlab代码，只改变初始条件。最终，三次自然样条插值函数s(x)表达式为：6.2652( x 0.25)310( 0.3 x)10.9697 ( x0.25),x0.25,0.30 ,3.4806(0.39 x)31.5993(x 0.3) 36.1137(0.39x)6.9518( x0.3) ,x0.30 ,0.39 ,s( x)x) 32.859(x 0.39)310.417(0.45x)11.1903( x0.39) ,x0.3

22、9, 0.45 ,2.399(0.452.1442(0.53x)38.3987(0.53x)9.1( x0.45) ,x0.45, 0.53 .曲线的图像如图所示：例 4 课后习题与思考题第6 题求yx 的二次插值式P2(x)，使：P2 (100)10,P2 (121)11,P2 (144)12(144);并计算115 的近似值并估计误差。本题采用拉格朗日二次插值法进行计算：以下为 Matlab代码：%=% 本段代码解决课本第 2 章习题与思考题第 6 题%=clearallclc% 自变量 x 与因变量 y 的取值IndVar = 100, 121, 144;DepVar = 10, 11,

23、 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 = DepVar(2) /.( IndVar(2) - IndVar(1) ) * ( IndVar(2) - IndVar(3) )Part_2 = conv(Coefs

24、_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