Richardson外推加速技术（含Romberg详细分析）的Matlab实现

一、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/8up_Bound = TwoPowerCount-1;% 计算新增节点的函数值之和for j = 1:2:up_Boundsum_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:countT_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))<AccuracyValuebreak;elsecount = count + 1;end % if条件结束
end % for循环结束
end % 函数结束

计算函数此处依然选用f(x)=sinx/x举例，读者可以根据自己的需要做相应的修改。代码如下：

%% 计算函数
function [note_value] = CalcuNoteValue(x)
if x == 0note_value = 1;
elsenote_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、收敛速度远远快于自适应步长复化梯度积分。