Adams 线性多步积分器（一）

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

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

\dfrac{dy}{dx} = 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\Sigma{p_{i}\dot{y}_{n-i}}, i =0,1,...,k
其中 pip_{i}是其系数,h为步长。
类似的，K阶的Moulton改正公式一般形式为：

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

y_{n+1} = y_{n}+h\Sigma{q_{i}\dot{y}_{n-i+1}}, i =0,1,...,k
其中 qiq_{i}是其系数。
其所谓的预测评估改正(PECE)计算方式是先用存储的前K步的一阶微分值计算出当前预测函数值，然后用当前函数值计算出当前的微分值，最后再代入Moulton公式计算最后的函数值，此过程可以迭代，但一般情况一次计算即可。其重点是计算K阶的系数 pip_{i}和 qiq_{i}。忽略系数的相关推导公式，直接给出其代码。
其中m_order为相应的阶数，m_bC和m_mC分别是bashforth和moulton公式的系数。函数 GMath::nchoosek(n, i)是计算组合数 CknC^{k}_{n}的函数。

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 coefficientsfor( 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 coefficientsfor( 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 线性多步积分器（一）相关推荐

分析步中的时间以及线性摄动中的时间?
1.在ABAQUS中的分析中有两种时间尺度,分别是总体时间和分析步时间. 总体时间始终贯穿于分析步中,总是在增长,是每个常规分析步时间的总和. 每一个分析步有自己的时间尺度,称之为分析步时间,对于每一 ...
四阶显式Adams法求方程组C语言,第五讲第4章线性多步法(续
<第五讲第4章线性多步法(续>由会员分享,可在线阅读,更多相关<第五讲第4章线性多步法(续(15页珍藏版)>请在人人文库网上搜索. 1.第4章线性多步法4.1 线性多步法的一 ...
将SolidWorks中的装配体导入到Adams中
如何将SolidWorks中的装配体导入到Adams中第一步:先将solidworks零件或者装配体保存成x_t格式或者另存为x_t格式.将输出的文件放到非中文(注意是非中文路径)路径下,比如D盘f ...
[渝粤题库]西北工业大学电力系统自动装置
电力系统自动装置一. 单选题 (共21题,共42分) 发电机调速系统的失灵区( ). (2分) A.越大越好 B.越小越好 C.要合适同步发电机无功调节特性特性下倾,对应的调差系数( ). (2分 ...
android电力监控平台,电力系统监控实验平台（网络型）
电力系统监控实验平台(网络型) 一.概述 "电力系统监控实验平台"是根据实验教学的要求,并结合工业实际应用和发展而研发的综合型实验平台.它能实现数据的集中监控和设备的自动控制.本 ...
电力系统监控实验平台QY-PGD20
随着电力系统规模不断扩大,系统发生故障的影响也越来越大,尤其大区域联网背景下的电力系统故障将会给经济.社会造成重大经济损失,因此保证电力系统安全稳定运行是电力生产的首要任务.电力系统是一个复杂的动态系 ...
《ABAQUS有限元分析实例详解》学习笔记_51CAE_新浪博客
石亦平老师的<ABAQUS有限元分析实例详解>当属ABAQUS学习的经典著作,一边认真拜读一边在此写下点滴笔记,供自己参考. 1.ABAQUS/CAE并没有自己专用的量纲系统,用户建立的整 ...
《电力系统自动装置》
试写出至少4种电力系统自动装置. 从硬件组成看,电力系统自动装置的结构形式有哪四种? 某电力系统负荷功率为6000MW,负荷调节效应系数KL*=2.设系统可能出现最大功率缺额为2000MW.希望系统恢 ...
运算放大器积分电路及积分电路设计
运算放大器积分电路及积分电路设计积分电路 Integrator circuit 在运算放大器积分器电路中,电容器插入反馈环路中,并在反相输入端与R1一起产生一个RC时间常数. 积分的物理意义积分的 ...
Python Statistical model: 用Scipy做线性回归分析（一步实现p value、标准差计算）
点赞发Nature 关注中Science p-value表示线性回归中的斜率.截距的显著性,p-value越小,表示得到的结果越显著. 之前我用curve fit的方法实现了线性函数的拟合,但是这有些 ...

Adams 线性多步积分器（一）

Adams 线性多步积分器（一）

Adams 线性多步积分器（一）相关推荐

最新文章

热门文章