数值计算——MATLAB数值积分原理详讲

巴山学长

发布于 2020-03-25 14:40:48

3.5K0

发布于 2020-03-25 14:40:48

文章被收录于专栏：巴山学长

大家的日常学习是一个循序渐进的过程，随着对问题的不断深入简单问题也会有新的发现。这不我们再来回头讲讲过冷水之前学习过程中遇到的数值积分的问题。对以下图像进行积分：只知道到图像点不知道函数解析形式。

显然这是一个简单的数值积分问题，但是过冷水会给大家分享简单问题吗？其必有玄妙，且听我道来。

1:γ=f(x)的表达式我们是不知道的，只知道一系列序列点，最好我们能够根据这一些列点求积分。实际在已知f(x)表达式的情况下，原函数我们也求不出来，所以该问题应该采用数值积分解决而不是符号表达式integrate。

2:采用数值积分实际无法积分整个区间，在采用quad()命令解决问题时，其值也和我们已知的积分值有出入，quad()为何不完美？

针对上述问题的原因，刚开始过冷水表示一无所知，所幸过冷水有幸接触到一份含金量较高的资料，遂解答了我的疑惑，这份资料就是：

随手打开就给予我很多灵感,告诉我遇到的问题是什么原因。Matlab提供的数值积分函数并不是真的直接给出该函数的数值积分，而是对所求函数处理后的积分。Matlab称为近似计算，而我们在实际应用中会误以为是精确结算，概念理论的混淆是借助软件进行学习研究的同学的大忌，以为现成软件可以解决你大部分疑惑，你竭尽全力都不能解决的问题，进行软件设计的人也不可以。

Matlab中无论被积函数是解析形式还是数表形式，其基本原理都是用多项式函数近似代替被积函数，用对多项式的积分结果近似代替被积函数的积分。现和大家分享最常用的三种插值型数值积分方法：矩形法、梯形法、抛物线法，多项式法。

矩形法

定积分的集合意义是算曲面梯形的面积，如果将区间[a,b]分成n等分，每个小区间上都是一个小的曲边梯形，用一个个小矩形代替这些小曲边梯形，然后把所有小矩形的面积加起来就近似等于整个曲边梯形的面积，于是便求出了定积分的近似值，这就是矩形法的基本原理。

假如f(x)在[0,1]上可积，则定积分的可以写为：

可知，当n充分大时，可将yn视为积分y的近似值。

注意！n充分大的条件只使用于γ=f(x)有解析形式，如果是数表形式，n值大小就受制于变量分组，可以适当采用插值方法增加分组。我的问题是在[0,1]区间上有50分组，所以我的N就是五十，由于无法满足n充分大的条件，所以我用该方法必然会和理想情况存在误差。由于我的积分上限是变量，所以采取累加的方式求和。

梯形法

将积分区间[a,b]n等分，用线段依次连接各分点，每段都形成一个小的直角梯形，如果用这些小直角梯形面积之和代替原来的小曲边梯形面积之和，就可以求得定积分的近似值。

解释：因为用plot绘图命令就是点与点的直线连接，所以梯形法和图像重合，实际是不重合的。

在第k个子区间

小曲边梯形面积近似为

则得到近似公式

它的实际含义是利用逐段线性函数作为f(x)的近似。

抛物线法

为了提高计算精度，可以用分段二次插值函数Sk代替f(x)。由于每段都要用到相邻两个小区间端点的三个函数值，所以小区间的数目必须是偶数。记n=2m，(k=0,1,2,....,m-1)在第k段的两个小区间上用三个节点(x2k,f2k)、(x2k+1,f2k+1)、(x2k+2,f2k+2)做二次插值函数Sk，然后积分可得

求m段之和就得到整个区间上的近似积分，有

这些公式看的让人疑惑。不如转化为能看的懂的代码比较实在，文末附代码。

三种方法所求积分图像为：

从这幅图可以看出：

1.三种方法的趋势都一致，存在差别。梯形法和抛物线法完全重合。

2三种方法实际都只是求了49组数据，有一组数据没有求，就是第50组数据。程序会有说明。

3.梯形法和抛物线发完全重合不是说梯形法和抛物线发有多好，恰好说明其中一种方法是有问题的，实际是抛物线发我们应用时存在问题。

抛物线的问题是程序要求f(x):我解析式，这样就可以求f((a+b)/2)的具体值，实际处理过程中我使用的是f((a+b)/2)=f(a/2)+f(b/2)，所以和图像重合。我们所使用的quad命令就是抛物线法的封装函数，所以格式quad(fun,a,b)。

关于三种方法的优劣，且听过冷水下回分享，经过定积分数值积分理论分析，发现方法并不过如此，还以为有多好，还是一种近似方法，和我的多项式拟合殊途同归，读者会问我也没讲多项式求积分的方法啊？你需要查看我数值优化—三种复杂函数数值积分方法实例演示。过冷水和大家的分享就这些，有疑问或者感兴趣的问题需要解答，可在下方留言，过冷水均会热心解答。

