数值计算方法（三）——变步长梯形法与龙贝格算法

变步长梯形算法

提出背景：
复化求积公式虽然能提高精度，但需要给出步长，步长精度太大则精度低，步长太小则计算量大，难以找到一个合适的步长（划分成的小区间的个数）

算法描述：
1.对所有已存在的子区间进行二分化，区间数由n变为2n
2.利用区间数为n时的积分值Tn以及新增的节点（即原来各子区间的中点）递推出区间数为2n时的积分值T2n
3.利用两次计算结果的差来估计误差，直到满足精度
公式如下：

代码实现：
注意：
1.区间数是每次乘2，不是递增的
2.在递推的过程中，需要计算新增节点的函数和

/***@name Variable_step：利用变步长梯形公式求积*@param1 down：积分下限*@param2 upper：积分上限*@param3 limit_error：误差限
**/
void Variable_step(double down,double upper,double limit_error)
{Tn[0]=(function(down)+function(upper))/2;double h;       //定义小区间的长度double n=0.5;        //定义划分的区间数double S=0;      //新增和//注意每次区间二分，不是自加for (int i=0;i<19;i++){//从前先后推 Tn[0]->Tn[1] 每次对区间进行二分n=n*2;        //计算小区间的长度h=(upper-down)/n;    S=0;//计算小区间内所有中点的函数和 for (int j=0;j<n;j++){S+=function(down+j*h+h/2);}//迭代到下一项Tn[i+1]=Tn[i]/2+(h/2)*S;//如果检测到误差小于限制，则直接输出if(abs(Tn[i+1]-Tn[i])/3<limit_error)break;}
}

龙贝格算法

提出背景：
梯形公式收敛阶为2阶（余项为h的2次），利用变步长梯形公式时收敛速度较慢，故提出了龙贝格算法加速收敛

算法描述：
1.对得到的梯形序列Tn（1次代数精度）的相邻两项进行线性组合，即可得到抛物序列Sn（3次代数精度）
2.对抛物序列Sn的相邻两项进行线性组合，可得到柯斯序列Cn（5次代数精度）
3.对柯斯序列Cn的相邻两项进行线性组合，可得到龙贝格序列（7次代数精度）
示意图：

公式：

m对应为加速收敛而做的线性组合次数。

代码实现：
注意：
1.integral_table[10][4]为龙贝格积分表，从第0列到第3列分别为Tn，Sn，Cn，Rn
2.Tn（integral_table[i][0]）通过变步长梯形法求得
3.Sn，Cn，Rn由前一列的同行以及上一行的项进行线性组合得到

/***@name Romberg_arithmetic：利用龙贝格算法建立积分结果表*@param1 down：积分下限*@param2 upper：积分上限
**/
void Romberg_arithmetic(double down,double upper)
{double S=0;double n=0.5;double h=(upper-down);integral_table[0][0]=(function(down)+function(upper))*(h/2);//计算顺序  外循环由列开始迭代：1列 2列 3列（其中第0列由变步长梯形法求得）//                   内循环行迭代//变步长梯形法求第0列 Tnfor (int i=0;i<10;i++){n=n*2;h=(upper-down)/n;S=0;for (int j=0;j<n;j++){S+=function(down+j*h+h/2);}integral_table[i+1][0]=integral_table[i][0]/2+S*(h/2);}//计算Sn Cn Rnfor (int i=1;i<4;i++){for(int j=i;j<10;j++){//利用龙贝格求积公式integral_table[j][i]=(pow(4.0,i)*integral_table[j][i-1]-integral_table[j-1][i-1])/(pow(4.0,i)-1);}}
}

输出结果：