变步长梯形算法
提出背景:
复化求积公式虽然能提高精度,但需要给出步长,步长精度太大则精度低,步长太小则计算量大,难以找到一个合适的步长(划分成的小区间的个数)
算法描述:
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列 Tn
for (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 Rn
for (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);
}
}
}
输出结果: