计算动力学学习.pptx

《计算动力学学习.pptx》由会员分享，可在线阅读，更多相关《计算动力学学习.pptx（133页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、11.1 1.1 振动的概念1.2 1.2 单自由度系统自由振动1.3 1.3 单自由度系统强迫振动1.4 1.4 两个自由度系统的振动1.5 1.5 非线性振动概述第第1章章 绪论绪论第1页/共133页21.1 振动的概念振动：就是物体在静平衡位置附近所作的往复运动。振动系统：在振动问题中所研究的对象。如机器或结构物等。激励：外界对振动系统的作用或引起机器运动的力。响应：机器或结构在激励作用下产生的动态行为。第2页/共133页3振动的概念振动分析：研究振动系统、激励(输入)和响应(输出)三者之间的关系。第3页/共133页4力学基本模型振动系统的力学基本模型中包括三个基本“元件”：质量块、弹簧

2、和阻尼器。质量块:是物体惯性大小的度量。弹簧:表示振动系统弹性的理想模型。阻尼器:任何振动在没有外界干扰(激励)时都会逐渐消失，因此，系统存在一种阻碍振动持续进行的阻力，这种阻力称为阻尼。第4页/共133页5振动机理 任何结构，之所以能产生振动，是因为它本身具有质量（惯性力）和弹簧（恢复力）。从能量关系看,质量可以储存动能,弹簧可以储存势能（变形能）。振动就是动能和势能不断地转换。第5页/共133页61.2 单自由度系统单自由度系统:可以用一个独立坐标来确定系统的位置及其运动规律的振动系统;单自由度线性系统的振动是最简单的振动系统;许多实际问题可以足够精确地简化为单自由度振动系统;单自由度振动

3、系统的一些概念、特征和研究方法，是研究复杂振动系统的基础。第6页/共133页7m-k系统已知质量为m，弹簧的刚度系数为k。取质量的静平衡位置为坐标原点,当重物偏离 x 时,利用牛顿定律可得到运动微分方程：第7页/共133页8梁的横向振动质质量量为为m的的重重物物放放在在简简支支梁梁的的中中部部，不不计计梁梁的的质质量量。设设梁梁长长为为l，材材料料的的弹弹性性模模量量为为E，截截面面惯惯性性矩为矩为I。则利用材料力学的概念可得到：。则利用材料力学的概念可得到：d dst第8页/共133页9m-c-k系统已知质量为m，弹簧的刚度系数为k，粘性阻尼系数为c。运动微分方程为：第9页/共133页10m

4、-c-k系统令阻尼比为令阻尼比为则方程可写为则方程可写为令其解为令其解为代入方程得到代入方程得到此特征方程的两个根是此特征方程的两个根是第10页/共133页11大阻尼情况不同的阻尼比不同的阻尼比，对应的解的形式不同，运动，对应的解的形式不同，运动性质也不同。性质也不同。（1 1）11（大阻尼情况）（大阻尼情况）此时特征方程有两个不同的实根此时特征方程有两个不同的实根,通解为通解为第11页/共133页12大阻尼情况给出初始条件：给出初始条件：t0时时则可确定系数则可确定系数B和和D第12页/共133页13大阻尼情况 这这种种情情况况对对应应的的运运动动是是一一种种衰衰减减运运动动，但但不不是是我

5、我们们所所关关心心的的振振动动形形式式。设设x00，v00，则则运运动动图图形形大大致致如如下。下。第13页/共133页14临界阻尼情况（2）1（临界阻尼情况）（临界阻尼情况）此时特征方程有重根此时特征方程有重根利用初始条件确定常数为利用初始条件确定常数为 此时的阻尼系数称为临界阻尼系数，记为此时的阻尼系数称为临界阻尼系数，记为cc通解为通解为第14页/共133页15临界阻尼情况 临临界界阻阻尼尼情情况况也也是是一一种种非非振振动动形形式式的的衰衰减减运运动，按不同的初始条件其运动图形如下。动，按不同的初始条件其运动图形如下。第15页/共133页16小阻尼情况（3）0 1（小阻尼情况）（小阻尼

