一、Romberg求积公式的演化
根据前面一篇博客《基于复化梯度求积的求积步长自适应matlab实现 》,我们知道可以通过对积分区间不断进行二分,然后采用复化梯形公式计算得到新的求积结果。但是其存在的主要问题是,收敛速度太慢,前述的例子需要进行10次二分才能达到10^-7的精度,达到该精度的区间数为1024,共需节点1025个。示例为一简单函数,因此计算量相对还算不大,但是如果计算函数复杂,要求精度更高,那么基于复化梯度求积的求积步长自适应的方法也难以满足需求。为此,数学家提出了Romberg求积公式。
设Tn和T2n分别为二分前后利用复化梯形求积公式求得的积分近似值,根据式(1.1)可以计算得到式(1.2):
式(1.1)
式(1.2)
当Tn与T2n非常接近时,则可以保证T2n的误差非常小。这种直接利用计算结果来估计误差的方法称为误差的事后估计法。其思想为误差补偿思想。记
式(1.3)
而根据经验证得:
式(1.4)
即用复化梯形法求出的二分前后两个积分近似值Tn与T2n按照式(1.3)作线性组合,所得到的结果就是复化simpson公式求得的积分近似值Sn。
类似地,对Sn与S2n作线性组合可以得到Cn:
式(1.5)
对Cn和C2n作线性组合得到Romberg公式:
式(1.6)
二、Richardson外推技术的成型
由于Richardson推理过程比较复杂,此处只给出相关链接(https://en.wikipedia.org/wiki/Richardson_extrapolation),如果有兴趣,请读者自行查阅相关资料。下面给出通项公式:
式(1.7)
后面将会用到上式。
三、Romberg算法描述
1)赋初值,计算
式(1.8)
并将1→k(k)为二分次数。(注:这句话的意思是其实就是指计算T数表的第一列,即前面的自适应步长复化梯度求积)
2)求梯形值,按照参考I的式(1.3)计算
3) 外推计算,求加速值。按照式(1.7)依次求出T数表中第k行的其余元素
4)精度控制,对给定误差,若对角线上相邻元素满足下式(1.9),则停止计算,取
作为满足精度要求的近似值;否则k+1→k,转2)继续计算,直到满足要求。
式(1.9)
四、Romberg的Matlab实现
实现代码如下:
%% Romberg求积公式
function [T_result] = RichardsonExtrapolated(x_LowBound,x_UpBound,AccuracyValue)
% 2017-11-03 xh_scu 1270978696@qq.com
% inputs:
% ------x_LowBound:计算区间下届
% ------x_UpBound:计算区间上届
% ------AccuracyValue:精度要求
% outputs:
% ------result:近似积分结果,为T数表
% step_1:计算区间[x_LowBound,x_UpBound]上的复化梯度积分
step_length = x_UpBound - x_LowBound;
T_result(1,1) = step_length * (CalcuNoteValue(x_LowBound) + CalcuNoteValue(x_UpBound))/2;
count = 1;
while 1
% 步长减半
step_length = step_length/2;
% 初始化新插入值的和
sum_new_value = 0;
% 计算新插入值的和
%区间个数
TwoPowerCount = 2.^count;
%累加新增的节点函数值之和
% | |
% | | | 新增节点为:a+(b-a)/2
% | | | | | 新增节点为:a+(b-a)*1/4, a+(b-a)*3/4
% | | | | | | | | | 新增节点为:a+(b-a)*1/8, a+(b-a)*3/8, a+(b-a)*5/8,a+(b-a)*7/8
up_Bound = TwoPowerCount-1;
% 计算新增节点的函数值之和
for j = 1:2:up_Bound
sum_new_value = sum_new_value + CalcuNoteValue(x_LowBound + step_length*j);
end
% 更新T_result(k+1,1)(注:实际为式1.9的T(k+1,0)
T_result(count+1,1) = T_result(count,1)/2 + step_length*sum_new_value;
% 循环计算第k次二分后的第k行的所有元素
for m = 1:count
T_result(count+1,m+1) = (4^m)*T_result(count+1,m)/((4^m)-1) - T_result(count,m)/((4^m)-1);
end
% 判断精度是否满足要求,满足则跳出循环;否则进行继续计算
if abs(T_result(count+1,count+1) - T_result(count+1,count))<AccuracyValue
break;
else
count = count + 1;
end % if条件结束
end % for循环结束
end % 函数结束
计算函数此处依然选用f(x)=sinx/x举例,读者可以根据自己的需要做相应的修改。代码如下:
%% 计算函数
function [note_value] = CalcuNoteValue(x)
if x == 0
note_value = 1;
else
note_value = sin(x)/x;
end
end
五、测试与分析
设积分区间为[0,1],精度为0.0000000001,测试结果如下表示。
二分次数 | T0 | T1 | T2 | T3 | T4 |
0 | 0.920735492403948 | 0 | 0 | 0 | 0 |
1 | 0.939793284806177 | 0.946145882273587 | 0 | 0 | 0 |
2 | 0.944513521665390 | 0.946086933951794 | 0.946083004063674 | 0 | 0 |
3 | 0.945690863582701 | 0.946083310888472 | 0.946083069350917 | 0.946083070387222 | 0 |
4 | 0.945985029934386 | 0.946083085384948 | 0.946083070351379 | 0.946083070367260 | 0.946083070367181 |
六、结论
根据测试结果可以得出:
1、测试函数f(x)=sinx/x只需要4次二分外推就可以达到0.00000001的精度。
2、收敛速度远远快于自适应步长复化梯度积分。