三次样条插值作业题.docx

1、三次样条插值作业题例1 设为定义在0,3上的函数，有下列函数值表：xi0123yi00.521.5且，试求区间0,3上满足上述条件的三次样条插值函数本算法求解出的三次样条插值函数将写成三弯矩方程的形式：其中，方程中的系数，将由Matlab代码中的变量Coefs_1、Coefs_2、Coefs_3以及Coefs_4的值求出。以下为Matlab代码：%=% 本段代码解决作业题的例1%=clear allclc % 自变量x与因变量y，两个边界条件的取值IndVar = 0, 1, 2, 3;DepVar = 0, 0.5, 2, 1.5; LeftBoun = 0.2;RightBoun = -1

2、; % 区间长度向量，其各元素为自变量各段的长度h = zeros(1, length(IndVar) - 1);for i = 1 : length(IndVar) - 1 h(i) = IndVar(i + 1) - IndVar(i);end % 为向量赋值mu = zeros(1, length(h);for i = 1 : length(mu) - 1 mu(i) = h(i) / (h(i) + h(i + 1);endmu(i + 1) = 1; % 为向量赋值lambda = zeros(1, length(h);lambda(1) = 1;for i = 2 : length

3、(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);for i = 2 : length(h) a = ( DepVar(i) - DepVar(i - 1) ) / ( IndVar(i) - IndVar(i - 1) ); b = ( DepVar(i + 1) - DepVar(i) ) / ( Ind

4、Var(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) ) ) / h(i); % 为矩阵A赋值% 将主对角线上的元素全部置为2A = zeros( length(d), length(d) );for i = 1 : length(d) A(i, i) = 2;end% 将向量的各元素赋给主对角

5、线右侧第一条对角线for i = 1 : length(d) - 1 A(i, i + 1) = lambda(i);end% 将向量d的各元素赋给主对角线左侧第一条对角线for i = 1 : length(d) - 1 A(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, I

6、ndVar(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), . 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(

7、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 = 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) :

8、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)3+, . num2str(Coefs_2), (x-, num2str( IndVar(i) ), )3+, num2str(Coefs_3), . (, num2str( IndVar(i + 1) ), -x)+, num2str(Coefs

9、_4), (x-, num2str( IndVar(i) ), ), . FontName, Times New Roman, FontSize, 14) hold onend % 过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, Times New Roman, . FontSize, 14, FontWeight, bold);ylabel( its(x), FontName, Times

10、New Roman, .Rotation, 0, FontSize, 14, FontWeight, bold);最终，三次样条插值函数s(x)表达式为：曲线的图像如图所示：例2 已知函数值表：xi1245yi1342试求在区间1,5上满足上述函数表所给出的插值条件的三次自然样条插值函数本算法求解出的三次样条插值函数将写成三弯矩方程的形式：其中，方程中的系数，将由Matlab代码中的变量Coefs_1、Coefs_2、Coefs_3以及Coefs_4的值求出。以下为Matlab代码：%=% 本段代码解决作业题的例2%=clear allclc % 自变量x与因变量y的取值IndVar = 1,

11、 2, 4, 5;DepVar = 1, 3, 4, 2; % 区间长度向量，其各元素为自变量各段的长度h = zeros(1, length(IndVar) - 1);for i = 1 : length(IndVar) - 1 h(i) = IndVar(i + 1) - IndVar(i);end % 为向量赋值mu = zeros(1, length(h);for i = 1 : length(mu) - 1 mu(i) = h(i) / (h(i) + h(i + 1);endmu(i + 1) = 0; % 为向量赋值lambda = zeros(1, length(h);lamb

12、da(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) - IndVar(i - 1) ); b = ( DepVar(i + 1) - DepVar(i) ) / ( IndVar(i + 1) - IndVar(i) ); c = (b - a) / ( IndVar

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

14、d % 求解向量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, IndVar(i) : 0.01 : IndVar(i + 1); Coefs_2 = M(i + 1)/(6 * h(i); Part_2 = conv( Coefs_2, . c

15、onv( 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 (Part_3, IndVar(i) : 0.01 : IndVar(i + 1); Coefs_4 = (DepVar(i + 1) - M(i

16、+ 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) % 在样条插值曲线的相应位置标注该段曲线的函数表达式 text(i, polyval(Part_1, 5), . itS, num2str(i), (x)=, num2str(Coef

17、s_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) hold onend % 过x=2和x=4两个横轴点作垂线 %line(2, 2, 4.5, 0.5, LineStyl

18、e, -);line(4, 4, 4.5, 0.5, LineStyle, -); % 为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)表达式为：曲线的图像如图所示：例3 课后习题与思考题第7题xi0.250.300.390.450.53yi0.50000.54770

19、.62450.67080.7280试求在区间0.25,0.53上满足上述函数表所给出的插值条件的三次自然样条插值函数s(x)求解出的三次样条插值函数将写成三弯矩方程的形式：本题采用和例2基本相同的Matlab代码，只改变初始条件。最终，三次自然样条插值函数s(x)表达式为：曲线的图像如图所示：例4 课后习题与思考题第6题求的二次插值式，使：并计算的近似值并估计误差。本题采用拉格朗日二次插值法进行计算：以下为Matlab代码：%=% 本段代码解决课本第2章习题与思考题第6题%=clear allclc % 自变量x与因变量y的取值IndVar = 100, 121, 144;DepVar = 1

20、0, 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 = DepVar(2) /. ( ( IndVar(2) - IndVar(1) ) * ( IndVar(2) - IndVar(3) )

21、);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

22、.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)=, num2str(Coefs_1), (x-, num2str(IndVar(2), .)(x-, num2str(IndVar(3), )+, num2str(Coefs_2),

23、 (x-, num2str(IndVar(1), .)(x-, num2str(IndVar(3), )+, num2str(Coefs_3), (x-, num2str(IndVar(1), .)(x-, num2str(IndVar(2), ), .FontName, Times New Roman, FontSize, 14) % 为x轴和y轴添加标注xlabel( itx, FontName, Times New Roman, . FontSize, 14, FontWeight, bold);ylabel( itP_2(x), FontName, Times New Roman, . Rotation, 0, FontSize, 14, FontWeight, bold); % 求115平方根的估计值fun115 = polyval(Part_1, 115) + polyval(Part_2, 115) + polyval(Part_3, 115) % 求115平方根的准确值sqrt(115) % 求误差，单位为%est = 100 * (fun115 - sqrt(115) / fun115最终，拉格朗日插值函数的表达式为：从计算结果可以看出，的准确值为10.7238，而通过拉格朗日插值法求出的近似值为10.7228，误差为0.0098%的图像为：