数值积分方法之2——梯形法与外推法求近似积分

说明

Matlab的版本为Matlab R2019b；这篇笔记的全部内容是基于上课时老师布置的小小小作业；个人能力有限笔记当中难免会有错误，如有发现烦请指正，谢谢！。

要求

采用外推法和梯形法求$f(x)=x^3$、$f(x)=x^5$两个函数的积分${\int }_2^{4}f(x)dx$，其中$h$选用$\frac{1}{2}$。

1. 所谓的梯形法：$I\approx T(h)=\frac{h}{2}{\sum }_{i=0}^{n-1}[f(x_i)+f(x_{i+1})]$
2. 所谓的外推法：$I\approx T(h)+\frac{1}{3}[T(h)-T(2h)]$

Matlab实现

fun1 =@(x) x.^3;
fun2 =@(x) x.^5;
resFun1 = integral(fun1,2,4)
resFun2 = integral(fun2,2,4)
resFun1Tra = TraFun(fun1,0.5,2,4)
resFun2Tra = TraFun(fun2,0.5,2,4)
resFun1Exp = resFun1Tra + (resFun1Tra - TraFun(fun1,1,2,4))/3
resFun2Exp = resFun2Tra + (resFun2Tra - TraFun(fun2,1,2,4))/3

% 真正的计算已经算完了，下面是绘图表示
fun1 =@(x) x.^3;
fun2 =@(x) x.^5;
resFun1 = integral(fun1,2,4)
resFun2 = integral(fun2,2,4)
resFun1Tra = TraFun(fun1,0.5,2,4)
resFun2Tra = TraFun(fun2,0.5,2,4)
resFun1Exp = resFun1Tra + (resFun1Tra - TraFun(fun1,1,2,4))/3
resFun2Exp = resFun2Tra + (resFun2Tra - TraFun(fun2,1,2,4))/3

[xs,ys] = getXY(fun1,0.5,2,4);
myFigure = figure(1);
clf(myFigure)
myAxes = axes("Parent",myFigure);
myLine1 = line(myAxes,2:0.01:4,fun1(2:0.01:4));
myLine1.Color = 'r';
myLine1.LineStyle = '-';
myLine1.LineWidth = 0.5;
myLine2 = line(myAxes,xs,ys, ...
    'color','b','linestyle','--','linewidth',0.5);
text([2 2 2],[50 55 60],{['梯形法面积：' num2str(resFun1)],...
    ['外推法面积：' num2str(resFun1Tra)],['直接积分面积：' num2str(resFun1Exp)]},...
    'b',"FontSize",12,'FontName','宋体');
myAxes.FontName = "宋体";
myAxes.FontSize = 12;
myAxes.XLim =  [1.8 4.2];
myAxes.YLim =  [0 70];
myAxes.XLabel.String = "x";
myAxes.YLabel.String = "f(x)";
myAxes.Box = "on";
myAxes.LineWidth = 0.5;

[xs,ys] = getXY(fun2,0.5,2,4);
myFigure = figure(3);
clf(myFigure)
myAxes = axes("Parent",myFigure);
myLine1 = line(myAxes,2:0.01:4,fun2(2:0.01:4));
myLine1.Color = 'r';
myLine1.LineStyle = '-';
myLine1.LineWidth = 0.5;
myLine2 = line(myAxes,xs,ys, ...
    'color','b','linestyle','--','linewidth',0.5);
text([2 2 2],[760 830 900],{['梯形法面积：' num2str(resFun2)],...
    ['外推法面积：' num2str(resFun2Tra)],['直接积分面积：' num2str(resFun2Exp)]},...
    'b',"FontSize",12,'FontName','宋体');
myAxes.FontName = "宋体";
myAxes.FontSize = 12;
myAxes.XLim =  [1.8 4.2];
myAxes.YLim =  [0 1050];
myAxes.XLabel.String = "x";
myAxes.YLabel.String = "f(x)";
myAxes.Box = "on";
myAxes.LineWidth = 0.5;

% 自定义的梯形法积分函数
% fun：目标积分函数
% h：梯形间隔
% dow：积分下限
% upp：积分上限
% res：积分结果
function res = TraFun(fun,h,dow,upp)
x = dow:h:upp;
y = fun(x);
len = length(x);
res = 0;
for i=1:1:len-1
    res = res + h/2*(y(i)+y(i+1));
end
end

% 自定义函数用于获取绘图坐标
function [xx,yy] = getXY(fun,h,dow,upp)
x = dow:h:upp;
y = fun(x);
len = length(x);
xx = zeros(100,len*2-1);
yy = zeros(100,len*2-1);
for i=1:1:len*2-1
    if i<=len
        xx(:,i) = repelem(x(i),100)';
        yy(:,i) = linspace(0,y(i),100)';
    else
        linFun = @(t) y(i-len)+(t-x(i-len)).*(y(i-len+1)-y(i-len))./(x(i-len+1)-x(i-len));
        xx(:,i) = linspace(x(i-len),x(i-len+1),100)';
        yy(:,i) = linFun(xx(:,i));
    end
end
end

运行结果

${\int }_2^4{x}^{3}dx$	${\int }_2^4{x}^{5}dx$

梯形法与外推法对比

从图中可以看出使用外推公式极大地提高了精度。