第３６卷 ２Ｏ１４年４月 湖州师范学院学报 Ｊｏｕｒｎａｌ ｏｆ Ｈｕｚｈｏｕ Ｕｎｉｖｅｒｓｉｔｙ Ｖｏ１．３６ Ａｐｒ．。２０１４ 两类基于ＭＡＴＬＡＢ的 非线性微分方程数值解的算法研究 薛亚宏 （甘肃工业职业技术学院电信学院，甘肃天水７４１０２５） 摘要：通过对工程动态控制及计算机仿真中有重要应用的两类非线性微分方程数值解的数学算法分析，建立四 阶定步长Ｒｕｎｇｅ—Ｋｕｔｔａ及Ｌｏｒｅｎｚ模型数值解的ＭＡＴＬＡＢ算法结构，讨论了变步长情形下的误差控制，绘制了 基于ＭＡＴＬＡＢ的Ｌｏｒｅｎｚ系统数值解在二维和三维空间下的图形，提出了在可接受误差限内的数值解检验的基 本思路． 关键词：ＭＡＴＬＡＢ；非线性微分方程；数值解；算法；误差 中图分类号：Ｇ７１ ＭＳＣ（２０００）：６５Ｌ１４ 文献标识码：Ａ 文章编号：１００９—１７３４（２０１４）０４—００１９—０５ ０ 引言 “数值分析”是计算数学的主要内容，其研究范围几乎涉及数学科学的一切分支，反过来很多数学分支 都需要探讨适用计算机的数值方法．在数值计算理论中，一般非线性微分方程的解析求解相当困难，在实 际应用和理论分析中往往求其数值解．同时，非线性微分方程的数值解问题及算法实现在工程动力系统控 制、人口统计、经济学、计算机仿真、数学建模等领域有十分重要的应用．本文将研究一般非线微分方程问 题数值求解的算法、误差、步长等问题，以及四阶定步长Ｒｕｎｇｅ—Ｋｕｔｔａ及其改进型的算法和Ｌｏｒｅｎｚ模型 数值求解等较为系统的方法，给出基于实变量运算结果的图形界面． １ 非线性微分方程的数值算法概述及误差与步长 １．１ 一般微分方程数值解概述 一般微分方程的数值解法多数是关于微分方程初值问题的数值解法，这类问题用一阶显式微分方程 主（￡）一ｆ（ｔ，ｚ（￡））， 组可描述为： 其中：ｚ （￡）＝［ （￡），ｚ ｚ（ｔ），…，ｚ （￡）］为状态向量；ｆ （・）一Ｅｆ （・），ｆｚ（・），…，ｆ （・）］为任意非线性 函数ｕＪ． 在初值问题中，若已知初始状态ｚ。一Ｉｘ （Ｏ），… （０）］ ，用数值求解法求出在ｔ∈ＥＯ，ｔｆ］内各个变 量状态 （￡）的数值，ｔ，又称为终止变量． １．２ 多元非线性微分方程数值算法及误差与步长 １．２．１ 初值问题的数值算法 对于多元非线性微分方程初值问题，Ｅｕｌｅｒ算法是较其它算法更易理解和控制的一类，在解决复杂微 分方程数值求解中被广泛采用．下面基于ＭＡＴＬＡＢ以Ｅｕｌｅｒ算法为例，给出非线性微分方程初值问题的 数值算法． ＊ 收稿日期：２０１４—０２—１０ 通讯作者：薛亚宏，讲师，硕士，研究方向：数学与应用数学．Ｅ—ｍａｉｌ：２３５３８６６３６０＠ｑｑ．ｃｏｒｎ ２０ 湖州师范学院学报 第３６卷 设已知在ｔ ｏ时刻系统状态向量的值为 （￡。），若选择一个较小的计算步长ｈ，则微分方程左侧的导数 近似为： ｚ（ｔ ｏ＋ｈ）一ｚ（￡０） （ ｏ＋ｈ）～ｔ。 。 代入微分方程得时刻方程的近似解为： 主（ ０＋ｈ）一ｚ（￡０）一ｈｆ（ｔ ０，．７Ｃ（￡０））． 由于近似解误差的存在，可将ｔ。＋ｈ时刻系统状态向量的值表示为： ｚ（ｔ ０＋ｈ）一ｚ（ｔ ｏ）＋ｈｆ（￡，ｚ（￡。））＋Ｒ ０， 或记为．７Ｃ 一ｚ（￡。＋ｈ），则主 一王（ 。＋ ）为系统状态向量在ｔ。＋＾时刻的近似值，即数值解 引．其中Ｒ。为 数值解的舍入误差． 一般地，设在ｔ 时刻系统的状态向量为ｚ ，则在ｔ ＋ｈ时刻Ｅｕｌｅｒ算法的数值解可写为： ＋１一ｚ ＋ｈｆ（￡， ＾）． 对于上述一般数值的求解表达式，用迭代法可由给定的初值问题逐步求出在ｔ∈Ｅｏ，Ｔ３内各时刻 ｔ。＋ ＾（忌一１，２，…， ）处的原始数值解． １．２．２数值解的精度控制及变步长影响因素 变步长对精度与速度会产生直接的影响．提高数值解精度的一种有效方法是减小步长值，但步长值连 续递减会产生以下问题： （１）降低运算速度．对选定的求解时限，减少步长等价于增加在固定时间区段内的计算点的数量，直 接导致了计算速度的急剧下降． （２）扩大累积误差．由于数值解固有误差，既使选择很小的步长，所得的数值解都将伴随一个舍入误差． 同时，减小步长会增加循环次数，从而使舍入误差的叠加的传递次数增多，最终形成较大的累积误差 ］． 舍入误差、累积误差和总误差关系如图１所示． 计算 误差 图１误差示意图 Ｆｉｇ．１ Ｓｃｈｅｍａｔｉｃ ｅｒｒｏｒ 步长ＪｊＩ 图２变步长示意图 Ｆｉｇ．２ Ｓｃｈｅｍａｔｉｃ ｖａｒｉａｂｌｅ ｓｔｅｐ ｓｉｚｅ 在多元非线性微分方程中，步长与误差是一对矛盾，对因变步长产生的误差进行有效控制也是非线性 微分方程数值求解的一项重要内容． 根据实际问题的算法和精度要求，可采用以下途径： （１）选择适当的步长．采用类似Ｅｕｌｅｒ的简单算法时，步长要适中，遵循宁小勿大的原则。 （２）改进近似算法的精度．由于Ｅｕｌｅｒ算法只是将原积分问题转为梯形法近似求解，因其精度较低而 不能有效逼近原始问题．可采用复合Ｓｉｍｐｓｏｎ法或更精确的样条插值法取代Ｅｕｌｅｒ算法，其中以Ｒｕｎｇｅ— Ｋｕｔｔａ法、Ａｄａｍｓ法最为常见． （３）采用变步长法．在多元微分方程数值解运算中，很多方法可通过变步长求解，若误差较小则增大 步长，若误差较大则自动减小步长，从而精确、有效地求解微分方程的初值问题． 厶 一般变步长算法的原理如图２所示．已知ｔ 时刻的状态变量ｚ ，分别计算步长为ｈ及 下的状态变 厶 量二川．如果两种步长下的误差ｅ—Ｉｌ主川一ｚ川ｌｌ小于给定的误差限，则可采用该步长 ；如果误差较 第４期 薛亚宏：两类基于ＭＡＴＬＡＢ的非线性微分方程数值解的算法研究 ２１ 大，则逐步减小步长，直到误差减小到允许范围内．这种自适应变步长算法能同时解决两个问题，即运算速 度和计算精度． ２ 两类一阶微分方程组数值解的算法与函数 求解一阶显式微分方程组的关键是用ＭＡＴＬＡＢ语言编写一个函数，描述需要求解的微分方程组，该 函数在程序头文件人口应表示为： ｆｕｎｃｔｉｏｎ ｘｄ—ｆｕｎｎａｍｅ（ｔ，ｘ）； ｆｕｎｃｔｉｏｎ Ｘｄ—ｆｕｎｎａｍｅ（ｔ，Ｘ，ｆｌａｇ，Ｐｌ，Ｐ２，…）． 在非线性微分方程求解中，有时需要对算法及控制条件进一步设置，这可以通过求解过程中的ｏｐｔｉｏｎｓ 变量进行修改，其方式初始变量可由ｏｄｅｓｅｔ（）函数获取．若相对误差设置较小时，则还需要给出以下声明： ｏｐｔｉｏｎｓ＝ｏｄｅｓｅｔ（’Ｒｅｌｔｏｌ’，ｌｅ一７）； ｏｐｔｉｏｎｓ＝－ｏｄｅｓｅｔ；ｏｐｔｉｏｎｓ．ＲｅｌＴｏｌ＝ｌｅ一７． 在实际应用中，往往需要定义一些附加参数，可由ｍ ，ｍ ，…，ｍ 表示．下面对数值分析中的两类重 要算法进行分析，并建立一般化的程序结构． ２．１ 四阶定步长Ｒｕｎｇｅ—Ｋｕｔｔａ算法及改进 四阶定步长Ｒｕｎｇｅ—Ｋｕｔｔａ算法是数值分析和系统仿真理论中的重要算法．其数学算法需先定义４个 附加向量： ｋ１一ｈｆ（ｔ ，ｚ ）， 是 一 ｃ ＋ ，ｚ ＋ ， 志。一 ＋ ，ｚ ＋ ， ｋ ４一ｈｆ（ｔ ＋ｈ，ｚ ＋忌３）． 其中：Ｓｉｍｐｓｏｍ为计算步长，在实际应用中是一个常数．然后由Ｒｕｎｇｅ—Ｋｕｔｔａ算法求解出下一个步长的 状态变量值为： １ ｚ ＋ｌ—ｚ ＋÷（愚１＋２ｋ ２＋２ｋ ３＋ｋ ４）． ｏ 对于上述表达式，用迭代法求出ｔ∈Ｅｔ。，ｔ ｒ］内各时刻点ｔ。＋ｈ，ｔ。＋２ｈ，…处的数值解． 基于数学上的算法，在ＭＡＴＬＡＢ中，可通过一组循环语句实现求解，其主体如下： ｆｕｎｃｔｉｏｎ［ｔｏｕｔ，ｙｏｕｔ］：ｒｋ一４（ｏｄｅｆｉｌｅ，ｔｓｐａｎ，ｙ０） ｔｓ＝ｔｓｐａｎ；ｔＯ＝ｔｓ（１）；ｔｆ＝ｔｓ（２）；ｙｏｕｔ：［］；ｔｏｕｔ＝［］；ｙ０＋ｙ０（：）； ｉｆ ｌｅｎｇｔｈ（ｔｓ）＝一３，ｈ—ｔｓ（３）；ｅｌｓｅ，ｈ＿－（ｔｓ（２）一ｔｓ（１））／１００；ｔｆ＝ｔｓ（２）；ｅｎｄ ｆｏｒ ｔ一［ｔＯ：ｈ：ｔｆ］ ｋｌ—ｈ＊ｆｅｖａｌ（ｏｄｅｆｉｌｅ，ｔ，ｙ０）；ｋ２，ｋ３，ｋ４一＊＊＊ ｙ０：ｙ０＋（ｋｌ＋２＊ｋ２＋２＊ｋ３＋ｋ４）／６；ｙｏｕｔ一［ｙｏｕｔ；ｙＯ’］；ｔｏｕｔ＝［ｔＯｕｔ；ｔ］； ｅｎｄ 以上ＭＡＴＬＡＢ算法中，变量ｔｓｐａｎ有两种构成方法，即等距时间向量和基于ｔ。，ｔ，，ｔ 的变步长向量． 在函数调用结束后，时间向量与状态变量构成的矩阵分别由程序中的ｔｏｕｔ和ｙｏｕｔ返回． 上述算法的优点是结构简洁、运算速度快，局限性是精度不高．为保证更高的精度和数值的稳定性，将 上述算法在步长递增方面进行改进，即在每个计算步长内对ｆ （・）函数进行６次求值，该算法称为四阶五 级ＲＫＦ算法． 对当前步长ｈ ，可以定义６个变量： ｉ－－１ 忌 一ｈ ｋｆ（ｔ ｋ＋ａｌｈ，，Ｘｋ＋ ｔ３．ｋ，），ｉ一１，２，…，６． ２２ ６ 湖州师范学院学报 第３６卷 式中：ｚ川一 ＋∑ｙ 是 为当前计算时刻，中间参数ｏ，ｉ， 可由ＲＫＦ算法系数表给出，此时状态向量可 ＝１ 表示为： ｚ川＝ｚ ＋∑ｙ 愚 ． ＝１ 如 如加 ０ ｍ加 ６ 上述算法中，若定义一个误差向量ｅ 一∑（ｙ 一ｙ ） 嘲，则可由￡ 大小变换步长，从而转化为自适 ｉ一１ 应变步长算法，使程序整体更易控制． ２．２ Ｌｏｒｅｎｚ模型的数值算法 Ｌｏｒｅｎｚ模型是描述混沌现象的一例著名方程，Ｌｏｒｅｎｚ系统无论从数学还是物理的角度都值得详细地 研究． 一 一 一 ．＼ 、，＼ 其中：Ｉ９，』Ｄ， 为定值 ；ｚ （０）一ｚ ２（０）一ｚ。（Ｏ）一￡为初始值．对于￡＝＝＝１０叫。，求方程组的数值解． ＋＋ 由于该方程为非线性微分方程，不存在解析解，可按以下步骤完成基于ＭＡＴＬＡＢ的数值算法程序编写： ｝；ｌ ㈤ （ ３ ＋ 首先，用Ｌｏｒｅｎｚｅｐ．ｍ函数来描述系统的动态模型，格式为： 加 ｍ Ｏ 一 加加 一彻 ＞＞ｆ一夺（ｔ，ｘ）Ｅｔ￣＊ｘ（１）＋ｘ（２）＊ｘ（３）； ； ；＊］．Ｏ 其次，调用数值解函数ｏｄｅ４５（）对Ｌｏｒｅｎｚｅｐ（）函数描述的系统进行数值求解，并将结果用图形显示．／．： 一＞＞ｔｆｉｎａ１—１００；ｘ０一［Ｏ；０；ｅ～１０］；Ｉｔ，ｘ］￣ｏｄｅ４５（ｆ，［Ｏ，ｔ＿ｆｉｎａ１］，ｘ０）；ｐｌｏｔ（ｔ，ｘ）， —ｆｉｇｕｒｅ；ｐｌｏｔ３（ｘ（：，１）；ｄ（：，２）；ｘ（：，３））； ＼， ａｘｉｓ（［ｍｌ ｍ２ ｍ ｍ４ ｍ５ ｍ６］）； 其中：ｔ ｆｉｎａｌ为设定的仿真终止时间；ｘ０为初始状态．两个命令分别绘制出系统各个状态的时间关系的二 维曲线图（见图３）、三个状态的相空间曲线图（见图４）． 图４相空间三维图 ４ Ｔｈｒｅｅ—ｄｉｍｅｎｓｉｏｎａｌ Ｐｈａｓｅ ｓｐａｃｅ ｄｉａｇｒａｍ 上面程序中，对于观察相空间轨迹还可采用ｃｏｍｅｔ３（）函数绘制动画式轨迹，只需对最后一个语句进行适 当改写即可．要建立一个基于ＭＡＴＬＡＢ微分方程组描述函数，必须调用ｏｄｅ４５（），所以编写ＭＡＴＬＡＢ函数描 述微分方程是初值问题数值解的关键，对一阶非线性微分方程尤为重要． ３ 结语 微分方程数值解的检验是数值分析理论的重要环节．若仿真算法和控制参数选择不当，则可能得出不可 信的结果，甚至是错误的结果．对数值解的验证可通过两种方式进行：一是修改仿真控制参数，如可以接受的 第４期 薛亚宏：两类基于ＭＡＴＬＡＢ的非线性微分方程数值解的算法研究 ２３ 误差限，将选项设置成一个更小的值，观察所得的结果，看看是否与第一次运行时结果一致，如果不存在不能 接受的差异，则应考虑再进一步减小误差限；二是选择微分方程求解算法检验运算结果的正确性． 科学或工程问题的求解和模拟最终往往都要归结为对数学模型的求解和对实验数据的拟合，而非线 性微分方程作为一种重要的数学终端模型，本身只能用数值求解．本文通过对四阶定步长Ｒｕｇｅ—Ｋｕｔｔａ算 法和Ｌｏｒｅｎｚ模型算法的系统分析，建立基于非线性微分方程的变步长自适应算法及函数显式描述下 ＭＡＴＬＡＢ程序结构，对步长、误差、数值解验证等关键要素进行分析，并得出相关结论，为以微分方程为 数学模型的工程问题的解决提供了一种有效的数值求解方法，对偏微分方程数值解的进一步讨论有重要 意义． 参考文献： ［１］商妮娜，秦惠增．非线性常微分方程数值解的渐近表示以及应用［Ｊ］．山东理工大学学报（自然科学版），２０１０，３３（５）：３０—３５ ［２］张春梅．一类二阶非线性微分方程边值问题的有效解法［Ｊ３．新疆大学学报（自然科学版），２００５，２２（４）：４０３—４０９． ［３］薛定宇，陈阳泉．基于的系统仿真技术与应用ＥＭ］．北京：清华大学出版社，２００２． ［４］邵新平．非线性微分方程的样条函数求解方法Ｉ－Ｊ］．浙江大学学报（自然科学版），２０１０（３）：１０—１２． ［５］蔡建平，王桂芳．非线性微分方程渐近解的阶的数值验证［Ｊ］．太原理工大学学报，２０１１，２０（５）：１５—２１． ［６］张化光，王智良，黄伟．混沌系统的控制理论［Ｍ］．沈阳：东北大学出版社，２００３． Ｔｈｅ Ａｌｇｏｒｉｔｈｍｓ Ａｂｏｕｔ ｔｈｅ Ｔｗｏ Ｔｙｐｅｓ ｏｆ Ｎｕｍｅｒｉｃａｌ Ｓｏｌｕｔｉｏｎ ｏｆ Ｎｏｎｌｉｎｅａｒ Ｄｉｆｆｅｒｅｎｔｉａｌ Ｅｑｕａｔｉｏｎｓ Ｂａｓｅｄ ｏｎ ＭＡＴＬＡＢ ＸＵＥ Ｙａｈｏｎｇ （Ｓｃｈｏｏｌ ｏｆ Ｔｅｌｅｃｏｍｍｕｎｉｃａｔｉｏｎｓ，Ｇａｎｓｕ Ｉｎｄｕｓｔｒｉａｌ Ｖｏｃａｔｉｏｎａｌ＆Ｔｅｃｈｎｏｌｏｇｉｃａｌ Ｃｏｌｌｅｇｅ，Ｔｉａｎｓｈｕｉ ７４１０２５，Ｃｈｉｎａ） Ａｂｓｔｒａｃｔ：Ｂａｓｅｄ ｏｎ ｔｈｅ ｍａｔｈｅｍａｔｉｃａｌ ａｌｇｏｒｉｔｈｍ ａｎａｌｙｓｉｓ ｆｏｒ ｔｈｅ ｔｗｏ ｋｉｎｄｓ ｏｆ ｎｏｎｌｉｎｅａｒ ｄｉｆｆｅｒｅｎｔｉａｌ ｅ— ｑｕａｔｉｏｎ ｎｕｍｅｒｉｃａｌ ｓｏｌｕｔｉｏｎｓ ｉｎ ｅｎｇｉｎｅｅｒｉｎｇ ｄｙｎａｍｉｃ ｃｏｎｔｒｏｌ ａｎｄ ｔｈｅ ｉｍｐｏｒｔａｎｔ ａｐｐｌｉｃａｔｉｏｎｓ ｉｎ ｃｏｍｐｕｔｅｒ ｓｉｍｕｌａｔｉｏｎ，ｔｈｅ ｅｓｓａｙ ｅｓｔａｂｌｉｓｈｅｓ ａ ｆｏｕｒｔｈ—ｏｒｄｅｒ ｆｉｘｅｄ ｓｔｅｐ ｌｅｎｇｔｈ Ｒｕｎｇｅ—Ｋｕｔｔａ ａｎｄ ｔｈｅ ａｌｇｏｒｉｔｈｍ ｓｔｒｕｃｔｕｒｅ ｏｆ Ｌｏｒｅｎｚ ｍｏｄｅｌ ｎｕｍｅｒｉｃａｌ ｓｏｌｕｔｉｏｎ ＭＡＴＬＡＢ，ｄｉｓｃｕｓｓｅｓ ｔｈｅ ｅｒｒｏｒ ｃｏｎｔｒｏｌ ｖａｒｉａｂｌｅ ｓｔｅｐ ｓｉｚｅ ｃａｓｅ，ｄｒａｗｓ ａ ｎｕｍｅｒｉｃａｌ ｓｏｌｕｔｉｏｎ ｏｆ ｔｈｅ Ｌｏｒｅｎｚ ｓｙｓｔｅｍ ｕｎｄｅｒ ｔｈｅ ｔｗｏ—ｄｉｍｅｎｓｉｏｎａｌ ａｎｄ ｔｈｒｅｅ—ｄｉｍｅｎ— ｓｉｏｎａｌ ｇｒａｐｈｉｃｓ ｂａｓｅｄ ｏｎ ｔｈｅ ＭＡＴＬＡＢ，ａｎｄ ｆｉｎａｌｌｙ ｐｏｉｎｔｓ ｏｕｔ ｔｈｅ ｎｕｍｅｒｉｃａｌ ｓｏｌｕｔｉｏｎ ｉｎ ｔｅｓｔｉｎｇ ｔｈｅ ｂａｓｉｃ ｔｒａｉｎ ｏｆ ｔｈｏｕｇｈｔ ｗｉｔｈｉｎ ａｃｃｅｐｔａｂｌｅ ｅｒｒｏｒ ｌｉｍｉｔｓ． Ｋｅｙ ｗｏｒｄｓ：ＭＡＴＬＡＢ；ｔｈｅ ｎｏｎｌｉｎｅａｒ ｄｉｆｆｅｒｅｎｔｉａｌ ｅｑｕａｔｉｏｎ；ｎｕｍｅｒｉｃａｌ ｓｏｌｕｔｉｏｎ；ａｌｇｏｒｉｔｈｍ；ｅｒｒｏｒ ＭＳＣ ２０００：６５Ｌ１４ ［责任编辑高俊娥］