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作线性组合得到Romb
c04b
erg公式：

    

  式（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、收敛速度远远快于自适应步长复化梯度积分。