已给sin0.32=0.314567,sin0.34=0.333487,sin0.36=0.352274,用线性插值及抛物插值计算sin0.3367的值并估计截断误差。
1. 线性插值
clc;clear; y=sin_L(0.32,0.314567,0.34,0.333487,0.3367); function y=sin_L(x0,y0,x1,y1,x) % sin_L输出sin(x)使用线性插值计算得到的函数值 % 例如: y=sin_L(0.32,0.314567,0.34,0.333487,0.3367) % y = % 0.3304 % R = % 9.1892e-06 % 以下为判断输入值是否合法的代码 if nargin~=5 %nargin表示number of function input arguments,判断函数输入参数的个数 error('请输入一次插值的2个插值节点和待插值点') end if ~( isnumeric(x0)&&isnumeric(y0)&&isnumeric(x1)&&isnumeric(y1)&&isnumeric(x) ) %isnumeric()检测字符串是否只由数字组成,如果字符串中只包括数字,就返回Ture,否则返回False。 error('输入参数必须是数') end % 核心计算的代码 %y=y0+(y1-y0)*(x-x0)/(x1-x0);%已知一次多项式,即直线上两点,可以用点斜式写出直线 y=y0*(x-x1)/(x0-x1)+y1*(x-x0)/(x1-x0); %两点式 % 以下为求解截断误差的代码 syms M2; % 因为sin(x)的二阶导数是本身,所以只需要挑出最大的y值,即可得到M2 %syms函数用于创建符号对象 if y0>y1 M2=y0; else M2=y1; end R=M2*abs((x-x0)*(x-x1))/2; disp(['插值结果是:y=',num2str(y)]); disp(['最大误差是:R=',num2str(R)]); end
2. 抛物插值
clc;clear; %y1=sin_L(0.32,0.314567,0.34,0.333487,0.3367);%一次性插值 y2=sin_T(0.32,0.314567,0.34,0.333487,0.36,0.352274,0.3367)%二次抛物插值 function y=sin_L(x0,y0,x1,y1,x) % sin_L输出sin(x)使用线性插值计算得到的函数值 % 例如: y=sin_L(0.32,0.314567,0.34,0.333487,0.3367) % y = % 0.3304 % R = % 9.1892e-06 % 以下为判断输入值是否合法的代码 if nargin~=5 %nargin表示number of function input arguments,判断函数输入参数的个数 error('请输入一次插值的2个插值节点和待插值点') end if ~( isnumeric(x0)&&isnumeric(y0)&&isnumeric(x1)&&isnumeric(y1)&&isnumeric(x) ) %isnumeric()检测字符串是否只由数字组成,如果字符串中只包括数字,就返回Ture,否则返回False。 error('输入参数必须是数') end % 核心计算的代码 %y=y0+(y1-y0)*(x-x0)/(x1-x0);%已知一次多项式,即直线上两点,可以用点斜式写出直线 y=y0*(x-x1)/(x0-x1)+y1*(x-x0)/(x1-x0); %两点式 % 以下为求解截断误差的代码 syms M2; % 因为sin(x)的二阶导数是本身,所以只需要挑出最大的y值,即可得到M2 %syms函数用于创建符号对象 if y0>y1 M2=y0; else M2=y1; end R=M2*abs((x-x0)*(x-x1))/2; disp(['插值结果是:y=',num2str(y)]); disp(['最大误差是:R=',num2str(R)]); end function y=sin_T(x0,y0,x1,y1,x2,y2,x) % sin_T输出sin(x)使用抛物插值计算得到的函数值 % 例如: y=sin_T(0.32,0.314567,0.34,0.333487,0.36,0.352274,0.3367) % y = % 0.3304 % R = % 2.0315e-07 % 以下为判断输入值是否合法的代码 if nargin~=7 error('请输入线性插值的插值节点和插值点') end if ~( isnumeric(x0)&&isnumeric(y0)&&isnumeric(x1)&&isnumeric(y1)&&isnumeric(x2)&&isnumeric(y2)&&isnumeric(x) ) error('输入参数必须是数') end % 核心计算的代码 y=y0*(x-x1)*(x-x2)/((x0-x1)*(x0-x2))+y1*(x-x0)*(x-x2)/((x1-x0)*(x1-x2))+y2*(x-x0)*(x-x1)/((x2-x0)*(x2-x1)); %{ % 以下为求解截断误差的代码 y_0=cos(x0); % 因为sin(x)的三阶导数是cos(x),那么只要求出x0,x1,x2的cos值,然后去最大即可得到M3 y_1=cos(x1); y_2=cos(x2); syms M3; if y_0>y_1 M3=y_0; else M3=y_1; end if y_2>M3 M3=y_2; end R=M3*abs((x-x0)*(x-x1)*(x-x2))/6 end %} % 以下为求解截断误差的代码 %sysm M3; M3=cos(x0);%因为在第一象限,cosx是减函数,所以节点处的最大值在第一个节点那里取得 R=M3*abs((x-x0)*(x-x1)*(x-x2))/6 end
参考来源:https://blog.csdn.net/m0_37395228/article/details/80874393
五,优点和缺点
拉格朗日插值法的公式结构整齐紧凑,在理论分析中十分方便,然而在计算中,当插值点增加或减少一个时,所对应的基本多项式就需要全部重新计算,
于是整个公式都会变化,非常繁琐。这时可以用重心拉格朗日插值法或牛顿插值法来代替。此外,当插值点比较多的时候,拉格朗日插值多项式的次数
可能会很高,因此具有数值不稳定的特点,也就是说尽管在已知的几个点取到给定的数值,但在附近却会和“实际上”的值之间有很大的偏差。
这类现象也被称为龙格现象,解决的办法是分段用较低次数的插值多项式。(https://www.cnblogs.com/fengzhiyuan/p/8645246.html)