6、情况）此时特征方程有一对共轭复根，通解为此时特征方程有一对共轭复根，通解为或写为或写为利用初始条件确定出常数利用初始条件确定出常数第16页/共133页17小阻尼情况 解解中中有有两两个个因因子子，一一个个是是衰衰减减的的指指数数函函数数 ,它它将将使使振振幅幅越越来来越越小小，直直至至振振动最终消失动最终消失;另另一一个个是是正正弦弦函函数数 ，它表示系统以相同的周期通过平衡位置。它表示系统以相同的周期通过平衡位置。第17页/共133页18小阻尼情况因因此此系系统统呈呈现现为为一一种种衰衰减减形形式式的的等等周周期期振振动形式。动形式。第18页/共133页19小阻尼情况 单单自自由由度度粘粘性

7、性阻阻尼尼系系统统在在小小阻阻尼尼情情况况下下的的衰衰减减振振动动是是我我们们最最为为关关心心的的振振动动形形式式。这种衰减振动具有下列特性：这种衰减振动具有下列特性：（1）振幅衰减）振幅衰减 由由前前面面的的解解可可以以看看出出，振振幅幅不不再再是是常常量量，而是以几何级数而是以几何级数 快速衰减快速衰减;（2）等时性）等时性 系统仍以相同的周期通过平衡位置系统仍以相同的周期通过平衡位置;第19页/共133页20小阻尼情况（3）振动频率变小，周期变长）振动频率变小，周期变长 此时系统振动的频率和周期为：此时系统振动的频率和周期为：因因此此：衰衰减减振振动动的的固固有有频频率率比比无无阻阻尼尼

8、系系统统的的固固有有频频率率小小，振振动动周周期期变变大大，但但影影响响不不大大，特特别别是是当当阻阻尼尼很很小小（0zxM0第76页/共133页77梁单元单位体积变形能单位体积变形能胡克定律胡克定律第77页/共133页78梁单元单位体积变形能单位体积变形能单元总变形能单元总变形能第78页/共133页79梁单元轴向位移轴向位移应变应变第79页/共133页80梁单元第80页/共133页81梁单元横截面对横截面对y y轴的惯性矩轴的惯性矩第81页/共133页82梁单元插值插值第82页/共133页83梁单元刚度阵第83页/共133页84弹簧单元刚度阵刚度阵第84页/共133页85两自由度系统第85页

9、/共133页86刚度阵的建立第86页/共133页87刚度阵的建立第87页/共133页881.5 非线性振动概述非线性特性非线性特性 材料非线性材料非线性振幅过大超出材料线弹性范围几何非线性几何非线性位移或变形过大使结构几何形状显著变化非线性阻尼非线性阻尼材料内摩擦阻尼、流体阻尼等都是非线性阻尼负刚度负阻尼负刚度负阻尼有些情况下会存在负刚度和负阻尼非线性系统 当真实系统弹性元件的力与位移之间的关系超出线性范围，或阻尼元件当真实系统弹性元件的力与位移之间的关系超出线性范围，或阻尼元件的力与运动速度之间的关系不满足作线性关系时，系统的运动微分方程不能的力与运动速度之间的关系不满足作线性关系时，系统的

10、运动微分方程不能用线性微分方程描述，称系统为非线性系统。当真实系统作小运动时，可忽用线性微分方程描述，称系统为非线性系统。当真实系统作小运动时，可忽略系统的高阶微小量，近似地将系统看作线性系统。略系统的高阶微小量，近似地将系统看作线性系统。第88页/共133页89非线性振动概述非线性振动的研究方法 非线性振动研究的方法有：定性分析、定量分析和数值分析方法。非线性振动研究的内容 非线性振动研究的基本内容之一就是建立对真实振动系统的计算方法，改进计算精度，探索某些特殊现象的规律。定性法定性法 研究已知解的领域内系统的一般稳定性特征，而不是运动的时研究已知解的领域内系统的一般稳定性特征，而不是运动的

11、时间历程。通常采用间历程。通常采用几何方法几何方法描述系统的运动特征。描述系统的运动特征。定量法定量法 通过一些渐近的通过一些渐近的解析方法解析方法研究系统运动的时间历程。研究系统运动的时间历程。数值法数值法 通过通过数值计算数值计算方法研究系统非线性振动的规律和现象。方法研究系统非线性振动的规律和现象。第89页/共133页90非线性振动与线性振动的区别线性振动线性振动 非线性振动非线性振动 自由振动频率与初始条件无关自由振动频率与初始条件无关 自由振动频率与振幅有关自由振动频率与振幅有关 强迫振动频率与激励力频率相强迫振动频率与激励力频率相等等 强迫振动频率成分复杂，有时强迫振动频率成分复杂

12、，有时与激励频率不相等的频率成分与激励频率不相等的频率成分突出突出稳定平衡位置附近的运动是稳稳定平衡位置附近的运动是稳定的定的 稳定平衡位置附近具有多种稳定平衡位置附近具有多种稳定和不稳定运动稳定和不稳定运动强迫振动中每个激励频率强迫振动中每个激励频率有一个对应的振幅有一个对应的振幅 强迫振动中幅频与相频曲线强迫振动中幅频与相频曲线发生弯曲，产生多值性发生弯曲，产生多值性 叠加原理成立叠加原理成立 叠加原理不成立叠加原理不成立第90页/共133页91典型微分方程类型 单摆方程单摆方程库仑（库仑（Coulomb)Coulomb)摩擦振动方程摩擦振动方程第91页/共133页92典型微分方程类型 单

13、摆方程单摆方程库仑（库仑（Coulomb)Coulomb)摩擦振动方程摩擦振动方程第92页/共133页93典型微分方程类型 范德波（范德波（van der Polvan der Pol）方程）方程希尔希尔(HiIl)(HiIl)方程方程第93页/共133页94单摆 由牛顿第二定律：由牛顿第二定律：非线性方程非线性方程式中角频率：第94页/共133页95单摆 线性化处理线性化处理忽略3次以上的高次项得线性方程第95页/共133页96单摆 令代入方程得得特征方程：特征根：得通解为：式中 为复常数。由于描述单摆振动的应为实函数，所以常数 必须满足条件：将 写成指数形式后得：该式是振幅为P，角频率为

14、的简谐振动，其振动波形为正弦曲线。角角频频率只与摆线率只与摆线 l 得长度有关，与摆锤质量无关得长度有关，与摆锤质量无关，称为固有角频率。第96页/共133页97单摆周期与摆角无关？看看实验结果：定性结论：1.周期随摆角增加而增加周期随摆角增加而增加2.随摆角增加波形趋于矩形随摆角增加波形趋于矩形第97页/共133页98单摆周期数学表达式 对方程乘以 后积分其中 积分设t=0时，周期为 T，在 时应有 ，故有：最后得：第98页/共133页991.6 常微分方程初值问题数值解法一阶常微分方程初值问题一阶常微分方程初值问题 在区间在区间a x ba x b上的数值解法上的数值解法。(1)第99页/

15、共133页100数值方法的基本思想数值方法的基本思想100对常微分方程初值问题对常微分方程初值问题(1)式的数值解法，就是要算式的数值解法，就是要算出精确解出精确解y(x)y(x)在区间在区间 a,b 上的一系列离散节点上的一系列离散节点 处的函数值处的函数值 的近似值的近似值 。相。相邻两个节点的间距邻两个节点的间距 称为步长，步长可以称为步长，步长可以相等，也可以不等。本章总是假定相等，也可以不等。本章总是假定h为定数，称为为定数，称为定定步长步长，这时节点可表示为，这时节点可表示为第100页/共133页101数值方法的基本思想数值方法的基本思想101数值解法需要把连续性的问题加以离散化，

16、从而求出离散节点的数值解。对常微分方程数值解法的基本出发点就是离散化。其数值解法有两个基本特点，它们都采用“步进式”，即求解过程顺着节点排列的次序一步一步地向前 推 进，描 述 这 类 算 法，要 求 给 出 用 已 知 信 息 计算 的递推公式。建立这类递推公式的基本方法是在这些节点上用数值积分、第101页/共133页102中的导数中的导数 进行不同的离散化处理进行不同的离散化处理。对于初值问题对于初值问题数值方法的基本思想数值方法的基本思想102数值微分、泰勒展开等离散化方法，对初值问题第102页/共133页103的数值解法，首先要解决的问题就是如何对微分方程的数值解法，首先要解决的问题就

17、是如何对微分方程进行离散化，建立求数值解的递推公式。递推公式通进行离散化，建立求数值解的递推公式。递推公式通常有两类，一类是计算常有两类，一类是计算yi+1时只用到时只用到xi+1,xi 和和yi，即前一步的值，因此有了初值以后就可以逐步往下计即前一步的值，因此有了初值以后就可以逐步往下计算，此类方法称为算，此类方法称为单步法单步法；其代表是；其代表是龙格龙格库塔法。库塔法。另一类是计算另一类是计算y yi+1i+1时，除用到时，除用到x xi+1i+1,x,xi i和和y yi i以外，还要以外，还要用到用到 ，即前面，即前面k步的值，此类方法步的值，此类方法称为称为多步法多步法；其代表是亚