1，图像代码：

gamma=[-14.1452369500343,-13.1579895070512,-12.2306514479636,-11.3653820829428,-10.5658087276651,-9.83710288991493,-9.18572811603586,-8.61872641845847,-8.14252651948296,-7.76147200944519,-7.47648188269286,-7.28430036283952,-7.17758040223437,-7.14569504407307,-7.17591621789159,-7.25457011872977,-7.36791565061951,-7.50266354433915,-7.64617021000172,-7.78638692993272,-7.91164509498800,-8.01033924480252,-8.07054944046791,-8.07963157995146,-8.02380266713948,-7.88776093211056,-7.65441199039263,-7.30482488630092,-6.81861012074670,-6.17496361282653,-5.35457485825240,-4.34233090696472,-3.13020757038409,-1.71917397115223,-0.118944012655777,1.65463867483629,3.58349134739212,5.65084429796858,7.84342138202411,10.1521617348977,12.5719581415545,15.1009536468031,17.7397635162559,20.4907991346686,23.3577463593264,26.3451908130469,29.4583617836340,32.7029642427086,36.0850732253904,39.6113276349503];
x=[1.00000000000000e-06,0.0204091428367347,0.0408172856734694,0.0612254285102041,0.0816335713469388,0.102041714183673,0.122449857020408,0.142857999857143,0.163266142693878,0.183674285530612,0.204082428367347,0.224490571204082,0.244898714040816,0.265306856877551,0.285714999714286,0.306123142551020,0.326531285387755,0.346939428224490,0.367347571061225,0.387755713897959,0.408163856734694,0.428571999571429,0.448980142408163,0.469388285244898,0.489796428081633,0.510204570918367,0.530612713755102,0.551020856591837,0.571428999428572,0.591837142265306,0.612245285102041,0.632653427938776,0.653061570775510,0.673469713612245,0.693877856448980,0.714285999285714,0.734694142122449,0.755102284959184,0.775510427795918,0.795918570632653,0.816326713469388,0.836734856306122,0.857142999142857,0.877551141979592,0.897959284816327,0.918367427653061,0.938775570489796,0.959183713326531,0.979591856163265,0.999999999000000];
%矩形法
figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
plot(x,gamma,'LineWidth',2);
bar1 = bar(x,gamma,'ShowBaseLine','off','BarWidth',1);
baseline1 = get(bar1,'BaseLine');
set(baseline1,'Visible','off');
xlabel('$x$','FontWeight','bold','Interpreter','latex');
ylabel('$\gamma$','FontWeight','bold','Interpreter','latex');
xlim(axes1,[0 1]);
set(axes1,'FontSize',14,'FontWeight','bold','LineWidth',2);
%梯形法
figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
for i=1:49
    fill([x(i),x(i+1),x(i+1),x(i),x(i)],[0,0,gamma(i+1),gamma(i),0],'b')
    hold on
end
plot(x,gamma,'LineWidth',2);
stem(x,gamma,'MarkerFaceColor',[0.850980401039124 0.325490206480026 0.0980392172932625],'MarkerSize',3,'LineWidth',1);
xlabel('$x$','FontWeight','bold','Interpreter','latex');
ylabel('$\gamma$','FontWeight','bold','Interpreter','latex');
xlim(axes1,[0 1]);
set(axes1,'FontSize',14,'FontWeight','bold','LineWidth',2);

2，算法代码：

format long
gamma=[-14.1452369500343,-13.1579895070512,-12.2306514479636,-11.3653820829428,-10.5658087276651,-9.83710288991493,-9.18572811603586,-8.61872641845847,-8.14252651948296,-7.76147200944519,-7.47648188269286,-7.28430036283952,-7.17758040223437,-7.14569504407307,-7.17591621789159,-7.25457011872977,-7.36791565061951,-7.50266354433915,-7.64617021000172,-7.78638692993272,-7.91164509498800,-8.01033924480252,-8.07054944046791,-8.07963157995146,-8.02380266713948,-7.88776093211056,-7.65441199039263,-7.30482488630092,-6.81861012074670,-6.17496361282653,-5.35457485825240,-4.34233090696472,-3.13020757038409,-1.71917397115223,-0.118944012655777,1.65463867483629,3.58349134739212,5.65084429796858,7.84342138202411,10.1521617348977,12.5719581415545,15.1009536468031,17.7397635162559,20.4907991346686,23.3577463593264,26.3451908130469,29.4583617836340,32.7029642427086,36.0850732253904,39.6113276349503];
x=[1.00000000000000e-06,0.0204091428367347,0.0408172856734694,0.0612254285102041,0.0816335713469388,0.102041714183673,0.122449857020408,0.142857999857143,0.163266142693878,0.183674285530612,0.204082428367347,0.224490571204082,0.244898714040816,0.265306856877551,0.285714999714286,0.306123142551020,0.326531285387755,0.346939428224490,0.367347571061225,0.387755713897959,0.408163856734694,0.428571999571429,0.448980142408163,0.469388285244898,0.489796428081633,0.510204570918367,0.530612713755102,0.551020856591837,0.571428999428572,0.591837142265306,0.612245285102041,0.632653427938776,0.653061570775510,0.673469713612245,0.693877856448980,0.714285999285714,0.734694142122449,0.755102284959184,0.775510427795918,0.795918570632653,0.816326713469388,0.836734856306122,0.857142999142857,0.877551141979592,0.897959284816327,0.918367427653061,0.938775570489796,0.959183713326531,0.979591856163265,0.999999999000000];
inum=0;
for i=1:49
    a=x(i);
    b=x(i+1);
    fa=gamma(i);
    inum=inum+fa*(b-a);
    S(1,i)=inum;
