Adams 线性多步积分器（一）

Adams 线性多步积分器（一）

在动力学系统仿真和计算中经常非线性微分方程的解算问题，其基本问题可以描述为：

常用的解算方法主要可以分为单步法和多步法，在轨道积分中常用的单步法有4阶Runge-Kutta和嵌套的RKF(7,8) 等方法，单步法相对来说解算结果稳定，但是对于步长要求严格，而且每一步都需要多次计算其右函数 f(x,y) ，以RKF78为例，每一步积分需要计算13次右函数值，而通常在轨道积分中，函数 f(x,y) 非常复杂计算也特别耗时，因此对于轨道积分而言，单步法效率不高。
而对于多步法而言，每一步只需要计算一次右函数值，因此其效率相对来说较高。但是多步法需要前K步的一阶微分值作为起步数据，通常可以借助RKF单步法起步。在多步法中较为常用的有Adams系列，BDF系列。本文重点介绍Adams的代码实现，其相关理论公式推导可以参考相关书籍。Adams积分器分为预测公式(Bashforth)和修正公式(Moulton)，K阶的Bashforth 的公式可以表达为：

其中 pi 是其系数,h为步长。
类似的，K阶的Moulton改正公式一般形式为：

其中 qi 是其系数。
其所谓的预测评估改正(PECE)计算方式是先用存储的前K步的一阶微分值计算出当前预测函数值，然后用当前函数值计算出当前的微分值，最后再代入Moulton公式计算最后的函数值，此过程可以迭代，但一般情况一次计算即可。其重点是计算K阶的系数 pi 和 qi 。忽略系数的相关推导公式，直接给出其代码。
其中m_order为相应的阶数，m_bC和m_mC分别是bashforth和moulton公式的系数。函数 GMath::nchoosek(n, i)是计算组合数 Ckn 的函数。

void GAdams::getCoef()
{
double *c = new double[m_order];
double *g = new double[m_order];
c[0] = 1.0;
for( int n = 1 ; n< m_order; n++ )
{
double s = 0.0;
for( int i =0 ; i<= n-1; i++)
{
s += c[i]/(n+1-i);
}
c[n] = -s;
}
for( int n = 0 ; n< m_order ; n++)
{
double s = 0.0;
for(int k = 0 ; k<=n; k++)
{
s += c[k];
}
g[n] = s;
}
//starting the bashforth coefficients
for( int n = m_order-1 ;n >= 0 ; n-- )
{
int sign = -1;
for( int i = 0 ; i<=n; i++ ) // m_order
{
sign *= -1;
m_bC[i] += sign*g[n]*GMath::nchoosek(n, i);
}
}
//starting the moulton coefficients
for( int n = m_order-1 ;n >= 0 ; n-- )
{
int sign = -1;
for( int i = 0 ; i<=n; i++ ) // m_order
{
sign *= -1;
m_mC[i] += sign*c[n]*GMath::nchoosek(n, i);
}
}
if( c!= NULL) {delete[] c; c = NULL;}
if( g!= NULL) {delete[] g; g = NULL;}
int testc = 0;
}
//calculate combination : Cnm = n!/m!/(n-m)!
double GMath::nchoosek(int n, int m)
{
double res = 0.0;
if( m == 0 )
{
return 1;
}
if( m == n )
{
return 1;
}
if( m > n )
{
return 0.0;
}
if ( m > n/2.0 )
{
m = n - m ;
}
double s1 =0.0, s2 =0.0;
for(int i = m+ 1; i<=n; i++)
{
s1 += log((double)i);
}
for( int
i = 2; i<= n-m ; i++)
{
s2 += log((double)i);
}
res = exp(s1-s2);
if( res < (numeric_limits<long>::max)() )
{
res =
static_cast<long>(res + 0.5);
}
return res;
}

最后

以上就是自觉期待为你收集整理的Adams 线性多步积分器（一）的全部内容，希望文章能够帮你解决Adams 线性多步积分器（一）所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。