18、当斯法。；其代表是亚当斯法。数值方法的基本思想数值方法的基本思想103第103页/共133页104Euler公式公式1041 Euler公式公式 欧拉（欧拉（Euler）方法是解初值问题的最简单的）方法是解初值问题的最简单的数值方法。初值问题数值方法。初值问题的解的解y=y(x)代表通过点代表通过点 的一条称之为微分方的一条称之为微分方程的积分曲线。积分曲线上每一点程的积分曲线。积分曲线上每一点 的切线的切线的斜率的斜率 等于函数等于函数 在这点的值。在这点的值。第104页/共133页105Euler法的求解过程是法的求解过程是:从初从初始点始点P0(即点即点(x0,y0)出发出发,作积分曲线

19、作积分曲线y=y(x)在在P0点上点上切线切线 (其斜率为其斜率为 ),与与x=x1直线直线相交于相交于P1点点(即点即点(x1,y1),得到得到y1作为作为y(x1)的近似值的近似值,如上图所示。过点如上图所示。过点(x0,y0),以以f(x0,y0)为斜率的切线方为斜率的切线方程为程为 当当 时时,得得 Euler公式公式105第105页/共133页106同样同样,过过点点P1(x x1 1,y,y1 1),),作积分曲线作积分曲线y=y(x)y=y(x)的切线的切线交直线交直线x=x2于于P2点点,切线切线 的斜率的斜率 =直线方程为直线方程为Euler公式公式这样就获得了这样就获得了P

20、1点的坐标点的坐标。第106页/共133页107由此获得了由此获得了P P2 2的坐标。重复以上过程的坐标。重复以上过程,就可获得一系就可获得一系列的点列的点:P P1 1,P P1 1,P Pn n。对已求得点对已求得点 以以 =为斜率作直线为斜率作直线 Euler公式公式当 时,得 第107页/共133页108这样这样,从从x x0 0逐个算出逐个算出对应的数值解对应的数值解 Euler公式公式当 时,得取取第108页/共133页109Euler公式公式 从图形上看从图形上看,就获得了一条近似于曲线就获得了一条近似于曲线y=y(x)y=y(x)的折的折线线 。通常取通常取 (常数常数),)

21、,则则Euler法的计算格式法的计算格式 i=0,1,n (2)还可用数值微分、数值积分法和泰勒展开法推还可用数值微分、数值积分法和泰勒展开法推Euler格式。以数值积分为例进行推导。格式。以数值积分为例进行推导。将方程将方程 的两端在区间的两端在区间 上积分得，上积分得，第109页/共133页110 选择不同的计算方法计算上式的积分项选择不同的计算方法计算上式的积分项 ,就会得到不同的计算公式。就会得到不同的计算公式。(3)Euler公式公式用左矩形方法计算积分项用左矩形方法计算积分项 第110页/共133页111 代入代入(3)式式,并用并用yi近似代替式中近似代替式中y(xi)即可得到向

22、前欧拉即可得到向前欧拉（Euler）公式）公式 由于数值积分的矩形方法精度很低，所以欧拉由于数值积分的矩形方法精度很低，所以欧拉（Euler）公式当然很粗糙。）公式当然很粗糙。Euler公式公式第111页/共133页1122 梯形公式梯形公式 为了提高精度为了提高精度,对方程对方程 的两端在区间上的两端在区间上 积分得，积分得，改用梯形方法计算其积分项，即改用梯形方法计算其积分项，即 (4)代入代入(4)(4)式式,并用近似代替式中即可得到梯形公式并用近似代替式中即可得到梯形公式 梯形公式梯形公式第112页/共133页113 (5)式式的的右右端端含含有有未未知知的的yi+1,它它是是一一个个

23、关关于于yi+1的的函函数数方方程程,这这类类数数值值方方法法称称为为隐隐式式方方法法。相相反反地地,欧欧拉拉法法是是关关于于yi+1的的一一个个直直接接的的计计算算公公式式，这这类类数数值值方法称为显式方法。方法称为显式方法。梯形公式梯形公式 由由于于数数值值积积分分的的梯梯形形公公式式比比矩矩形形公公式式的的精精度度高高，因因此此梯梯形形公公式式（5）比比欧欧拉拉公公式式(2)的的精精度度高高一一个个数数值值方法。方法。(5)第113页/共133页1143 改进的欧拉公式改进的欧拉公式 显显式式欧欧拉拉公公式式计计算算工工作作量量小小,但但精精度度低低。梯梯形形公公式式虽虽提提高高了了精精

24、度度,但但为为隐隐式式公公式式,需需用用迭迭代代法法求求解解,计计算算工工作作量量大大。综综合合欧欧拉拉公公式式和和梯梯形形公公式式便便可可得得到到改改进进的的欧拉公式。欧拉公式。先先用用欧欧拉拉公公式式(2)求求出出一一个个初初步步的的近近似似值值 ,称称为预测值为预测值,它的精度不高它的精度不高,再用梯形公式再用梯形公式(5)对它对它改进的欧拉公式改进的欧拉公式第114页/共133页115改进的欧拉公式改进的欧拉公式预测预测 校正校正 （10）可可以以证证明明,公公式式(10)的的精精度度为为二二阶阶。这这是是一一种种一一步显式格式步显式格式,它可以表示为嵌套形式它可以表示为嵌套形式。校正

25、一次,即迭代一次,求得yi+1,称为校正值,这种预测-校正方法称为改进的欧拉公式：第115页/共133页116（11）或者表示成下列平均化形式（12）改进的欧拉公式改进的欧拉公式第116页/共133页1174 龙格龙格-库塔方法库塔方法1 龙格龙格-库塔库塔(Runge-Kutta)法的基本思想法的基本思想 Euler公式可改写成公式可改写成 则则yi+1的的表表达达式式y(xi+1)与与的的Taylor展展开开式式的的前前两两项项完完全相同全相同,即局部截断误差为即局部截断误差为 。改进的改进的Euler公式又可改写成公式又可改写成 第117页/共133页118 上述两组公式在形式上有一个共

26、同点:都是用f(x,y)在某些点上值的线性组合得出y(xi+1)的近似值yi+1,而且增加计算的次数f(x,y)的次数,可提高截断误差的阶。如欧拉公式:每步计算一次f(x,y)的值,为一阶方法。改进欧拉公式需计算两次f(x,y)的值，它是二阶方法。它的局部截断误差为 。龙格龙格-库塔法的基本思想库塔法的基本思想第118页/共133页119 于于是是可可考考虑虑用用函函数数f(x,y)在在若若干干点点上上的的函函数数值值的的线线性性组组合合来来构构造造近近似似公公式式，构构造造时时要要求求近近似似公公式式在在(xi,yi)处处的的Taylor展展开开式式与与解解y(x)在在xi处处的的Taylo

27、r展展开开式式的的前前面面几几项项重重合合，从从而而使使近近似似公公式式达达到到所所需需要要的的阶阶数数。既既避避免免求求偏偏导导,又又提提高高了了计计算算方方法法精精度度的的阶阶数数。或或者者说说,在在 这这一一步步内内多多预预报报几几个个点点的的斜斜率率值值，然然后后将将其其加加权权平平均均作作为为平平均均斜斜率率，则则可可构构造造出出更更高高精精度度的的计计算算格格式式，这这就就是是龙龙格格库库塔塔（Runge-Kutta）法法的的基基本本思想。思想。龙格龙格-库塔法的基本思想库塔法的基本思想第119页/共133页1202 二阶龙格二阶龙格库塔法库塔法 在在 上上取取两两点点xi和和 ,

28、以以该该两两点点处处的的斜斜率率值值k1和和k2的的加加权权平平均均(或或称称为为线线性性组组合合)来来求求取取平平均均斜率斜率k*的近似值的近似值K，即，即 式中式中:k1为为xi点处的切线斜率值，点处的切线斜率值，k2为为 点处的切线斜率值点处的切线斜率值,比照改进的欧拉比照改进的欧拉 法法,将将 视为视为 ，即可得，即可得 二阶龙格二阶龙格-库塔法库塔法第120页/共133页121二阶龙格二阶龙格-库塔法库塔法对常微分方程初值问题(1)式的解 y=y(x),根据微分中值定理，存在点 ，使得 式中式中 也即也即（13）K可看作是可看作是y=y(x)在区间在区间 上的平均斜率。所以上的平均斜

29、率。所以第121页/共133页122（14）将将y(xi)在在x=xi处进行二阶处进行二阶Taylor展开：展开：（15）二阶龙格二阶龙格-库塔法库塔法可得计算公式为：将将 在在x=xi处进行一阶处进行一阶Taylor展开：展开：第122页/共133页123将以上结果代入（将以上结果代入（14）得：）得：（16）二阶龙格二阶龙格-库塔法库塔法对式对式(15)和和(16)进行比较系数后可知进行比较系数后可知,只要只要 第123页/共133页124式式(17)中具有三个未知量中具有三个未知量,但只有两个方程但只有两个方程,因而有无因而有无穷多解。若取穷多解。若取 ,则则p=1，这是无穷多解中的，这

30、是无穷多解中的一个解，将以上所解的值代入式一个解，将以上所解的值代入式(14)并改写可得并改写可得 二阶龙格二阶龙格-库塔法库塔法（17）成立成立,格式格式(14)的局部截断误差就等于的局部截断误差就等于有2阶精度第124页/共133页125二阶龙格二阶龙格-库塔法库塔法 不难发现，上面的格式就是改进的欧拉格式。不难发现，上面的格式就是改进的欧拉格式。凡满足条件式（凡满足条件式（17）有一簇形如上式的计算格式，）有一簇形如上式的计算格式，这些格式统称为二阶龙格这些格式统称为二阶龙格库塔格式。因此改进的库塔格式。因此改进的欧拉格式是众多的二阶龙格欧拉格式是众多的二阶龙格库塔法中的一种特殊库塔法中

31、的一种特殊格式。格式。第125页/共133页126二阶龙格二阶龙格-库塔法库塔法 此计算公式称为变形的二阶龙格此计算公式称为变形的二阶龙格库塔法。式中库塔法。式中 为区间为区间 的中点。的中点。若取若取 ,则则 ，此时二阶龙格，此时二阶龙格-库塔库塔法的计算公式为法的计算公式为 第126页/共133页1273 三阶龙格三阶龙格-库塔法库塔法 为了进一步提高精度，设除为了进一步提高精度，设除 外再增加一点外再增加一点 并并用用三三个个点点 ,的的斜斜率率k1,k2,k3加加权权平平均均得得出出平平均均斜斜率率k*的的近近似似值值，这这时时计计算算格格式式具具有有形形式式:三阶龙格三阶龙格-库塔法

32、库塔法第127页/共133页128三阶龙格三阶龙格-库塔法库塔法于是可得于是可得 为了预报点为了预报点 的斜率值的斜率值k3,在区间在区间 内有两内有两个斜率值个斜率值k1和和k2可以用可以用,可将可将k1,k2加权平均得出加权平均得出 上的平均斜率上的平均斜率,从而得到从而得到 的预报值的预报值（18）第128页/共133页129三阶龙格三阶龙格-库塔法库塔法运运用用Taylor展展开开方方法法选选择择参参数数 ,可可以以使使格格式式(18)的的局局部部截截断断误误差差为为 ,即即具具有有三三阶阶精精度度，这这类类格格式式统统称称为为三三阶阶龙龙格格库库塔塔方方法法。下下列列是是其其中中的的

33、一种，称为库塔（一种，称为库塔（Kutta）公式。）公式。（19）第129页/共133页1304 四阶龙格四阶龙格库塔法库塔法 如如果果需需要要再再提提高高精精度度，用用类类似似上上述述的的处处理理方方法法，只只需需在在区区间间 上上用用四四个个点点处处的的斜斜率率加加权权平平均均作作为为平平均均斜斜率率k*的的近近似似值值，构构成成一一系系列列四四阶阶龙龙格格库库塔塔公公式式。具具有有四四阶阶精精度度，即即局局部部截截断断误误差差是是 。由于推导复杂，这里从略，只介绍最常用的一种由于推导复杂，这里从略，只介绍最常用的一种四阶经典龙格四阶经典龙格库塔公式。库塔公式。四阶龙格四阶龙格-库塔法库塔法第130页/共133页131四阶龙格四阶龙格-库塔法库塔法（20）第131页/共133页132第132页/共133页133感谢您的观看。第133页/共133页