Adams线性多步积分器应用详解

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

dydx=f(x,y);y(x0)=y0     d y   d x   =f(x,y);y( x 0 )= y 0  

  常用的解算方法主要可以分为单步法和多步法，在轨道积分中常用的单步法有4阶Runge-Kutta和嵌套的RKF(7,8) 等方法，单步法相对来说解算结果稳定，但是对于步长要求严格，而且每一步都需要多次计算其右函数 f(x,y)  f(x,y) ，以RKF78为例，每一步积分需要计算13次右函数值，而通常在轨道积分中，函数 f(x,y)  f(x,y) 非常复杂计算也特别耗时，因此对于轨道积分而言，单步法效率不高。
 

  而对于多步法而言，每一步只需要计算一次右函数值，因此其效率相对来说较高。但是多步法需要前K步的一阶微分值作为起步数据，通常可以借助RKF单步法起步。在多步法中较为常用的有Adams系列，BDF系列。本文重点介绍Adams的代码实现，其相关理论公式推导可以参考相关书籍。Adams积分器分为预测公式(Bashforth)和修正公式(Moulton)，K阶的Bashforth 的公式可以表达为：
 

yn+1=yn+hΣpiy˙n−i,i=0,1,...,k   y  n + 1  = y  n  +hΣ  p  i      y ˙    n − i   ,i=0,1,...,k

  其中 pi   p  i   是其系数,h为步长。
 

  类似的，K阶的Moulton改正公式一般形式为：
 

yn+1=yn+hΣqiy˙n−i+1,i=0,1,...,k   y  n + 1  = y  n  +hΣ  q  i      y ˙    n − i + 1   ,i=0,1,...,k

  其中 qi   q  i   是其系数。
 

  其所谓的预测评估改正(PECE)计算方式是先用存储的前K步的一阶微分值计算出当前预测函数值，然后用当前函数值计算出当前的微分值，最后再代入Moulton公式计算最后的函数值，此过程可以迭代，但一般情况一次计算即可。其重点是计算K阶的系数 pi   p  i   和 qi   q  i   。忽略系数的相关推导公式，直接给出其代码。
 

  其中m_order为相应的阶数，m_bC和m_mC分别是bashforth和moulton公式的系数。函数 GMath::nchoosek(n, i)是计算组合数 Ckn   C  n   k   的函数。
 

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;
    }

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删