end
%% 梯形法
%所有左点的数组
a=x(1,1:49);
%所有右点的数组
b=x(1,2:50);
%所有左点值
fa=gamma(1,1:49);
%所有右点值
fb=gamma(1,2:50);
%梯形面积
f=(fb+fa)./2.*(b-a);
%加和梯形面积
 S2=cumsum(f);
for i=1:49
 S(2,i)=S2(1,i);
end
%% 抛物线法
imum=0;
gamma=[-14.1452369500343,-13.1579895070512,-12.2306514479636,-11.3653820829428,-10.5658087276651,-9.83710288991493,-9.18572811603586,-8.61872641845847,-8.14252651948296,-7.76147200944519,-7.47648188269286,-7.28430036283952,-7.17758040223437,-7.14569504407307,-7.17591621789159,-7.25457011872977,-7.36791565061951,-7.50266354433915,-7.64617021000172,-7.78638692993272,-7.91164509498800,-8.01033924480252,-8.07054944046791,-8.07963157995146,-8.02380266713948,-7.88776093211056,-7.65441199039263,-7.30482488630092,-6.81861012074670,-6.17496361282653,-5.35457485825240,-4.34233090696472,-3.13020757038409,-1.71917397115223,-0.118944012655777,1.65463867483629,3.58349134739212,5.65084429796858,7.84342138202411,10.1521617348977,12.5719581415545,15.1009536468031,17.7397635162559,20.4907991346686,23.3577463593264,26.3451908130469,29.4583617836340,32.7029642427086,36.0850732253904,39.6113276349503];
x=[1.00000000000000e-06,0.0204091428367347,0.0408172856734694,0.0612254285102041,0.0816335713469388,0.102041714183673,0.122449857020408,0.142857999857143,0.163266142693878,0.183674285530612,0.204082428367347,0.224490571204082,0.244898714040816,0.265306856877551,0.285714999714286,0.306123142551020,0.326531285387755,0.346939428224490,0.367347571061225,0.387755713897959,0.408163856734694,0.428571999571429,0.448980142408163,0.469388285244898,0.489796428081633,0.510204570918367,0.530612713755102,0.551020856591837,0.571428999428572,0.591837142265306,0.612245285102041,0.632653427938776,0.653061570775510,0.673469713612245,0.693877856448980,0.714285999285714,0.734694142122449,0.755102284959184,0.775510427795918,0.795918570632653,0.816326713469388,0.836734856306122,0.857142999142857,0.877551141979592,0.897959284816327,0.918367427653061,0.938775570489796,0.959183713326531,0.979591856163265,0.999999999000000];
inum=0;
for i=1:49
    a=x(i+1);
    b=x(i);
    c=(a+b)./2;
    fa=gamma(i+1);
    fb=gamma(i);
    fc=(fa+fb)./2;
    inum=inum+(fb+4*fc+fa)*(a-b)/6;
    S(3,i)=inum;
end
figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
plot1 = plot(x(2:50),S,'LineWidth',1,'Parent',axes1);
set(plot1(1),'DisplayName','矩形法','LineWidth',2);
set(plot1(2),'DisplayName','梯形法','MarkerSize',8,'Marker','o');
set(plot1(3),'DisplayName','抛物线法','Marker','x','LineStyle','none');
xlabel('$x$','Interpreter','latex');
ylabel('$y$','Interpreter','latex');
title('不同数值积分方法比较');
set(axes1,'FontSize',20,'LineWidth',2);
legend1 = legend(axes1,'show');
set(legend1,'Position',[0.720839435964634 0.73500000461936 0.111273789999251 0.166111106491751]);

推荐指数：★★★☆ (7/10分)

好不好用只有用了才知道！若觉得好，别忘分享给和您一样爱学习研究的小伙伴哦！

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2020-03-24，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法