前　言 
PREFACE
科学运算问题是每个理工科学生和科技工作者在课程学习、科学研究与工程实践中常常遇到的问题，不容回避。对于非纯数学专业的学生和研究者而言，从底层全面学习相关数学问题的求解方法并非一件简单的事情，也不易得出复杂问题的解。所以，利用当前最先进的计算机工具，高效、准确、创造性地求解科学运算问题是一种行之有效的方法，尤其能够满足理工科人士的需求。
作者曾试图在同一部著作中叙述各个数学分支典型问题的直接求解方法，由清华大学出版社出版了《高等应用数学问题的 MATLAB求解》。该书从 2004年出版之后多次重印再版，并于 2018年出版了第 4版，还配套发布了全新的 MOOC课程 ①，一直受到广泛的关注与欢迎。首次 MOOC开课的选课人数接近 14 000人，教材内容也被数万篇期刊文章和学位论文引用。
从作者首次使用 MATLAB语言算起，已经有 30余年了，通过相关领域的研究、思考与一线教学实践，积累了大量的实践经验资料。这些不可能在一部著作中全部介绍，所以与清华大学出版社策划与编写了这套“薛定宇教授大讲堂”系列著作，系统深入地介绍基于 MATLAB语言与工具的科学运算问题的求解方法。
本系列著作不是原来版本的简单改版，通过 10余年的经验和资料积累，全面贯穿 “再认识 ”的思想写作此书，深度融合科学运算数学知识与基于 MATLAB的直接求解方法与技巧，力图更好地诠释计算机工具在每个数学分支的作用，帮助读者以不同的思维与视角了解工程数学问题的求解方法，创造性地得出问题的解。
本系列著作卷 I可以作为学习 MATLAB入门知识的教材与参考书，也为读者深入学习与熟练掌握 MATLAB语言编程技巧，深度理解科学运算领域 MATLAB的应用奠定一个坚实的基础。后续每一卷试图对应一个数学专题或一门数学课程进行展开。全套系列著作的写作贯穿“计算思维”的思想，深度探讨该数学专题的问题求解方法。本系列著作既适合于学完相应的数学课程之后，深入学习利用计算 
. MOOC网址：https://www.icourse163.org/learn/NEU-1002660001。
·ii·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
机工具的科学运算问题求解方法与技巧，也可作为相应数学课程同步学习的伴侣，在学习相应课程理论知识的同时，侧重学习基于计算机的数学问题求解方法，从另一个角度观察、审视数学课程所学的内容，扩大知识面，更好地学习、理解和实践相应的数学课程。
本书是系列著作的卷 V。本书系统介绍基于 MATLAB的微分方程求解方法。首先介绍微分方程的解析解方法，然后介绍各种微分方程的数值求解方法。对微分方程的初值问题，有针对性地给出一般常微分方程、特殊微分方程以及延迟微分方程和分数阶微分方程的数值求解方法，探讨微分方程的某些性质与行为，并介绍基于框图的初值问题求解方法。本书还介绍微分方程边值问题与偏微分方程的数值求解方法。
值此系列著作付梓之际，衷心感谢相濡以沫的妻子杨军教授，她数十年如一日的无私关怀是我坚持研究、教学与写作的巨大动力。
薛定宇 2019年 10月


目　录 
CONTENTS
第 1章微分方程简介 ······································· 1 
1.1微分方程建模简介···································· 1 
1.1.1电路的建模 ··································· 1 
1.1.2力学问题的建模································ 3 
1.1.3社会系统的建模································ 3 
1.2微分方程发展简史···································· 5 
1.3本书主要内容 ······································· 8本章习题 ················································ 9第 2章常微分方程的解析解 ·································· 10 
2.1一阶微分方程的解析解 ································ 10 
2.1.1可由简单积分求解的微分方程····················· 11 
2.1.2齐次线性方程·································· 11 
2.1.3非齐次线性方程································ 12 
2.1.4可分离变量的非线性微分方程····················· 13 
2.2特殊函数与二阶线性微分方程 ·························· 14 
2.2.1 Gamma函数 ·································· 15 
2.2.2超几何函数 ··································· 16 
2.2.3 Bessel微分方程 ································ 17 
2.2.4 Legendre微分方程与 Legendre函数 ················· 19 
2.2.5 Airy函数 ····································· 20 
2.3常系数线性微分方程的求解 ···························· 21 

2.3.1线性常系数微分方程解析解的数学描述 ············· 21 
2.3.2基于 Laplace变换的求解方法······················ 22 
2.3.3非齐次微分方程的求解 ·························· 24 
2.3.4非零初值的微分方程求解 ························ 25 
·iv·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 
2.4一般微分方程的解析解 ································ 27 
2.4.1简单微分方程的解析解 ·························· 27 
2.4.2常系数高阶线性微分方程的解析解 ················· 29 
2.4.3线性时变微分方程的解析解······················· 31 
2.4.4线性时变微分方程组的求解······················· 32 
2.4.5边值问题的计算机求解 ·························· 33 
2.5线性矩阵微分方程的求解 ······························ 34 
2.5.1线性状态空间方程的解析解······················· 35 
2.5.2状态方程的直接求解 ···························· 36 
2.5.3 Sylvester微分方程的求解 ························ 37 

2.5.4基于 Kronecker乘积的 Sylvester微分方程直接求解 ····· 38 
2.6特殊非线性微分方程的解析解 ·························· 39 
2.6.1可解的非线性微分方程 ·························· 39 
2.6.2解析解不存在的非线性微分方程 ··················· 41本章习题 ················································ 41第 3章微分方程的初值问题 ·································· 45 
3.1一阶显式微分方程组的初值问题························· 45 
3.1.1初值问题的数学形式 ···························· 45 
3.1.2初值问题解的存在性与唯一性····················· 46 
3.2定步长数值算法与实现 ································ 46 
3.2.1 Euler算法 ···································· 47 
3.2.2二阶 Runge–Kutta算法 ·························· 50 
3.2.3四阶 Runge–Kutta算法 ·························· 51 
3.2.4 Gill算法 ······································ 52 
3.2.5 m阶 Runge–Kutta算法 ·························· 53 
3.2.6定步长多步算法与实现 ·························· 56 
3.3变步长数值算法与实现 ································ 58 
3.3.1提高求解效率的措施 ···························· 58 
3.3.2变步长方法简介································ 59 
3.3.3四级五阶 Runge–Kutta变步长算法 ················· 60 
3.3.4基于 MATLAB的微分方程求解函数 ················ 61 

3.3.5基于 MATLAB的带有附加参数的微分方程求解 ······· 65 
3.3.6避免附加参数的方法 ···························· 67 
3.4微分方程数值解的验证 ································ 68 
3.4.1计算结果的验证································ 68 
3.4.2中间计算结果的动态处理 ························ 70 
3.4.3更高精度的数值计算函数 ························ 71 
3.4.4计算步长与定步长显示 ·························· 72 
3.4.5高阶非线性微分方程的求解实例 ··················· 74本章习题 ················································ 75
第 4章微分方程的标准型变换 ································ 80 
4.1单个高阶常微分方程变换方法 ·························· 80 
4.1.1高阶显式微分方程的变换 ························ 81 
4.1.2时变微分方程的求解方法 ························ 84 
4.1.3微分方程的奇点································ 85 
4.1.4含有常数参数的状态增广方法····················· 87 
4.2复杂高阶微分方程的变换与求解························· 88 
4.2.1含有最高阶导数二次方的微分方程 ················· 88 

4.2.2含有最高阶导数奇数次方的微分方程 ··············· 90 
4.2.3含有最高阶导数的非线性运算····················· 91 
4.3高阶常微分方程组的变换 ······························ 92 
4.3.1简单的显式微分方程组 ·························· 92 
4.3.2定步长算法的局限性 ···························· 98 
4.3.3简单的隐式微分方程组 ·························· 100 
4.3.4更复杂的非线性方程组 ·························· 102 
4.4矩阵型微分方程的变换 ································ 104 
4.4.1矩阵型微分方程的变换与求解····················· 104 
4.4.2 Sylvester微分方程 ······························ 106 
4.4.3 Riccati微分方程 ······························· 107 
4.5一类 Volterra积分微分方程的变换 ······················· 109本章习题 ················································ 112
第 5章特殊微分方程 ······································· 116 
5.1刚性微分方程 ······································· 116 
5.1.1线性微分方程的时间常数 ························ 117 
5.1.2刚性现象 ····································· 117 
5.1.3刚性微分方程的直接求解 ························ 119 
5.1.4微分方程刚性的检测 ···························· 122 
5.1.5刚性微分方程的定步长求解······················· 126 
·vi·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 
5.2隐式微分方程 ······································· 127 
5.2.1隐式微分方程的一般数学描述····················· 127 
5.2.2隐式微分方程相容初值的变换····················· 129 
5.2.3隐式微分方程的直接求解 ························ 131 
5.2.4多解隐式微分方程的求解 ························ 134 
5.3微分代数方程 ······································· 135 
5.3.1微分代数方程的一般形式 ························ 135 
5.3.2微分代数方程的指数类型 ························ 136 
5.3.3半显式微分代数方程的直接求解 ··················· 136 

5.3.4微分代数方程直接求解方法的局限性 ··············· 139 
5.3.5一般微分代数方程的隐式微分方程求解 ············· 140 
5.3.6微分代数方程的指数降型方法····················· 145 
5.4切换微分方程 ······································· 147 
5.4.1线性切换微分方程 ······························ 147 
5.4.2过零点检测与事件设置 ·························· 148 
5.4.3非线性切换微分方程 ···························· 151 
5.4.4不连续微分方程································ 152 
5.5线性随机微分方程···································· 154 
5.5.1线性随机微分方程的传递函数····················· 154 
5.5.2连续随机系统仿真的误区 ························ 155 
5.5.3随机线性系统的离散化 ·························· 156本章习题 ················································ 160
第 6章延迟微分方程 ······································· 164 
6.1带有延迟常数的延迟微分方程数值解 ····················· 164 
6.1.1从普通微分方程到延迟微分方程 ··················· 164 
6.1.2零历史函数的延迟微分方程求解 ··················· 166 
6.1.3非零历史函数的延迟微分方程····················· 170 
6.2变延迟的微分方程···································· 172 
6.2.1变延迟的微分方程模型 ·························· 172 
6.2.2基于时间延迟的延迟微分方程····················· 173 
6.2.3基于状态延迟的微分方程 ························ 176 
6.2.4带有广义延迟的延迟微分方程····················· 177 
6.3中立型延迟微分方程的求解 ···························· 179 
6.3.1中立型延迟微分方程 ···························· 179 
6.3.2变延迟中立型微分方程 ·························· 182 
6.4带有延迟的 Volterra积分微分方程 ······················· 183本章习题 ················································ 184
第 7章微分方程的性质与行为 ································ 187 
7.1微分方程的稳定性···································· 187 
7.1.1常系数线性微分方程的稳定性····················· 187 
7.1.2 Routh–Hurwitz稳定性判据 ······················· 189 

7.1.3 Lyapunov函数与 Lyapunov稳定性 ················· 192 
7.1.4时变微分方程的自治化 ·························· 193 
7.1.5一般非线性系统的稳定性判定····················· 193 

7.1.6基于数值仿真的复杂系统稳定性判定 ··············· 195 
7.2微分方程的特殊行为 ·································· 197 
7.2.1极限环 ······································· 198 
7.2.2周期解 ······································· 201 
7.2.3混沌与吸引子·································· 204 
7.2.4 Poincaré映射 ·································· 208 
7.3微分方程的线性化近似 ································ 210 
7.3.1平衡点 ······································· 210 
7.3.2非线性微分方程的线性化 ························ 213 
7.3.3平衡点的性态·································· 216 
7.4微分方程的分岔 ····································· 217本章习题 ················································ 218
第 8章分数阶微分方程······································ 219 
8.1分数阶微积分的定义与数值计算························· 220 
8.1.1分数阶微积分的定义 ···························· 220 

8.1.2不同分数阶微积分定义的关系与性质 ··············· 221 
8.1.3 Grünwald–Letnikov定义的数值计算 ················ 222 
8.1.4 Caputo微积分定义的数值计算 ···················· 223 
8.2同元次线性分数阶微分方程的解析解 ····················· 224 
8.2.1 Mittag-Leffler函数 ······························ 224 
8.2.2同元次线性分数阶微分方程······················· 225 
8.2.3一个重要的 Laplace变换公式······················ 226 
8.2.4基于部分分式展开的解析解方法 ··················· 227 
8.3常系数线性分数阶微分方程的数值求解 ··················· 231 
8.3.1线性方程的闭式解法 ···························· 231 

8.3.2 Riemann–Liouville微分方程 ······················ 233 
·viii·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 
8.3.3 Caputo微分方程 ······························· 235 
8.3.4等效初值的计算································ 237 
8.3.5微分方程的高精度算法 ·························· 239 
8.4非线性分数阶微分方程的求解 ·························· 242 
8.4.1预估方程 ····································· 243 
8.4.2校正求解方法·································· 246 
8.4.3隐式 Caputo微分方程的高精度矩阵算法············· 247本章习题 ················································ 249
第 9章常微分方程的框图求解 ································ 251 
9.1 Simulink必备知识 ···································· 252 
9.1.1 Simulink简介 ·································· 252 
9.1.2 Simulink相关模块 ······························ 252 
9.2微分方程的框图建模思想 ······························ 254 
9.2.1积分器链与关键信号生成 ························ 254 
9.2.2微分方程的框图描述方法 ························ 255 
9.2.3微分方程的求解································ 257 
9.2.4算法与参数设定································ 258 
9.3微分方程建模举例···································· 260 
9.3.1一般微分方程组································ 260 
9.3.2微分代数方程·································· 263 
9.3.3切换微分方程·································· 265 
9.3.4不连续微分方程································ 267 
9.3.5延迟微分方程·································· 267 
9.3.6非零历史函数的延迟微分方程····················· 269 
9.3.7随机微分方程·································· 271 
9.4分数阶微分方程的 Simulink求解························· 272 
9.4.1分数阶算子的模块逼近 ·························· 273 

9.4.2 Riemann–Liouville分数阶微分方程的建模与求解 ······ 274 
9.4.3 Caputo导数的模块计算·························· 276 
9.4.4 Caputo分数阶微分方程的建模与求解··············· 277 
9.4.5分数阶延迟微分方程 ···························· 279
本章习题 ················································ 280第 10章微分方程的边值问题 ································· 283 
10.1微分方程标准边值问题 ······························· 283 10.2二阶微分方程两点边值问题的打靶求解 ·················· 284 
10.2.1线性时变方程边值问题的打靶算法 ················ 285 
10.2.2线性微分方程的有限差分算法 ···················· 287 
10.2.3非线性方程边值问题的打靶算法 ·················· 289 
10.3高阶微分方程两点边值问题 ··························· 293 
10.3.1 MATLAB的直接求解函数······················· 293 
10.3.2简单边值问题的求解 ··························· 294 
10.3.3复杂边值条件的描述与求解······················ 298 
10.3.4带有待定参数的边值问题 ······················· 298 
10.3.5半无穷区间的边值问题 ························· 301 
10.3.6带有浮动边值的多解微分方程 ···················· 302 
10.3.7积分微分方程的边值问题 ······················· 303 

10.4基于最优化技术的微分方程边值问题求解 ················ 304 
10.4.1简单边值问题的最优化求解······················ 304 
10.4.2隐式微分方程的边值问题 ······················· 305 
10.4.3延迟微分方程的边值问题 ······················· 308 
10.4.4多点已知值的微分方程问题······················ 309 
10.4.5浮动边值问题的重新求解 ······················· 311 
10.4.6基于框图的边值问题求解方法 ···················· 312 
10.4.7分数阶微分方程的边值问题······················ 313本章习题 ················································ 314
第 11章偏微分方程入门····································· 316 
11.1扩散方程的数值求解 ································· 317 
11.1.1一维扩散方程的数学形式与解析解 ················ 317 
11.1.2扩散方程的离散化方法 ························· 318 
11.1.3非齐次扩散方程 ······························· 321 
11.1.4高维扩散方程的数学形式 ······················· 323 
11.2几种特殊形式的偏微分方程 ··························· 323 
11.2.1偏微分方程的分类 ····························· 323 
11.2.2特征值型偏微分方程 ··························· 325 
11.2.3边界条件的分类 ······························· 325 
11.3典型二维偏微分方程求解界面·························· 326 
11.3.1偏微分方程求解程序概述 ······················· 326 
11.3.2偏微分方程几何区域绘制 ······················· 327 
11.3.3偏微分方程边界条件描述 ······················· 328 
·x·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 
11.3.4偏微分方程求解举例 ··························· 329 
11.3.5解的其他显示方法 ····························· 330 
11.3.6函数参数的偏微分方程求解······················ 332 
11.4一般偏微分方程的求解 ······························· 333 
11.4.1创建空白的偏微分方程对象模型 ·················· 333 
11.4.2几何区域的语句描述 ··························· 333 
11.4.3边界条件与初始条件描述 ······················· 336 
11.4.4偏微分方程的描述 ····························· 337 
11.4.5偏微分方程的数值求解 ························· 338本章习题 ················································ 343参考文献··················································· 345 
MATLAB函数名索引 ········································· 350术语索引··················································· 354 
微分方程简介  
第 1章 
MATLAB微分方程求解 


方程是含有未知变量的等式。通常，方程分为代数方程与微分方程。本丛书卷 III和卷 IV介绍过代数方程及其求解方法。代数方程主要描述未知变量之间的静态关系，微分方程则是动态系统的数学基础。
本丛书卷 II全面介绍了函数与微积分的理论与计算方法，常微分方程主要描述未知变量之间的动态关系。一般情况下，除了描述未知变量之间的关系，常微分方程还描述未知变量与时间 t之间的关系，同时描述未知变量当前时刻的值与其过去的值之间的关系。因此，微分方程的描述与求解通常比代数方程更复杂。
本章 1.1节以电路、力学问题与社会系统为例，介绍从现象建立微分方程模型的方法，1.2节介绍微分方程的发展简史，1.3节给出本书内容的概略介绍。 

1.1微分方程建模简介
在科学技术研究中，甚至在日常生活与社会科学领域，通常可以将某一种现象或某一个系统用微分方程描述，本节将举几个简单的例子演示电路的建模、力学系统的建模与社会系统建模等。 
1.1.1电路的建模
在电路理论中，由欧姆定律可知，电阻中通过的电流与其两端的电压之间的关系是静态的，即 u(t)= Ri(t)，由当前时刻的 i(t)可以直接计算出 u(t)。对一个纯电阻电路而言，整个电路可以由静态方程表示。静态方程可描述整个电路中任意一个支路的电流与任意元件两端的电压的关系。
对于另外两种常用的电气元件电容与电感，则需要引入微分方程描述其瞬时电流与瞬时电压之间的关系。
（1）如果电路中含有电容元件 C，则可以由下面的微分方程描述其瞬时电压 u(t)与瞬时电流 i(t)之间的关系，即 i(t)= R dud(tt) （1-1-1）

·2
·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解


（2）如果电路中含有电感元件 L，则可以由下面的微分方程描述其瞬时电压 
u(t)与瞬时电流 i(t)之间的关系，即 u(t)= C did(tt) （1-1-2）
下面将通过例子演示简单电路模型的微分方程建模方法。
电路理论基本知识可知，电感元件 两端的电压与回路电流 之间满足如下的微分Li()t. .. .. 
　　例 1-1考虑图 1-1给出的简单电路，试建立该电路的数学模型。 

Li(t) R+ + u(t) C 
uc(t) 
. . 
图 1-1电阻、电容与电感的串联电路
　　解该电路主要描述给定电压 u(t)与电容 C两端电压 uc(t)之间的动态关系，前者称为电路的输入信号，后者称为电路的输出信号。回路电流 i(t)为电路的中间信号。由
方程： ul(t)= L did(tt) （1-1-3）
由著名的 Kirchhooff定律可以直接写出输入信号与输出信号之间的方程，即 
u(t)= Ri(t)+ L did(tt)+ uc(t)（1-1-4）
另外，由电路理论可知，电容两端的电压 uc(t)与回路电流 i(t)满足下面的微分方
程模型。 i(t)= C dudct(t) （1-1-5）
将方程（1-1-5）代入方程（1-1-4），则得出到整个电路的微分方程模型为 
d2u(t)= RC dudct (t)+ LC dutc2 (t)+ uc(t)（1-1-6）
另外，如果选择 uc(t)与 i(t)为电路的状态，则可以写出下面的状态方程模型（有关状态方程模型的概念将在后面介绍）。 
dudct (t)= C 1 i(t) （1-1-7）
di(t)1 R 1 
= u(t) . i(t)+ uc(t)
dtL LL 
可以看出，由简单的微积分运算可以建立电路的微分方程模型。如果有多个回路，则有可能建立高阶微分方程。如果电路中存在非线性元件（例如二极管、三极管等电子元件），则电路模型可能为非线性微分方程模型。可以看出，微分方程模型是描述电路中动态现象的最常用的也是不可或缺的数学模型。
这里描述的关系是集中参数（ lumped parameter）电路经常使用的关系。在真实的电路中，由于电感存在漏电阻等因素，单纯使用常微分方程描述系统模型不够准确，需要引入偏微分方程描述，这时采用的模型为分布参数（ distributed parameter）模型。 

1.1.2力学问题的建模
先考虑中学课程中学习的 Newton第二定律，即 F (t)= Ma(t)（1-1-8）其中，物体的质量 M是常数。为了更好地刻画变量的动态性能，这里有意将外力与加速度都写成时间 t的函数，可以理解为瞬时施加的外力 F (t)和物体的瞬时加速度 a(t)。在实际系统中，物体的加速度一般是不能通过直接测量得到的，可以测量的相关物理量是物体的瞬时位置信息 x.(t)，则物体的瞬时位移可定义为 x(t)= x.(t) . x.(t0)。位移与外力之间的关系是什么呢？因为加速度是位移的二阶导数，因此由 Newton第二定律可以直接写出系统的微分方程模型为 F (t)= Mx ′′ (t)（1-1-9）
下面给出一个简单的例子介绍力学系统的微分方程建模问题。
　　例 1-2考虑图 1-2给出的简单力学系统示意图，如果以外力 F (t)拖动物体，那么物体的位移 x(t)与外力 F (t)之间的数学模型是什么？


　　解 Newton第二定律需要考虑的是物体受力的总和。假设外力 F (t)是正方向的力，则反方向的力包括弹簧的力和阻尼器的力。由 Hooke定律可知，弹簧的力与位移成正比，记作 .Kx(t)；阻尼器的力与重物的速度成正比，记作 .R dx(t)/dt。这样，由 Newton第二定律建立该力学系统的微分方程模型为 
F (t) . Kx(t) . R dx(t) dt  = M d2x(t) dt2 （1-1-10）  
1.1.3 社会系统的建模 

在自然科学与工程领域，微分方程模型是不可或缺的，因为所有连续动态过程都需要通过微分方程进行描述。另外，在社会科学等领域，也常常需要使用微分方
·4·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
程描述某些现象与规律。下面给出几个社会领域的微分方程建模实例。
　　例 1-3人口数量模型。 Malthus人口模型是英国牧师、经济学家 Thomas Robert Malthus（1766.1834）分析了以往英国的人口数据之后提出的。如果 t时刻的世界人口数用 x(t)表示，且已知人口的自然增长率常数为 r，则 Malthus的人口模型可以用微分
方程表示为 x ′ (t)= rx(t),且已知 x(0) = x0
该微分方程的解为 x(t)= x0er(t.t0)，即人口是按指数规律增长的。如果 t很大，则人口数将趋于无穷大，这是不现实的，因为人口达到一定数量时，空间与资源等约束条件将抑制人口数量的无限增长。
比利时数学生物学家 Pierre-Fran.ois Verhulst构造了一个函数 r(t)，后来被称为 Logistic函数（又因其形状被称为 Sigmoid函数或 S型函数）。这样，人口数量模型变成
了下面的时变微分方程。 x ′ (t)= r(t)x(t),且已知 x(0) = x0
　　例 1-4 Richardson军备竞赛模型 [1]。Richardson军备竞赛模型是英国数学家、物理学家 Lewis Fry Richardson（1881.1953）提出的 [2]。Richardson军备竞赛模型有三个前提。假设有两个国家 X和 Y，如果一个国家发现另一个在购买武器上花费了大量的财富，则它自己也会花更多的钱去购买武器；但是，军费开支本身是社会的经济负担，更高的军费支出水平将抑制未来支出的增长，这是第二个前提；第三个前提是，受本国文化或国家领导者意志的影响，增加军费会受到鼓励或压制。在这三个前提下， Richardson提出了描述军备竞赛的微分方程模型。 
{ 
x ′ (t)= ay(t) . mx(t)+ g （1-1-11）
y ′ (t)= bx(t) . ny(t)+ h
　　例 1-5 Lotka–Volterra捕食者与猎物模型。假设在某个封闭的区域内，只有狐狸和兔子两种动物，初始时间各有一定的数量。由于狐狸多了，会捕食更多的兔子，兔子数量就会减少，而兔子减少，狐狸就会因饿死而减员；狐狸减少，由于缺乏天敌，兔子数量将增加，这种现象会周而复始。假设 t时刻兔子与狐狸数量分别为 x(t)和 y(t)，则可以由下面的微分方程描述这两种生物的数量变化。 
{ 
x ′ (t)=叫
x(t) . Δx(t)y(t) y ′ (t)= αx(t)y(t) . Ψy(t)
其中，叫、Δ、α和 Ψ为正实数。该模型是由美国数学家 Alfred James Lotka（1880.1949）首先提出的，后来意大利数学家 Vito Volterra（1860.1940）也给出了类似的模型。该模型也常应用于经济学等领域。
　　例 1-6流行病传染模型。英国流行病学家 William Ogilvy Kermack（1898. 1970）和英国数学家 Anderson Gray McKendrick（1876.1943）在 1927年提出了流行病传染
模型 [3]
，
。
后来其模型中的变量名被替换并写成如下的形式 [4]
S ′ (t)= .ΔI(t)S(t),S(0) = S0 
I ′ (t)= ΔI(t)S(t) . ΨI(t),I(0) = I0 

. .. 
..
R ′ (t)= ΨI(t),R(0) = R0


. .. .. 1.2微分方程发展简史
所以，该模型又称 SIR模型。其中，S(t)、I(t)和 R(t)分别是易感（susceptible，表示未患病者）人群的人数、感染者（infective）的人数、移除者（removal，指隔离或死亡）的人数，且 S(t)+ I(t)+ R(t)= N(t)，N(t)为人口总数。为方便处理，对这些数量作了归一化处理，使得 S(t)+ I(t)+ R(t)=1。在模型参数中，Δ是单位时间内被感染者接触的次数，又称为感染率；Ψ为移除率，表示隔离或死亡的人数与感染的人数的比例。
如果该病症有 ξ天的潜伏期和感染期，则 SIR模型可改写成下面的延迟微分代数
[4]方程。 ′ .S ()= ΔI()S()ttt√ 
I ′ (t)= ΔI(t)S(t) . ΔI(t . ξ )S(t . ξ) 
S(t)+ I(t)+ R(t)=1
从以上例子可以看出，在社会科学研究中，若掌握了微分方程这样强有力的数学工具，可以在研究中引入新的视角，探讨很多难以研究的问题，得出定量的结果。 
英国科学家 Isaac Newton（1643.1727，图 1-3（a）
最早研究了微分方程，他采用的方法是幂级数法。例如，Newton在 1671年的著作“Methodus Fluxionum et Serierum Infinitarum”（流数法与无穷级数）中，研究了微分方程 [5] 
y ′ (x)=1 . 3x + y(x)+ x 2 + xy(x)（1-2-1）
利用幂级数的方法，将 y(x)设置成一个系数未知的幂级数，代入方程并令等号两端 x同次幂的系数相等，再求解代数方程，得出的结果为 
111 1
66
y(x)= x . xx + x 3 . x + x 5 . x 
3 6 30 45 
其中，xx是 Newton原文使用的记号，其实该项即 x2。
当然，从现代数学对微分方程的认知看，这并不是一个完美的解。后面将介绍该问题的解，并介绍借助于计算机的求解方法。
德国数学家 Gottfried Wilhelm Leibniz（1646.1716，图 1-3（b）
在 1693年求解了一个几何学上的微分方程 [5] 
y ′ (x)= . a2 y(. x) y(x)（1-2-2）
他利用积分方法得出了该微分方程的解为 
√ 
√
∫
a 
2
a2 . y
√ 

2
a2 . y
2 . a ln a .
dy = . a2 . y
x = 

yy
y 
·6·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 
Leibniz还提出了可分离变量的微分方程求解方法，可以求解一大类微分方程。
其实，早在 Newton或 Leibniz建立微积分学说之前，意大利物理学家、天文学家 Galileo Galilei（1564.1642，图 1-3（c）
就研究了一些通常需要微分方程研究的问题，但当时尚未建立微积分理论，所以他采用的方法自然不是微分方程，而是传统的代数与几何方法。

（a）Isaac Newton（b）Gottfried Leibniz（c）Galileo Galilei图 1-3 Newton、Leibniz、Galileo的画像注：图像均来源于维基百科
法国数学家 Alexis Claude Clairaut（1713.1765，图 1-4（a）
对一类特别的应用问题感兴趣，研究了最早的隐式微分方程 [5]。 
y(x) . xy ′ (x)+ f(y ′ (x)) = 0（1-2-3）瑞士数学家 Jacob Bernoulli（又名 James或 Jacques，1654.1705，图 1-4（b）
和 John Bernoulli（又名 Jean，1667.1748）研究了微分方程问题，并提出可以“解析”求解的一类非线性微分方程，称为 Bernoulli方程。

（a）Alexis Clairaut（b）Jacob Bernoulli（c）Leonhard Euler
图 1-4 Clairaut、Bernoulli、Euler的画像注：图像均来源于维基百科 
1768年，瑞士数学家 Leonhard Euler（1707.1783，图 1-4（c）
出版了其著名的
三卷著作“Institutionum Calculi Integralis”（积分学基础），其中有很多关于微分方程的原创性成果，包括著名的积分因子法与全微分方程的求解，也包括后来称为 Euler算法的最早的微分方程数值解法。
早期微分方程的研究侧重于解析解与性质的研究，因为没有计算机这类研究工具的支持，所以微分方程的数值求解并没有得到足够的重视。
德国数学家 Carl David Tolmé Runge（1856.1927，图 1-5（a）
在 Euler算法的基础上增加了一些中间插值点，旨在提高 Euler数值算法的精度，于 1895年提出了微分方程的求解算法，并将算法拓展成微分方程组的求解。1900年德国数学家 Karl Heun（1859.1929）在 Runge算法的基础上进一步整理了原始烦琐的算法，给出了基于 Gauss求积公式的更一般性的方法。1901年，德国数学家 Martin Wilhelm Kutta（1867.1944，图 1-5（b）
将 Runge与 Heun所提出的算法拓展到了更一般的高阶规范格式，并推导出 5阶算法。法国学者 A.Huˇta在 1956年发表了 6阶算法 [6]。 1963年，新西兰数学家 John Charles Butcher（1933.）在前人工作的基础上给出了更一般的 Runge–Kutta表格，又称 Butcher表格 [7]。1968年的一篇 NASA研究报告中，Erwin Felhberg提出了 7阶和 8阶的 Runge–Kutta算法计算公式 [8]。奥地利数学家 Ernst Hairer（1949.）在 1978年提出了 10阶 Runge–Kutta算法的计算公式 [9]。这些算法统称为 Runge–Kutta算法。

（a）Carl Runge（b）Martin Kutta（c）John Adams
图 1-5 Runge、Kutta和 Adams的照片注：图像均来源于维基百科
从算法的分类看，Runge–Kutta算法属于单步数值方法。除了单步法之外，还有例如 Adams算法的一大类线性多步方法。Adams算法是以英国数学家、天文学家 John Couch Adams（1819.1892，图 1-5（c）
命名的，与 Runge–Kutta算法相比有着更悠久的历史 [10]。Adams类算法是利用已知的几个点，逐步递推得出微分方
·8·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
程的数值解。但是，多步方法起始阶段的几个点处的解如何获得？如果这几个起始点问题不能被很好地解决，将影响整个求解过程的精度与效率。
陆续出现的计算机数学语言与工具都配备了强大的微分方程求解工具，可以很容易地求出微分方程的数值解甚至解析解。本书的目标是指导读者如何使用 MATLAB求解各种各样的微分方程。 

1.3本书主要内容
本书将侧重介绍如何利用 MATLAB直接求解微分方程问题。
第 1章给出了一些实际问题的微分方程建模的简单例子。
第 2章主要介绍微分方程的解析解方法。首先介绍微分方程求解的一些基础的理论知识，然后侧重介绍如何利用 MATLAB求出常系数线性微分方程、一般微分方程以及矩阵微分方程的解析解。此外，还探讨了某些极特别的非线性微分方程的解析解问题。
本书从第 3章开始介绍微分方程的数值解方法。第 3章及后面的章节介绍微分方程初值问题的数值求解方法。首先简要介绍一些常规的数值算法及 MATLAB实现，然后介绍 MATLAB中一阶显式微分方程组的数值求解方法，并介绍微分方程数值解的验证方法。
第 4章主要介绍将各类一般微分方程变换成一阶显式微分方程组的方法。介绍单个微分方程、复杂高阶微分方程的变换方法，还介绍微分方程组及矩阵型微分方程的变换方法。如果能够成功变换，则可以调用第 3章介绍的求解程序得出微分方程的数值解。
第 5章介绍几种不适合用普通微分方程求解函数进行求解的问题及其求解方法，包括刚性微分方程的求解方法、隐式微分方程的求解方法及微分代数方程的求解方法，还介绍开关微分方程与随机线性微分方程的数值求解方法。
第 6章介绍延迟微分方程的一般形式，并介绍一般延迟微分方程、变延迟微分方程及中立型延迟微分方程的数值解法。
第 7章介绍微分方程的性质与行为分析，系统地分析稳定性、周期性等微分方程的性质，探讨极限环、混沌、分岔等特殊的微分方程行为，并介绍非线性微分方程的平衡点与线性化方法。
第 8章介绍分数阶微分方程的求解方法。首先介绍分数阶微积分的基本定义与性质，然后分别介绍分数阶线性微分方程、一般分数阶非线性微分方程的数值解方法。特别地，还介绍这类微分方程的高精度数值解方法。
第 9章介绍基于 Simulink的微分方程框图求解方法。首先介绍 Simulink环境的入门知识，通过例子演示一般微分方程、微分代数方程、切换微分方程以及随机微分方程和延迟微分方程的建模与仿真方法，同时介绍各类分数阶微分方程的通用建模与仿真方法。
第 10章介绍微分方程的边值问题求解。首先介绍打靶方法的基本思路与求解算法，然后介绍利用 MATLAB现有函数的边值问题求解方法，最后介绍基于最优化技术的边值问题求解方法，得出用传统方法难以求解的边值问题的解。
第 11章简单介绍偏微分方程的数值解方法。首先介绍扩散偏微分方程的数值求解方法，然后介绍利用 MATLAB偏微分方程工具箱对几类特殊形式的偏微分方程进行求解的方法。
本书侧重介绍如何求解微分方程问题，如果能够得出微分方程的解析解，则尽量得出问题的解析解；如果不能得出解析解，则得出微分方程的数值解。总之，本书的宗旨是指导读者如何充分利用 MATLAB的强大功能，以最可靠的方法研究微分方程。
本章习题 
1.1如图 
1-1所示的电路中，如果原电阻元件被替换为电阻 R与另一个电容 C1的串联形式，试重新列写微分方程。如果 R与 C1并联，微分方程又该如何列写？ 

1.2在 
Malthus人口模型中，如果某年 t0处的初始人口为 60亿，且已知自然增长率为 2％，试求出 10年后的人口数量。100年之后呢？ 

1.
3试用幂级数法求解 Newton微分方程式（1-2-1）的前 8项。 

1.
4如果你有微分方程的基础，已知函数 f(c) = 5(c3 . c)/2，试求解 Clairaut的隐式微分方程式（1-2-3）。


常微分方程的解析解  
第 2章 
MATLAB微分方程求解 


第 1章介绍了几类问题的常微分方程建模。如果已经为某一个动态系统或某种现象建立了微分方程模型，则下一步需要求解该微分方程，得到系统的响应。求解微分方程比求解代数方程要复杂得多，通常需要更多的技巧和专门的求解方法。 
2.1节主要探讨一阶微分方程的解析解方法，介绍利用底层方法求解不同类型的一阶微分方程的解析解的思路，得出可解方程的解析解。 2.2节探讨二阶线性微分方程的解析解求解方法，底层方法同样需要各种各样的技巧，并介绍几种常用的特殊函数的定义与求解方法。 2.3节介绍基于 Laplace变换的微分方程解析解方法。利用 Laplace变换工具，将常系数线性微分方程映射成代数方程并直接求解，同时探讨非零初值线性微分方程的直接求解方法。 2.4节介绍符号运算工具箱的常微分方程直接求解方法，利用这类方法能够更容易、更方便地求解复杂的常微分方程、时变微分方程等一般线性微分方程。 2.5节介绍线性微分方程、线性状态方程与 Sylvester微分方程的求解方法。 2.6节尝试通过 MATLAB现成的求解函数直接求解非线性微分方程。
值得注意的是，只有极少数非线性微分方程是存在解析解的，绝大多数微分方程不存在解析解，应该考虑后面章节将介绍的数值求解方法。
本章前几节主要介绍微分方程求解过程的手工推导方法，为读者初步了解微分方程的底层求解方法作必要的铺垫。如果读者只关心如何利用计算机求解微分方程，则可以略过这些内容，从 2.4节开始阅读。 

2.1一阶微分方程的解析解
本节首先探讨一阶线性微分方程的解析解求解方法，先讨论几种特殊的形式，例如简单方程、齐次与非齐次线性微分方程，再给出可分离变量的非线性微分方
程的解法。
2.1.1可由简单积分求解的微分方程
一阶微分方程的种类千差万别，给方程的求解带来了巨大的挑战。某些微分方程的求解需要掌握并尝试各种各样的求解技巧，而很多方程即使掌握了某些技巧也不能一定真正找出其解析解。这里介绍一些简单的一阶微分方程，并给出可以尝试的求解方法。
　　定义 2-1最简单的一类一阶微分方程的数学形式为 
dyd(xx)= f(x)	（2-1-1）
也可以将其简写成 y ′ (x)= f(x)的形式。这时，微分方程的解析解实际上就是 
f(x)的不定积分，即 y(x)= ∫ f(x)dx + C	（2-1-2）
这类微分方程的解就是一元函数的不定积分。如果不定积分有解析解，则微分方程有解析解，否则微分方程是没有解析解的。
　　例 2-1试求解下面的一阶微分方程。 
d cos x sin x(2x + 4) 
y(x)= .
dxx2 +4x +3 (x2 +4x + 3)2　　解对照式（2-1-1）给出的数学形式可以看出，题目中方程等号右侧就是 f(x)。所以，可以输入原函数，并直接进行积分，由下面的语句求取微分方程的解析解。 
>> 	syms x 
f(x)=cos(x)/(x^2+4*x+3)-sin(x)*(2*x+4)/(x^2+4*x+3)^2; 
Y=simplify(int(f)) %求积分并化简

得出的解析解为 sin x/(x2 +4x + 3)。其实，该解只是不定积分的原函数（primitive function），微分方程的解应该是该解加上任意常数 C。对 C取几个不同的值，可以绘制出微分方程的解曲线，如图 2-1所示。可见，这样得出的微分方程的解是某一个解通过上下平移得到的。 
>> 	fplot([Y,Y+1,Y+2,Y+3],[0,10]) 


2.1.2齐次线性方程
前面提到的可以直接通过积分求解的微分方程不是真正意义上的微分方程。一般微分方程中还应该含有 y(x)函数本身，这样的微分方程就不能通过不定积分直接求解了。这里先讨论一种简单的一阶齐次线性微分方程，并给出其一般的求解方法。
　　定义 2-2一阶齐次（homogeneous）线性微分方程的一般形式为 
dyd(xx)+ f(x)y(x)=0	（2-1-3）
·12·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 

图 2-1微分方程的解曲线
经过简单变换，可以由下面的积分得出该微分方程的解析解。 
dyy((xx)) = .f(x)dx（2-1-4）
由此可以推导出 ln y(x)= . ∫ f(x)dx（2-1-5）
最终得出微分方程的解析解为 
∫ 
y(x)= Ce. f(x)dx（2-1-6）
对一阶微分方程而言，方程的解中存在一个待定系数 C，如果知道函数 y(x)在某一个时刻的确切值，则可以通过求解代数方程的方式唯一地确定 C的值。高阶微分方程可能含有多个待定系数。如果已知解函数 y(x)或其导数的若干值，则可以唯一地确定方程的解析解。
　　例 2-2试求解微分方程 y ′ (x)= .xy(x)。　　解对于这种可以分离的微分方程，可以将其变换为 y ′ (x)/y = .x。对其左侧与右侧分别积分计算 
>> syms x y(x); L=int(diff(y)/y), R=int(-x)

则可以得出方程的解为 
2 
2
ln y(x)= .x+ C1 ≈ y(x)= Ce.x /2 
2 

2.1.3非齐次线性方程
如果微分方程右边增加一项 x的函数 g(x)，则原始的齐次微分方程就变成一般的非齐次微分方程。本节将探讨一阶非齐次线性微分方程的解析解的求解思路。
　　定义 2-3一阶非齐次（inhomogeneous）线性微分方程的一般形式为 
dyd(xx)+ f(x)y(x)= g(x)（2-1-7）
手工求解这类方程比较烦琐，需要使用特殊的技巧。例如，找到一个辅助函数 
.(x)，使得 [dy(x) ] d ().(x) dx + f(x)y(x)= dx.(x)y(x)（2-1-8）
这样，原方程可以变换为 
[]
dd x	.(x)y(x)= .(x)g(x)（2-1-9）
对得到的方程两边作区间 [x0,x]的定积分运算，则 
∫
x 
.(x)y(x) . .(x0)y(x0)= .(t)g(t)dt（2-1-10）
x0 
最终，可推导出方程的解析解为 ∫
x 
y(x)= ..((xx0)) y(x0)+ .(1 x) .(t)g(t)dt（2-1-11）
x0 
可见，寻找辅助函数 .(x)是求解微分方程的关键一步。寻找辅助函数的方法是
与 f(x)表达式密切相关的，需要求解下面的积分问题。 
(∫ )
.(x)= exp f(x)dx	（2-1-12）
　　例 2-3试求解下面一阶微分方程非齐次方程的解析解 [11]。 
dy(t)+ eγt 
y(t)= keγt ,y(0) = y0
dt 　　解对比定义 2-3中的标准型可见， f(t)= eγt，g(t)= keγt，x0 =0。这样，就可以利
用下面的 MATLAB语句直接计算出微分方程的解析解。 
>> 	syms t x k lam y0 
f(t)=exp(lam*t); g(t)=k*exp(lam*t); phi=exp(int(f)) 
y=phi(0)*y0/phi+1/phi*int(phi(x)*g(x),x,0,t)

得出方程的解析解与辅助函数分别为 
/γ
y(t)= k +(y0 . k) e(.e>t.1)/γ ,.(t)= ee>t
当然，这里给出的解析解求解方法比较烦琐，所以不对该方法作过多的解释，
后面将介绍利用 MATLAB的更通用、直接的求解方法。 

2.1.4可分离变量的非线性微分方程
　　定义 2-4一类可以分离变量的一阶非线性微分方程可以写成 
dyd(xx)= g(y(x))f(x),y(x0)= y0（2-1-13）
　　定理 2-1定义 2-4中可分离变量的非线性微分方程的解析解为 
∫∫
y(x) x 
gd(ξξ )= f(t)dt	（2-1-14）
y(x0) x0 
·14·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
如果式（ 2-1-14）中的两个积分表达式可以求出，则可以得出原始方程的解析解。然而，绝大部分非线性微分方程是不能得出解析解的，一般不可分离变量的非线性微分方程也是没有解析解的。
　　例 2-4试求解下面的一阶非线性微分方程。 
dy(t) 
+8y(t)+ y 2(t)= .15,y(0) = 0
dt 　　解由定义 2-4中的非线性方程可见， g(y) = 15+8y + y2，f(x)=1。这样，由下面的语句可以直接求解微分方程。 
>> 	syms x t u y 
g(y)=15+8*y+y^2; f(t)=-sym(1); 
G=simplify(int(1/g(u),u,0,y)), F=int(f,0,t) 
eq=G==F, y=solve(eq,y) %构造并求解隐式方程

其中，式（2-1-14）左侧 G为 arctanh(4) . arctanh(y + 4)，右侧为 .t。由此，可以构造出隐式方程 arctanh(4) . arctanh(y +4) = .t，求解该方程则可以得出微分方程的解析解为 y(t)= tanh(t + arctanh(4)) . 4。
其实，从求取积分 G的结果可见，若手工求解，则可能对 1/g(y)采用部分分式展开的形式进行求解，与这里得到的 arctanh(·)不同。如果想得到这类结果，可以改写 G的形式。例如，可以使用函数 rewrite()改变其默认的显示形式。 
>> 	G=rewrite(G,'exp'), F=int(f,0,t) 

eq=G==F, y=simplify(solve(eq)) %构造并求解隐式方程这时得到的隐函数方程为 ln(y + 3)/2 . ln(y + 5)/2 . ln 3/2 + ln 5/2 = .t，则由该方程可以得出微分方程的解析解为 y(t) = 6/(5e2t . 3) . 3，这个结果与前面得出的结果是完全等效的。
可以看出，对一般一阶线性微分方程而言，单纯依赖手工求解是极其麻烦的。求解过程涉及复杂函数的积分，在很多情况下是不可求解的，需要引入特殊函数来描述方程的解析解。非线性方程的求解可能更加麻烦，且绝大多数非线性微分方程不存在解析解，在尝试获得某些微分方程的解析解时，数学家们发明了各种各样的特殊函数。尽管如此，仍有大部分微分方程不存在解析解，所以从第 3章开始将系统地探讨微分方程的数值解算法与 MATLAB应用。 


2.2特殊函数与二阶线性微分方程
与一阶微分方程相比，二阶微分方程的求解难度更大。本节将探讨几类常用二阶线性微分方程的数学模型，给出这些微分方程的“解析解”，并介绍几种常用的特殊函数。　　定义 2-5二阶线性微分方程的一般形式为 
d2a(x) dyx(2 x)+ p(x) dyd(xx)+ q(x)y(x)= f(x)（2-2-1）
其中，系数函数 a(x)、p(x)、q(x)与 f(x)均为已知函数。
随着各个系数函数的不同变化，二阶线性微分方程的求解方法也各不相同。文献 [12]列出了 1000多种不同的系数选择与方程的求解方法，绝大多数是假设 f(x)=0的齐次微分方程。对很多无解析解的线性微分方程，数学家们发明了许多特殊函数，用来描述某类微分方程的解析解。
为更好地理解某些微分方程的解析解，有必要介绍一下常用的特殊函数。比较常用且与微分方程密切相关的特殊函数包括 Gamma函数、超几何函数、各类 Bessel函数、 Legendre函数与 Airy函数等。有些特殊函数对应于特定的二阶线性微分方程，本节对这些方程与函数作一个简单的介绍。 
2.2.1 Gamma函数
　　定义 2-6 Gamma函数是由下面的无穷积分定义的： 
∫← 
f(z)= e.ttz.1dt	（2-2-2）
0 
　　定理 2-2 Gamma函数的一条重要的性质是 
f(z +1) = zf(z)	（2-2-3）
　　定理 2-3作为一个特例，对非负整数 z，可以由式（ 2-2-3）直接得出阶乘公式。 
f(z +1) = zf(z)= z(z . 1)f(z . 1) ··· = z!（2-2-4）
可以认为 Gamma函数是阶乘的插值，或可以认为 Gamma函数是阶乘在非整数 z上的推广。如果 z为负整数， f(z + 1)趋近于 ±∞。在 MATLAB下，可以由函数 y=gamma(x)直接计算 Gamma函数。如果 x为向量，得出的 y也是向量。此外， x也可以是矩阵或其他数据结构，得到的 y与 x结构相同。注意，这里的函数 gamma()只能处理实数输入变量。
　　例 2-5试绘制 x在区间 (.5, 5)的 Gamma函数曲线。　　解可以由下面的语句直接绘制 Gamma函数曲线，如图 2-2所示。由于 Gamma函数在 z =0，.1，.2，···各点上趋于无穷大，调用函数 ylim()限定 y轴范围，使得图形更有意义。 
>> 	a=-5:0.002:5; plot(a,gamma(a)), ylim([-15,15]) 
hold on, v=[1:4]; plot(v,gamma(v),'o'), hold off 

·16·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 

图 2-2 Gamma函数曲线
对一些特定的 z值，可以计算出 Gamma函数的值。 
(	) () ∨() ∨() ∨
3 	丘 
f	1 = ∨丘, f3= 丘, f5= 丘, f 7 =15 , ···（2-2-5）
2 2224 2 8 
　　例 2-6试用 MATLAB证明式（2-2-5）。　　解式（2-2-5）中的 Gamma函数值可以由下面的语句直接求出。 
>> 	syms t z; Gam=int(exp(-t)*t^(z-1),t,0,inf); 
I1=subs(Gam,z,sym(1/2)), I2=subs(Gam,z,sym(3/2)) 
I3=subs(Gam,z,sym(5/2)), I4=subs(Gam,z,sym(7/2))

　　定理 2-4更一般地，如果 z不是整数，则该通项可以拓展成 
()
f 	z + m1 =(2丘)(m.1)/2m1/2.mz f(mz)（2-2-6） 
2.2.2超几何函数
超几何函数是一种常用的特殊函数，也是很多特殊函数的基础。这里首先给出超几何函数的定义，然后介绍超几何函数的计算方法。
　　定义 2-7超几何函数（hypergeometric function）的一般形式 [13]为 

(	)
pFqa1,a2, ··· ,ap; b1,b2, ··· ,bq; z
←
f(b1)f(b2) ··· f(bq) ∑ f(a1 + k)f(a2 + k) ··· f(ap + k) zk （2-2-7）= 
f(a1)f(a2) ··· f(ap) f(b1 + k)f(b2 + k) ··· f(bq + k) k!
k=0 式中， bi不能是非正整数。如果 p . q，则该函数对所有的 z都是收敛的；如果 p = q +1，则函数在 |z| < 1时收敛；若 p>q +1，则对所有的 z函数都将发散。 
MATLAB符号运算工具箱提供的函数 hypergeom()可以用于计算超几何函数 pFq(a1,a2, ··· ,ap; b1,b2, ··· ,bq; z)，其调用格式为 y=hypergeom([a1,a2,··· ,ap],[b1,b2,··· ,bq],z) 
例如，超几何函数 1F1(a; b; z)可以用函数 y=hypergeom(a,b,z)直接计算，而函数 2F1(a, b; c; z)可以由 y=hypergeom([a,b],c,z)直接计算。
　　例 2-7试绘制 Gauss超几何函数 2F1 (1.5, .1.5; 1/2; (1 . cos x)/2)的曲线。　　解超几何函数 2F1 (1.5, .1.5; 1/2; (1 . cos x)/2)可以简化成 cos 1.5x。如果使用函数 hypergeom()，则可以绘制出该超几何函数的曲线如图 2-3所示。可以看出，该曲线与 cos 1.5x完全重合。 
>> 	syms x, y=hypergeom([1.5,-1.5],0.5,0.5*(1-cos(x))); 
fplot([cos(1.5*x),y],[-pi,pi]) 

f(x)

图 2-3例 2-7的超几何函数曲线
　　定理 2-5作为一个特例，超几何函数 2F1(a, b; c; z)为下面微分方程的解。 z(1 . z) d2dyz(2 z)+ [c . (a + b + 1)]dyd(zz) . aby(z)=0（2-2-8）这里特意将自变量写成 z，表示自变量可以为复数变量。 
2.2.3 Bessel微分方程
如果式（ 2-2-1）中的二阶线性微分方程各个系数选择不同的形式，则可以得到不同形式的特殊微分方程，例如常用的 Bessel微分方程。本节首先介绍 Bessel微分方程的数学形式，然后介绍各类 Bessel函数及其计算方法与曲线绘制。
　　定义 2-8如果二阶齐次线性微分方程的一般形式可以写成 
x 2 d2dyx(2 x)+ x dyd(xx) +(x 2 . η2)y(x)=0（2-2-9）则该方程称为 Bessel微分方程。　　定理 2-6 η阶 Bessel微分方程的“解析解”的通解为 
y(x)= C1Jδ(x)+ C2J.δ(x)	（2-2-10）
·18·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
其中，C1和 C2为任意常数，Jδ (x)为第一类 η阶 Bessel函数，定义为 
←
Jδ(t)= ∑ (.1)m 2δ+2mtδ+2m （2-2-11）
m! f(η + m + 1)
m=0 
式中，f(·)为 Gamma函数。由于原方程真实意义下的解析解不存在，所以数学家发明了相应的特殊函数
定义微分方程的解析解。这就像 e.t2这类函数的不定积分不存在，所以数学家发明了特殊函数 erf(·)研究其解析解。　　定理 2-7若 η为正整数，第一类 Bessel函数有如下性质： 
Jδ(x)=(.1)δJ.δ (x)（2-2-12） 
Jδd(xx)= ηx Jδ (x) . Jδ+1(x)（2-2-13） 
∫
x 
tδJδ.1(t)dt = x δJδ(x)（2-2-14）
0 
　　定理 2-8 η阶第一类 Bessel函数是一种特殊的超几何函数。 
()
2 
Jδ(x)= (x/2)δ 0F1 η + 1; . x（2-2-15）
f(η +1) 4 
Bessel函数是以德国天文学家、数学家 Friedrich Wilhelm Bessel（1784.1846）命名的。 Bessel函数是由瑞士数学家 Daniel Bernoulli（1700.1782）首先提出并由 Bessel推广的特殊函数。除了前面介绍的第一类 Bessel函数之外，还可以使用第二类与第三类 Bessel函数。
　　定义 2-9 η阶第二类 Bessel函数 Yδ(t)的一般形式为 
Yδ (t)= Jδ(t) cossin ηt ηt. J.δ(t) （2-2-16）
第二类 Bessel函数又称为 Neumann函数，是以德国数学家 Carl Gottfried Neu-mann（1832.1925）命名的。
　　定义 2-10第三类 Bessel函数 Hδ(x)的数学形式为 
H(1) δ (x)= Jδ(x)+ jYδ(x), H(2) δ (x)= Jδ(x) . jYδ(x)（2-2-17） 
MATLAB提供了函数 besselj()、bessely()与 besselh()，分别用于求解第一类、第二类与第三类 Bessel函数。这些函数的调用格式为 
y=besselj(η,x); y=bessely(η,x); y=besselh(η,k,x)
其中，k =1，k =2对应于第三类 Bessel函数的两种情形。
　　例 2-8试绘制 η =0, 1, .1, 2时的第一类 Bessel函数曲线。　　解选择几个不同的阶次，则可以由下面的语句直接绘制 Bessel函数曲线，如图 2-4所示。可见，J1(x)= .J.1(x)。 
>> 	x=-5:0.1:5; 
y1=besselj(0,x); y2=besselj(1,x); y3=besselj(-1,x); 
y4=besselj(2,x); plot(x,y1,x,y2,x,y3,x,y4) 


图 2-4不同阶次的 Bessel函数曲线 



2.2.4 Legendre微分方程与 Legendre函数 
Legendre函数满足 Legendre微分方程。下面将给出 Legendre微分方程与 Leg-endre函数的定义，然后给出相应函数的计算方法。 Legendre方程与 Legendre函数是以法国数学家 Adrien-Marie Legendre（1752.1833）命名的。
　　定义 2-11 Legendre微分方程的数学形式为 
[	]
2 
(1 . z 2) d2dyz(2 z) . 2z dyd(zz)+ λ(1 + λ) . 1 . μz2 y(z)=0（2-2-18）
其中，复数变量 λ与 μ分别称为 Legendre函数的度数与阶次。
　　定义 2-12 Legendre函数的数学形式可以由超几何函数定义。 
μ/2
1 1+ z 	1 . z
Pμγ(z) = f(1 . λ) (1 . z)2F1 (.λ, λ + 1;1 . μ;2 )（2-2-19）
收敛域为 |1 . z| < 2，μ =0, 1, 2, ··· ,λ。
　　例 2-9试绘制 μ =3时 Legendre函数的曲线。　　解取 μ =3，则可以用下面的命令绘制出 4条不同的 Legendre函数的曲线，如图 2-5所示。当然，可以考虑使用下面的命令绘制不同的 Legendre函数曲线。 
>> 	x=-1:0.04:1; Y=legendre(3,x); plot(x,Y) 

·20·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 

P(x) 
P(x)
P(x) 
2333
P(x) 
图 2-5不同阶次的 Legendre函数曲线 
13

2.2.5 Airy函数 
Airy函数是以英国数学家、天文学家 George Biddell Airy（1801.1892）命名的一类特殊函数，该函数定义为 Airy微分方程的解析解。本节首先给出 Airy微分
03
方程的一般形式，再介绍 Airy函数的计算与曲线绘制。
　　定义 2-13 Airy微分方程是形式简单的方程，其数学形式为 
d2dyz(2 z) . zy(z)=0	（2-2-20）
其中，自变量 z可以是复数变量。
　　定义 2-14 Airy函数与第二类 Airy函数的定义如下： 
Ai(z)= 丘 1∫0 ← cos (t3 3 + zt)dt	（2-2-21） 
∫← 
Bi(z)= 丘 10 [exp (. t3 3 + zt) + sin (t3 3 + zt)]dt（2-2-22） 
MATLAB提供了函数 z=airy(k,z)，可以直接求取复数矩阵的 Airy函数值。 k =0，1时分别计算 Ai(z)及其一阶导数；k =2，3时分别计算 Bi(z)及其一阶导数。如果只想计算函数 Ai(z)，则可以省略变元 k。当然，若输入变元为实数矩阵，则得出的结果也为同维度的实数矩阵。
　　例 2-10 Airy函数是复数函数，试绘制 Airy函数实部在复平面的等高线图。　　解可以生成复数网格矩阵 z = x + jy，则可以得出 Airy函数的矩阵。为增加分辨力，对得出的复数矩阵作 ±4处理，最终绘制出如图 2-6所示的等高线图。 
>> 	[x,y]=meshgrid(-4:0.1:4); z=x+y*1i; y1=airy(z); 
y1=real(y1); ii=find(y1>=4); y1(ii)=4; 
ii=find(y1<=-4); y1(ii)=-4; contourf(x,y,y1,80) 


图 2-6 Airy函数实部的等高线图正如前面指出，对二阶线性微分方程的各个系数函数的不同取值，可能构造出不同的微分方程，致使二阶微分方程的手工求解变得极其困难，而二阶非线性微分方程或三阶以上的线性、非线性微分方程的求解更加复杂，只能求解极简单形式的微分方程。后续内容将利用成熟的微分方程解析求解工具，直接求解二阶或高阶的微分方程，可以利用统一的函数直接求解，尝试得出微分方程的解析解。如果微分方程的解析解不存在，则下一步应考虑如何得出所研究微分方程的数值解。 


2.3常系数线性微分方程的求解
高阶变系数非线性微分方程是不容易求得解析解的，绝大多数非线性微分方程不存在解析解，所以这里主要讨论常系数线性微分方程的求解方法。首先，给出常系数线性微分方程的一般数学形式，然后简要介绍 Laplace变换以及基于 Laplace变换的常系数线性微分方程的求解方法，最后介绍非零初值微分方程组的一般解析求解方法。 
2.3.1线性常系数微分方程解析解的数学描述
这里给出常系数线性微分方程的数学形式，在后续内容中，采用不同的方法探讨这类微分方程的求解方法。
　　定义 2-15假设已知常系数线性微分方程的一般描述为 
dny(t) dn.1y(t) dn.2y(t) dy(t)
+ a1 + a2 + ··· + an.1 + any(t)
dtn dtn.1 dtn.2 dt 
dmu(t) dm.1u(t) du(t) （2-3-1）= b1 + b2 + ··· + bm + bm+1u(t)
dtm dtm.1 dt 
其中，ai,bi均为常数，n称为微分方程的阶次。
·22·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
若把微分方程看成一个动态系统模型，则可以将信号 y(t)看成系统的输出信号，将信号 u(t)看成系统的输入信号。这样，微分方程描述的就是系统输入信号与输出信号之间的动态关系。 
2.3.2基于 Laplace变换的求解方法
在一般微分方程的著作中，提到线性微分方程的解析解法都建议将其变换为代数方程，求出特征根，再重构出微分方程的通解。然而，很少有著作明确说明为什么要先求解代数方程。这里主要介绍基于 Laplace变换的求解方法。
　　定义 2-16时域函数 f(t)的 Laplace变换的定义为 
∫← 
L [f(t)] = f(t)e.stdt = F (s)（2-3-2）
0 
其中，L [f(t)]为 Laplace变换的简单记号。
　　定义 2-17如果已知函数的 Laplace变换表达式 F (s)，则可以通过下面的反变换公式反演求出其 Laplace反变换。 
1 ∫η+j← 
f(t)= L .1[F (s)]= 2丘j η.j← F (s)estds（2-3-3）
其中，ν大于 F (s)的奇点实部。
　　定理 2-9 Laplace变换的微分性质是微分方程应用中的重要性质。 
[] L [dfd(tt) ] = sF (s) . f(0+)（2-3-4）
一般地，n阶微分可以由下式求出。 
L dndftn (t)= s nF (s) （2-3-5）n.2f
.s n.1f(0+).s ′ (0+) .· ··. f(n.1)(0+)
若假设函数 f(t)及其各阶导数的初值均为 0，则式（ 2-3-5）可以简化成 
L [dndftn (t) ] = s nF (s)（2-3-6）
利用式（ 2-3-6）介绍的性质，对零初值问题有 L [dmy(t)/dtm]= smL [y(t)]，可以对应得出下面的多项式代数方程。 
s n + a1s n.1 + a2s n.2 + ··· + an.1s + an =0（2-3-7）
　　定理 2-10假设代数方程的特征根 si均可以求出，且假设它们均相异，则可以得出原微分方程解析解的一般形式为 
y(t)= C1er1t + C2er2t + ··· + Cnernt + Ψ(t)（2-3-8）
其中， Ci为待定系数， Ψ(t)为满足 u(t)输入的特解。 si有重根的情况也有相应的解析解形式。
　　定理 2-11如果特征方程含有 m阶重根 ri，则可以按照下面的形式构造出通解对应的项 ()
C1 + C2t + ··· + Cmtm.1e.rit（2-3-9）
从以上介绍的两个定理可见，如果能得出特征方程的根，则可以手工构造出微分方程解析解的一般形式。
根据代数方程式（ 2-3-7），由著名的 Abel–Ruffini定理可知，四次及以下阶次的多项式方程是能求出根的解析解的。故可以得出结论，低阶常系数线性微分方程有一般意义下的解析解。结合多项式代数方程的准解析解法，可以得出一般高次多项式代数方程的高精度数值解，构造出高阶常系数线性微分方程的准解析解方法。
　　例 2-11试求解下面的齐次常系数线性微分方程。 
y(4)(t) + 10y ′′′ (t) + 35y ′′ (t) + 50y ′ (t) + 24y(t)=0
　　解如果微分方程输出信号 y(t)及其各阶导数的初值均为零，则依据 Laplace变换的性质，可以推导出如下方程。 
s4Y (s) + 10s3Y (s) + 35s2Y (s) +50sY (s) + 24Y (s) =(s4 + 10s3 + 35s2 + 50s + 24)Y (s)=0
求解特征方程 s4 + 10s3 + 35s2 + 50s +24 = 0，可以得出其解为 s1 = .1，s2 = .2， s3 = .3，s4 = .4。因此，可以写出微分方程的通解为 
y(t)= C1e.t + C2e.2t + C3e.3t + C4e.4t
　　例 2-12试求出下面齐次微分方程的解析解通解。 
y(5)(t) + 12y(4)(t) + 57y ′′′ (t) + 134y ′′ (t) + 156y ′ (t) + 72y(t)=0
已知，y(t)及各阶导数的初值为零。　　解可以手工写出特征方程为 s5 + 12s4 + 57s3 + 134s2 + 156s +72 = 0，运用下面的语句可以直接求解对应的代数方程。 
>>symss; %声明符号变量 
F=s^5+12*s^4+57*s^3+134*s^2+156*s+72; %描述多项式 
r=solve(F) %求解多项式方程

可以得出特征根为 s1 = .2，s2 = .2，s3 = .2，s4 = .3，s5 = .3，可以手工写出微分方程的通解为 y(t)=(C1 + C2 + C3t2)e.2t +(C4 + C5t)e.3t 
·24·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 


2.3.3非齐次微分方程的求解
前面给出的例子都假设微分方程为齐次方程，所以只需求解 y(t)各阶导数的线性组合，即可直接写出微分方程的通解。如果方程为非齐次的，则可以考虑对方程右侧直接求取 Laplace变换，就可以写出输出信号的 Laplace变换表达式。得出该表达式后，对其结果求取 Laplace反变换，可以直接得出微分方程的解析解。
如果不使用计算机工具，需要使用部分分式展开的方法，将得出的输出信号 Laplace变换表达式写成简单的分式形式，再写出对应的微分方程解析解。采用 MATLAB求解工具，就可以直接跳过部分分式展开过程，使用函数 ilaplace()，得到微分方程的解析解。下面通过例子演示直接求解的方法。
　　例 2-13重新考虑例 2-12中的微分方程，如果该方程为下面给出的非齐次方程，试得出该微分方程的解析解。 
y(5)(t) + 12y(4)(t) + 57y ′′′ (t) + 134y ′′ (t) + 156y ′ (t) + 72y(t)= e.t sin t
　　解考虑利用下面的语句直接写出输出信号的 Laplace变换表达式，直接使用求解语句，可以得出微分方程的解析解。根据解析解可以绘制出输入信号与输出信号的曲线，如图 2-7所示。由于输出信号的幅值较小，所以在图中绘制的是 20倍的 y(t)信号。 
>> 	syms s t; u(t)=exp(-t)*sin(t); 
F=s^5+12*s^4+57*s^3+134*s^2+156*s+72; 
U=laplace(u); Y=U/F; y=ilaplace(Y), y=simplify(y) 
fplot([u,20*y],[0,10]) %绘制输入与输出信号曲线 


图 2-7微分方程解的输入、输出曲线
由上面语句可以得出微分方程的解析解为 
31911 1 	1 
2e.2t
y(t)= e.2t . e.3t . te.2t . te.3t + e.t (cos t . 7 sin t)+ t
42525100 	4 
手工合并同类项并化简，可以得出方程的解析解为 
1 	11 
y(t)= (3 . 2t + t2)e.2t . (18 + 5t)e.3t + e.t (cos t . 7 sin t)
4 	25 100 

2.3.4非零初值的微分方程求解
如果微分方程的初值不是零，则式（ 2-3-6）中的表达式不成立。利用式（ 2-3-5）将微分方程变换成代数方程，求解所得的代数方程，对得到的结果进行 Laplace反变换，则可以得到原始方程的解析解。下面举例演示这样的求解方法。
　　例 2-14假设输入信号为 u(t)= e.5t cos(2t + 1) + 5，已知初值 y(0) = 3，y ′ (0) = 2， y ′′ (0) = y ′′′ (0) = y(4)(0) = 0，试求出下面微分方程的解析解。 
y(4)(t) + 10y ′′′ (t) + 35y ′′ (t) + 50y ′ (t) + 24y(t)= u(t)
　　解为简单起见，考虑用下面的语句直接写出方程左侧表达式的 Laplace变换形式。其中， eqn是 Y (s)的多项式， eqn1是与初始值有关的表达式，这样就可以得到 Y (s)的表达式。调用函数 ilaplace()，可得出微分方程的解析解。绘制出输出信号的曲线，如图 2-8所示。 
>> 	syms t s; y0=3; y1=2; u=exp(-5*t)*cos(2*t+1)+5; 
eqn=s^5+12*s^4+57*s^3+134*s^2+156*s+72; 
eqn1=(s^4*y0+s^3*y1)+12*(s^3*y0+s^2*y1)+... 

57*(s^2*y0+s*y1)+134*(s*y0+y1)+156*y0; 
Y(s)=(laplace(u)+eqn1)/eqn; y=ilaplace(Y); 
y=simplify(y), fplot(y,[0,10]) %化简并绘制解的曲线


图 2-8微分方程的输出曲线
求解微分方程，则可以得出方程的解析解为 
(	)()
1642 cos 1 1372 sin 1 4203 3 cos 15 sin 1 4702 
y(t)= e.2t . + . e.3t . + 
2197 21978 4 89 
()(	)
cos 1 sin 1 475 83 cos 1 64 sin 1 1425 
.te.3t . + . te.2t . + 
4 3 1691694 
(	)()
sin24t (9 cos 1 + 46 sin 1) 23 cos 19 sin 1 
+e.5t 
cos 2t .	. 
46 cos 1 . 9 sin 1 8788 17576 
·26·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 
t2e.2t 
3 cos 12 sin 1 451 5 
x ′′ (t) . x(t)+ y(t)+ z(t)=0 
x(t)+ y(t)+ z ′′ (t) . z(t)=0
　　解考虑给出的初值，将方程中的 x ′′替换成 s2X(s) . sx(0) = s2X(s) . s，并将 x、 
).+++5213132722-15　　例 试求解下面的高阶常系数线性微分方程组。 ′′ .()+ ()()+ ()=0 ttttxyyz′′′其中，。(0)=1(0)= (0)= (0)= (0)= (0)=0xyzxyz， .下面的代数方程。 2...X() X()+ Y ()+ Z()=0 ssssss. 2.X()+ Y () Y ()+ Z()=0 sssss..2.X()+ Y ()+ Z() Z()=0sssss可以使用下面的语句求解该代数方程。 eqns=[s^2*X-s-X+Y+Z==0,X+s^2*Y-Y+Z==0,X+Y+s^2*Z-Z==0]; sol=solve(eqns,[X,Y,Z]);sol.X,sol.Y,sol.ZLaplace得到的变换表达式如下： ss .X()= Y ()= Z()= sss, 
y和 z替换成 X(s)、Y (s)和 Z(s)，将 y ′′和 z ′′替换成 s2Y (s)和 s2Z(s)，可以直接列写出
(	XYZ >>symss 
. .. .. 
3 
(s2 + 1)(s2 . 2) 	(s2 + 1)(s2 . 2)
对得到的代数方程的解进行 Laplace反变换： 
>> 	x=ilaplace(sol.X), y=ilaplace(sol.Y) 
z=ilaplace(sol.Z); fplot([x,y],[0,0.5])

可以得出微分方程的解析解为 
∨	∨
21 	11 
x(t)= cosh 2t + cos t, y(t)= z(t)= cos t . cosh 2t
33 33 由微分方程的解还可以绘制出信号 x(t)和 y(t)的时间响应曲线，如图 2-9所示。 

图 2-9方程解的时间响应曲线


2.4一般微分方程的解析解
本节介绍使用 MATLAB语言及其符号运算工具箱求解常系数线性微分方程的解析解的方法。
前面介绍了低阶微分方程的底层求解方法，这些方法的局限性较大，不适合求解复杂的问题，不同的微分方程需要引入不同的求解方法，不利于微分方程问题的实际求解，所以需要统一的、更规范的求解方法。本节介绍基于 MATLAB符号运算工具箱的专门函数求解一般微分方程解析解的方法。 
2.4.1简单微分方程的解析解 
MATLAB符号运算工具箱提供了 dsolve()函数，可尝试求出一般微分方程
的解析解。如果某个微分方程被证实没有解析解，则应该考虑获得其数值解。 MATLAB的 dsolve()函数可以采用下面的形式调用： sols=dsolve(f1,f2,··· ,fm); %默认的自变量为 t y=dsolve(f1,f2,··· ,fm,varlist); %指明自变量 [x,y,··· ]=dsolve(f1,f2,··· ,fm,varlist); %指明返回函数的变量名 MATLAB的早期版本中，允许使用符号表达式和字符串描述各个微分方程 fi，
而在以后的版本中，有可能只支持符号表达式的描述方式。所以，在本书中，将减少介绍字符串的描述方式，侧重于符号表达式的描述与求解方式。
　　例 2-16试利用计算机工具求出式（1-2-1）中的微分方程。　　解有了强大的计算机工具，当时 Newton无法全面研究的问题就可用下面简单、规范的 MATLAB命令直接解出。 
>> 	syms x y(x) %申明符号变量，并申明符号函数 
y1=dsolve(diff(y)==1-3*x+y+x^2+x*y)

得出的解析解为 	∨ (∨ ())
y1(x)= Cex(x+2)/2 . x +3 2丘 erf 2 x 2 +21 e(x+1)2/2 +4
其中， C为任意常数， erf(·)是一个特殊函数。注意，在描述方程中的等号时，应该使用双等号 ==，否则函数调用语句会报错。实际上， Newton得到的解只是 y(0) = 0时的特解。用 Taylor幂级数展开该解： 
>> 	y2=dsolve(diff(y)==1-3*x+y+x^2+x*y,'y(0)==0'); 
y2=expand(taylor(y2,'Order',10))

得到下面的表达式，其中末尾 6项与 Newton的结果一致。 
9 8 7 6543
x13 xxxxxx
2 
y2(x) .. . + . + . + . x + x
9072 5040 630 45 30 6 3 
　　例 2-17试求解式（1-2-2）中的微分方程。
·28·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
　　解尝试利用下面的 MATLAB语句直接求解，但是将得出方程无解的提示。 
>> 	syms x y(x); syms a positive 
y1=dsolve(diff(y)==-y/sqrt(a^2-y^2)) %直接求解微分方程

为什么方程无解析解，而 Leibniz得出了解析解呢？分析 Leibniz的解可知，他得出的是以 y为自变量的，关于 x的显式解，而反过来确实求不出 y = f(x)式的解析解。仿照 Leibniz的思路，利用 dy/dx = 1/(dx/dy)这一事实，可以使用下面的语句： 
>> syms y x(y); x1=dsolve(1/diff(x)==-y/sqrt(a^2-y^2)) %微分方程的解析解得出的解析解为 
(√ )
√ 
x1(y)= C . a ln a2 . 1 . a 1 . √ a2 . y2
y2 y2 
可以看出，Leibniz当时得出的只是方程的一个特解。
　　例 2-18试用微分方程求解函数直接求解例 2-4中给出的一阶微分方程。　　解采用 dsolve()函数求解微分方程，需要申明符号变量 x与符号函数 y(x)，可以用符号表达式直接描述该微分方程，并直接求解。注意，在描述微分方程时，如果右侧的表达式为零，则可以略去 ==0。由下面的语句可以描述并求解微分方程，得到方程的解析解。 
>> 	syms x y(x) 
eqn=diff(y)+8*y+y^2==-15; y=dsolve(eqn, y(0)==0) 
y=simplify(y) %得出微分方程的解析解并化简

得到微分方程的解 y(x)= .(15e2x . 15)/(5e2x . 3)。实际上，这里得到的结果与例 2-4的结果虽然形式上不一致，但二者是相等的。
在 MATLAB的早期版本中，还允许使用字符串描述微分方程与初值。例如，采用下面的语句直接求解问题，可得出完全一致的结果。使用字符串描述微分方程时，通常需要使用记号 D3y描述 y的三阶导数。本书不建议使用字符串描述微分方程。 
>> 	y=dsolve('Dy+8*y+y^2=-15','y(0)=0'); y=simplify(y)

　　例 2-19试求解下面的线性微分方程组。 
{ 
′ 
y1(t)= y2(t) 
′ 
y2(t)= .λy1(t)
其中，微分方程的初值为 y1(0) = a，y2(0) = b，且 λ> 0。　　解首先申明正的符号变量 λ，可以使用下面的语句描述微分方程与已知初值，并直接求解该微分方程。 
>> syms t a b y1(t) y2(t); syms lam positive 
[y1,y2]=dsolve(diff(y1)==y2,diff(y2)==-lam*y1,... 
y1(0)==a,y2(0)==b) %微分方程组的求解

得出微分方程组的解析解为 
∨ ∨ 	∨∨∨ 
y1(t)= a cos λt + ∨ b sin λt, y2(t)= b cos λt . aλ sin λt λ 

2.4.2常系数高阶线性微分方程的解析解
前面介绍的 dsolve()函数可以直接用于高阶线性微分方程或其他复杂微分方程的直接求解。只需将微分方程及已知条件用符号表达式直接描述出来，就可以调用该函数直接求解微分方程了。
为简洁地描述微分方程与已知条件，可定义一些中间变量记录 y(t)信号的各阶导数，中间变量的定义方法将通过下面的例子演示。
　　例 2-20试重新求解例 2-14中的微分方程，并与前面的结果进行比较。　　解由于需要描述 y(x)函数的各阶导数及其初值，所以定义一些中间变量描述 y(x)的各阶导数。例如，可以用变量 d3y描述中间变量 y ′′′ (x)。从规范求解的角度看，采用这样的记号会使得方程描述更简洁易懂。用下面的语句可以直接求解该微分方程，得出的结果与例 2-14中得到的结果完全一致。可以看出，这里给出的语句更简洁易懂，更容易检验描述是否正确，便于查错。最后的绘图语句可以得到与图 2-8完全一致的结果。 
>> 	syms t y(t); u=exp(-5*t)*cos(2*t+1)+5; 
d1y=diff(y); d2y=diff(y,2); d3y=diff(y,3); d4y=diff(y,4); 
y=dsolve(diff(d4y)+12*d4y+57*d3y+134*d2y+156*d1y+72*y==u,... 

y(0)==3,d1y(0)==2,d2y(0)==0,d3y(0)==0,d4y(0)==0); y=simplify(y), fplot(y,[0,10]) %微分方程的求解与绘图
　　例 2-21假设输入信号为 u(t)= e.5t cos(2t + 1) + 5，试求出下面微分方程的通解。 
y(4)(t) + 10y ′′′ (t) + 35y ′′ (t) + 50y ′ (t) + 24y(t)=5u ′′ (t)+4u ′ (t)+2u(t)
　　解由于不涉及信号 y(t)的初值，所以可以按照例 2-20的方式申明一些中间变量，也可以不申明这些中间变量，直接用 diff()函数表示导数。 >> syms t y(t); u(t)=exp(-5*t)*cos(2*t+1)+5; d1y=diff(y); d2y=diff(y,2); d3y=diff(y,3); d4y=diff(y,4); y=dsolve(d4y+10*d3y+35*d2y+50*d1y+24*y==... 5*diff(u,t,2)+4*diff(u,t)+2*u); y=simplify(y) %求微分方程的解析解并化简结果通过上述语句进行求解，可以得出原微分方程的通解为 
5 343 	547 
y(t)= . e.5t cos(2t + 1) . e.5t sin(2t + 1) 
12 520 	520 
+C1e.4t + C2e.3t + C3e.2t + C4e.t其中， Ci为任意常数。若给出初值或边值，则可以通过已知条件建立代数方程，求出 Ci的值。这样的思路和高等数学中微分方程的求解是一致的。
·30·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
检验上述微分方程的通解，可以将其代入原方程（有时需要调用 simplify()函数化简得出的误差），可见误差为零，说明无论 Ci如何取值，得出的通解都满足原始方程。 
>> err=diff(y,4)+10*diff(y,3)+35*diff(y,2)+50*diff(y)+24*y-... (5*diff(u,t,2)+4*diff(u,t)+2*u) %解的检验
　　例2-22前面介绍的方程只含有实数极点，实际上符号运算工具箱提供的 dsolve()函数同样适用于有复数极点的微分方程的解析解。假设微分方程如下： 
y(5)(t)+5y(4)(t) + 12y ′′′ (t) + 16y ′′ (t) + 12y ′ (t)+4y(t)= u ′ (t)+3u(t)
假设输入信号为正弦信号 u(t)= sin t，且 y(0)= y ′ (0)= y ′′ (0)=y ′′′ (0)=y(4)(0)=0，试用解析方法求解该方程。　　解用下面的语句可以求出原微分方程的解析解。首先，定义一些中间变量描述并求解微分方程，得出的结果如图 2-10所示。 
>> 	syms t y(t); u(t)=sin(t); 
d1y=diff(y); d2y=diff(d1y); d3y=diff(d2y); d4y=diff(d3y); 
y=dsolve(diff(d4y)+5*d4y+12*d3y+16*d2y+12*d1y+4*y==... 

diff(u)+3*u,... y(0)==0,d1y(0)==0,d2y(0)==0,d3y(0)==0,d4y(0)==0); y=simplify(y), fplot(y,[0,30]) %微分方程的求解与绘图

图 2-10趋于等幅振荡的曲线
方程解析解的数学描述为 
124 111 
y(t)= e.t . cos t . sin t . e.t cos t + e.t sin t . te.t cos t
555 102 
或更简单地手工修改为 
()
12 4 t 11 
y(t)= e.t . cos t . sin t . + e.t cos t + e.t sin t
5552 10 
由得出的曲线可见，该解在 t值较大时近似等幅振荡的曲线。因为其解析解的前两项是等幅振荡的正余弦函数，当 t逐渐增大，后面各项逐渐消失。
　　例 2-23试重新求解例 2-15中给出的常系数线性微分方程组。
　　解可以用下面的命令直接求解该微分方程组，得出的结果与例 2-15完全一致。 
>> 	syms t x(t) y(t) z(t) 
d1x=diff(x); d1y=diff(y); d1z=diff(z); 
[x,y,z]=dsolve(diff(x,2)-x+y+z==0,x+diff(y,2)-y+z==0,... 

x+y+diff(z,2)-z==0, x(0)==1, d1x(0)==0,... 
y(0)==0, d1y(0)==0,z(0)==0, d1z(0)==0); 
x=simplify(x), y=simplify(y), z=simplify(z)

　　例 2-24试求解线性微分方程组。 
{ 
x ′′ (t)+2x ′ (t)= x(t)+2y(t) . e.t y ′ (t)=4x(t)+3y(t)+4e.t
　　解线性微分方程组也可以用 dsolve()函数直接求解。上述线性微分方程组可以
由下面的 MATLAB语句直接求解。 >> syms t x(t) y(t) %声明符号变量与函数，下面的语句直接求解微分方程 [x,y]=dsolve(diff(x,2)+2*diff(x)==x+2*y-exp(-t),... diff(y)==4*x+3*y+4*exp(-t)) %微分方程求解
得出的结果为 
{ 	叫叫 
+ C2e(1+6)t + C3e.(.1+6)t x(t)= .6te.t + C1e.t (∨ )(1+叫 6)t (∨ )C3e.(.1+叫 6)t y(t)=6te.t .C1e.t +22+ 6C2e+22. 6+ e.t/2 

2.4.3线性时变微分方程的解析解
前面研究的是常系数线性微分方程的直接求解方法，如果线性微分方程的系
数是自变量的函数，则可以将自变量看作时间变量，这类微分方程又称为时变微分
方程。时变微分方程是可以尝试由 dsolve()函数直接求解的，下面给出几个例子，
演示时变微分方程的直接求解方法。
　　例 2-25考虑下面给出的二阶时变线性微分方程。 
y ′′ (x)+ ay ′ (x)+(bx + c)y(x)=0
　　解将微分方程用符号表达式表示出来，使用下面的求解语句： 
>> 	syms x a b c y(x) 
y=dsolve(diff(y,2)+a*diff(y)+(b*x+c)*y==0)

可以得出微分方程的解析解如下所示，其中用到了 Airy特殊函数。 
(	)()
C1 .a2 +4c +4bx C2 .a2 +4c +4bx 
y(x)= ∨ Ai . 	+ ∨ Bi . 
eax 	1/3 eax 1/3
4b(1/b)	4b(1/b)
注意，对这个具体问题而言，不宜使用下面的字符串描述微分方程。 
>> 	y=dsolve('D2y+a*Dy+(b*x+c)*y') 

·32·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
2-26　　例 试求出时变线性微分方程的解析解。 　　解dsolve()对于某些时变线性微分方程，仍然可以尝试用函数直接求解。 y=simplify(y)%微分方程的求解与结果的化简处理)′2.(21)()+2()=0xxyxyxdsolve()中间变量，直接描述微分方程即可。可以由函数直接求解时变微分方程。 y=dsolve(x^2*(2*x-1)*diff(y,3)+... -2*x*diff(y)+2*y)%微分方程解的检验
√ 
>> syms x y(x) 

33 3 
y(x)= C1x + . 2C2 x ++ C3 (2x + 6) x + 
√ 222
　　例 2-27试求解时变微分方程。 
′′′ (x) . (4x . 3)xy ′′ (x) . 2xy 
>> syms x y(x) 

D2y用这样的记号描述 函数的二阶导数时，计算机将误以为自变量为 而不是实ty'x'际的自变量得出的解析解是错误的，应该在语句末尾加选项。x，′′′′3 .(2+3)()+3(2+3)()6()=0xyxxyxyxy=dsolve((2*x+3)^3*diff(y,3)+3*(2*x+3)*diff(y)-6*y==0); 0可以得出微分方程的解析解如下所示，将解析解代入原方程可知误差为。 (　　解如果把自变量 看成时间，则微分方程的系数都是时间的函数，所以这类微分x方程又可以认为是时变微分方程。对这个具体问题而言，无须描述初值，所以不必建立(4*x-3)*x*diff(y,2)-2*x*diff(y)+2*y==0); simplify(y)%微分方程的解析求解与化简 simplify(x^2*(2*x-1)*diff(y,3)+(4*x-3)*x*diff(y,2)... 0得出该方程的解析解如下所示，代入原方程误差为 。由于 是自由常数，所以 .C2C11 可以合并成一个常数则可以得出更简单的结果。CC，30
)
(
y(x)= . 12C1 . C3 +8C3x + 32C2x 2 +8C3x 2 ln x
16x
2.4.4线性时变微分方程组的求解
. .
. 
线性时变微分方程组也可以尝试由符号运算工具箱中的 dsolve()函数直接求解，如果能得出解析解更好，如果得不到解析解也正常，毕竟存在解析解的时变微分方程组不是很多。如果不能得出问题的解析解，可以考虑获得问题的数值解。下面给出一个微分方程组的解析解方法。
　　例 2-28试求出下面线性时变微分方程组的解析解。 
d2y(x) 
+2y(x)+4z(x)= ex 
dx2 
d2z(x) 
. y(x) . 3z(x)= .x
dx2 
.
. 
　　解求解之前需要申明 y和 z均为 x的函数。由下面的语句可以直接求解本例的线性微分方程组。 
>> 	syms x y(x) z(x) 
[y1,z1]=dsolve(diff(y,2)+2*y+4*z==exp(x),... 

diff(z,2)-y-3*z==-x) y1=simplify(y1), z1=simplify(z1) %微分方程组求解与化简 simplify(diff(y1,x,2)+2*y1+4*z1-exp(x)) simplify(diff(z1,x,2)-y1-3*z1+x) %微分方程解的检验
得出方程组的解如下所示，经检验该通解满足原始方程组。 
y1(x)= ex . 2x . 4C1 cos x . 4C2 sin x .∨ 2 C3e叫 2x + ∨ 2 C4e.叫 2x 22 z1(x)= x . ex + C1 cos x + C2 sin x + ∨ 2 C3e叫 2x .∨ 2 C4e.叫 2x 2 22 
. .
. .
. 

2.4.5边值问题的计算机求解
前面介绍的 dsolve()函数还可以用于求解边值问题。下面通过例子演示微分方程边值问题的求解方法。
　　例 2-29考虑例 2-20中给出的微分方程，假设已知边界处的函数值为 y(0) = 1/2， y ′ (0) = 1，y(2丘)=0，y ′ (2丘)=0，试求解该微分方程。　　解回顾 2.3节介绍的基于 Laplace变换的方法。Laplace变换需要使用函数的初值，而这里涉及 y(x)函数在不同时刻的值，所以该微分方程不能使用 Laplace变换方法进行求解，只能借助于 dsolve()函数求解该微分方程。下面给出的语句都是常规的，采用求解一般微分方程的统一框架。 
>> 	syms t y(t); u(t)=exp(-5*t)*cos(2*t+1)+5; 

d1y=diff(y); d2y=diff(y,2); d3y=diff(y,3); d4y=diff(y,4); 

y=dsolve(d4y+10*d3y+35*d2y+50*d1y+24*y==... 

5*diff(u,t,2)+4*diff(u,t)+2*u,... 

y(0)==1,d1y(0)==0,y(2*pi)==0,d1y(2*pi)==0); 

y=simplify(y), fplot(y,[0,10]) %微分方程的求解与绘图可以得出满足该条件的方程解析解。由于该解过于冗长，这里暂不给出。该解对应的函数曲线如图 2-11所示。 
dsolve()函数还可以用于直接求解更复杂的已知点问题。只要相互独立的已知条件个数等于方程解析解的待定系数个数，就可以调用该函数直接求解，得到微分方程的待定系数，从而得出微分方程的解析解。
　　例 2-30考虑例 2-20中给出的微分方程，假设已知一些点的函数值：y(0) = 1/2， y ′ (丘)=0，y ′′ (2丘)= y ′ (2丘)=0，试求解该微分方程。
·34·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解

图 2-11方程解的曲线
　　解和例 2-29相比，这个例子给出了 t =丘这个中间点的已知值，其他的没有变化，所以可以直接采用类似的语句求解方程。 
>> 	syms t y(t); u(t)=exp(-5*t)*cos(2*t+1)+5; 
d1y=diff(y); d2y=diff(y,2); d3y=diff(y,3); d4y=diff(y,4); 
y=dsolve(d4y+10*d3y+35*d2y+50*d1y+24*y==... 

5*diff(u,t,2)+4*diff(u,t)+2*u,... y(0)==1/2,d1y(pi)==0,d2y(2*pi)==0,d1y(2*pi)==0); y=simplify(y), fplot(y,[0,2*pi]) %求解微分方程并绘图该解对应的函数曲线如图 2-12所示。虽然已知若干零初值，该解析解的表达式仍过于复杂，此处不给出。

图 2-12方程解的曲线 


2.5线性矩阵微分方程的求解
高阶常系数线性微分方程也可以表示成线性状态空间方程的形式。本节将探讨一般线性状态空间方程的数学形式并给出闭式解的数学公式，然后介绍直接求解方法。此外，本节还将探讨 Sylvester微分方程的解析求解方法。
2.5.1线性状态空间方程的解析解
线性状态空间模型是很多领域经常使用的微分方程数学模型，是前面介绍的常系数线性高阶微分方程的另一种表现形式。例如，控制科学领域通常采用线性状态空间模型描述一般线性系统动态模型。本节首先给出状态空间模型的数学形式，然后给出线性状态方程的解析解方法。
　　定义 2-18假设线性状态空间模型的一般表示为 
{ 
x ′ (t)= Ax(t)+ Bu(t) （2-5-1）
y(t)= Cx(t)+ Du(t)其中， A, B, C, D为常数矩阵，向量 x(t)称为系统的状态变量向量，且已知状态向量初值 x(t0)= x0。u(t)与 y(t)可以分别认为是系统的输入与输出信号。式（ 2-5-1）中第一个式子可以改写成 x ′ (t) . Ax(t)= Bu(t)（2-5-2）该式左右两端同时乘以 e.At，可得 
e.At x ′ (t) . e.AtAx(t)= e.AtBu(t)（2-5-3）
由于 Ae.At = e.AtA，所以等式左边正好是导数公式。 
dd t[e.At x(t)] = e.AtBu(t)（2-5-4）
对方程两边进行定积分运算，可以推导出下面的公式。 
e.At x(t) . e.At0 x(t0)= ∫t e.Aλ Bu(ξ)dξ（2-5-5）
t0 
由该式可以直接推导出线性微分方程的解析解定理。　　定理 2-12定义 2-18中状态方程的解析解可以写成 
x(t)= eA(t.t0)x(t0)+ ∫t eA(t.λ)Bu(ξ) dξ（2-5-6）
t0 
其中， eA(t.t0)又称为状态转移矩阵，记作 i(t, t0)= eA(t.t0)，用于描述从 t0到 t时刻的状态变迁。
对于给定的输入信号 u(t)，通过矩阵积分运算可以直接得出该方程的解析解。从解析解表达式可见，涉及矩阵的 e指数求解，也涉及函数的积分运算，可以调用 MATLAB的符号运算函数进行计算。利用 MATLAB符号运算工具箱中的函数，式（ 2-5-6）可以由下面的语句直接求解。 
x=expm(A*(t.t0))*x(t0)+int(expm(A*(t.ξ))*B*subs(u,t,ξ),ξ,t0,t) 
·36·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
　　定义2-19如果定义 2-18中的外部输入信号 u(t)三 0，则系统在状态变量初值 x(t0)的驱动下运动。这时，系统的状态方程模型为 
x ′ (t)= Ax(t), x(t0)= x0（2-5-7）　　例 2-31假设输入信号为 u(t)=2+2e.3t sin 2t，试求出下面矩阵描述的状态空
间方程的解析解。  
  .19  .16  .16  .19  .  1  .  2  .  
A =..  21 20  16 17  17 16  19 20  ..,  B =.. 0 1 ..,  CT =.. 1 0 ..,  D = 0 
.20  .16  .16  .19  2  0  

其中，初始状态向量为 x = [0, 1, 1, 2]T。　　解将相关矩阵与向量输入 MATLAB工作空间，对照式（2-5-6），可以直接得出方程的解析解。 
>> 	syms t tau; u(t)=2+2*exp(-3*t)*sin(2*t); %定义输入信号 A=[-19,-16,-16,-19; 21,16,17,19; 20,17,16,20; -20,-16,-16,-19]; 
)
B=[1; 0; 1; 2]; C=[2 1 0 0]; x0=[0; 1; 1; 2]; y=C*(expm(A*t)*x0+int(expm(A*(t-tau))*B*u(tau),tau,0,t)); simplify(y) %求解式（2-5-6）的直接实现
得出的解析解的数学表示为 
119 127t 135 77 
y(t)= e.t + 57e.3t + e.t +4t2e.t . e.3t cos 2t + e.3t sin 2t . 54
8484 (
手工合并同类项，该方程的解析解可简化为 
119 127t 	135 77
2
y(t)=++4te.t + 57e.3t . e.3t cos 2t + e.3t sin 2t . 54 
84 	84 

2.5.2状态方程的直接求解
考虑利用 dsolve()函数直接求解状态方程，由前面的叙述可知，这样做的难点在于如何构造一个符号型函数矩阵 X(t)。作者编写了 any_matrix()函数，用于直接定义任意符号函数矩阵，该函数的程序如下所示，后面将通过例子演示该函数的使用方法。 
function A=any_matrix(nn,sA,varargin) %生成任意矩阵函数 v=varargin; n=nn(1); if length(nn)==1, m=n; else, m=nn(2); end s=''; k=length(v); K=0; if n==1 || m==1, K=1; end if k>0, s='('; for i=1:k, s=[s ',' char(v{i})]; end 
s(2)=[]; s=[s ')']; end for i=1:n, for j=1:m %采用循环结构对元素逐个处理
if K==0, str=[sA int2str(i),int2str(j) s]; 
else, str=[sA int2str(i*j) s]; end 
eval(['syms ' str]); eval(['A(i,j)=' str ';']); 

end, end

　　例 2-32试输入如下两个符号函数矩阵。 
. . ..
a11(x, y)  a12(x, y)  a13(x, y)  a14(x, y)  f11(t)  f12(t)  
A = .. a21(x, y) a31(x, y)  a22(x, y) a32(x, y)  a23(x, y) a33(x, y)  a24(x, y) a34(x, y)..,  B = .. f21(t) f31(t)  f22(t) f32(t)  
a41(x, y)  a42(x, y)  a43(x, y)  a44(x, y)  f41(t)  f42(t)  

　　解使用下面的语句可自动生成两个符号函数矩阵，其矩阵形式与预期完全一致。 
>> 	syms x y t 
A=any_matrix(4,'a',x,y), B=any_matrix([4,2],'f',t)

　　例 2-33试用微分方程描述的方法重新求解例 2-31中给出的状态方程。　　解有了 any_matrix()函数，可以用符号表达式描述原始微分方程，再直接求解，得出与前面完全一致的结果。 
>> 	syms t; u(t)=2+2*exp(-3*t)*sin(2*t); A=[-19,-16,-16,-19; 21,16,17,19; 20,17,16,20; -20,-16,-16,-19]; B=[1; 0; 1; 2]; C=[2 1 0 0]; x0=[0; 1; 1; 2]; x(t)=any_matrix([4,1],'x',t); %构造任意矩阵函数 y=dsolve(diff(x)==A*x+B*u,x(0)==x0); %求线性微分方程的解析解 y=C*[y.x1; y.x2; y.x3; y.x4] %构造输出信号解析解 
2.5.3 Sylvester微分方程的求解 
Sylvester微分方程是一类特殊形式的矩阵型微分方程。本节首先给出该方程的一般定义，然后探讨该方程的三种求解方法。
　　定义 2-20矩阵 Sylvester微分方程的一般形式为 [14] 
X ′ (t)= AX(t)+ X(t)B, X(0) = X0（2-5-8）
其中，A → Rn×n，B → Rm×m，X, X0 → Rn×m。
　　定理 2-13矩阵 Sylvester微分方程的解 [14]为 X(t)= eAtX0eBt。
　　例 2-34试求解下面给出的矩阵微分方程。 
..
]
[
. 	..
.. 
.1 .20 .10 .1 .1 .3 .1 .2 .21 11 
X ′ (t)=	X(t)+ X(t), X(0) =
.11 .20 0 .210 1211 01 
　　解将相关矩阵输入 MATLAB环境中，由定理 2-13可直接求解该微分方程。
. ..
·38·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 
>> 	A=[-1,-2,0,-1; -1,-3,-1,-2; -1,1,-2,0; 1,2,1,1]; 
B=[-2,1; 0,-2]; X0=[0,-1; 1,1; 1,0; 0,1]; 
A=sym(A); B=sym(B); X0=sym(X0); syms t 

到的误差矩阵都是零矩阵。 
(t2/2 . 2)e.3t +2e.4t (2t + 1)e.4t +(.2+ t2 + t3 . 4t)e.3t 2e.4t . (t + 1)e.3t (2t + 1)e.4t . (t2 +3t)e.3t X(t)= 
.e.3t(t2 . 2)/2 .te.3t(t2 +2t . 6)/2 (2 + t)e.3t . 2e.4t (2 + t2 +4t)e.3t . (2t + 1)e.4t 
示该方程的直接解析求解方法。
　　例 2-35试用 dsolve()函数直接求解例 2-34中的 Sylvester微分方程。　　解与前面的例子类似，先将 X定义为符号矩阵函数，然后采用下面的求解语句求解原方程，得出的结果与例 2-34完全一致。 
>> 	A=[-1,-2,0,-1; -1,-3,-1,-2; -1,1,-2,0; 1,2,1,1]; 

X=simplify(expm(A*t)*X0*expm(B*t))%求解微分方程 	simplify(diff(X)-A*X-X*B),subs(X,t,0)-X0%检验得到的解如下所示。将微分方程的解代回原方程，可以发现该解满足原方程与初值，得.... Sylvesterdsolve()也可以考虑利用 函数直接求解 微分方程。下面通过例子演A=sym(A);B=sym(B);X0=sym(X0);syms t 
B=[-2,1; 0,-2]; X0=[0,-1; 1,1; 1,0; 0,1]; 

X(t)=any_matrix([4,2],'x',t); %构造任意矩阵函数 x=dsolve(diff(X)==A*X+X*B, X(0)==X0) %求解线性微分方程 X=[x.x11, x.x12; x.x21 x.x22; x.x31 x.x32; x.x41 x.x42] 

2.5.4基于 Kronecker乘积的 Sylvester微分方程直接求解
考虑引入本丛书卷 III介绍的 Kronecker乘积，将矩阵型 Sylvester方程变换成向量型矩阵方程。引入 Kronecker乘积，则式（ 2-5-8）可变换为 
. ....
)
(
x ′ (t)=Im @ A + BT @ Inx(t)（2-5-9）其中， x(t)= vec(X(t))为矩阵 X(t)按列展开得到的列向量。将方程变换为一阶显
式微分方程的标准型，可采用 dsolve()函数直接求解，也可以通过 expm()函数求解。求解之后再由 reshape()函数将向量变换回矩阵形式。
　　例 2-36试利用 Kronecker乘积公式重新求解例 2-34中的 Sylvester微分方程。　　解利用 Kronecker乘积得出系数矩阵，由下面的语句可得出原 Sylvester微分方程的所有解，与例 2-34的结果完全一致。 
>> 	A=[-1,-2,0,-1; -1,-3,-1,-2; -1,1,-2,0; 1,2,1,1]; 
B=[-2,1; 0,-2]; X0=[0,-1; 1,1; 1,0; 0,1]; 

A0=kron(eye(2),A)+kron(B',eye(4)); syms t 

x=expm(A0*t)*X0(:); x1=reshape(x,4,2) %求解矩阵微分方程
若采用微分方程求解语句，则可以使用下面的语句并得出完全一致的结果。 
>> 	x(t)=any_matrix([8,1],'x',t); 
y=dsolve(diff(x)==A0*x, x(0)==X0(:)); %另一种求解方法 
x2=[y.x1 y.x5; y.x2 y.x6; y.x3 y.x7; y.x4 y.x8] 




2.6特殊非线性微分方程的解析解
部分非线性微分方程也是可以用 dsolve()函数求解析解的，这类方程描述方式和前面介绍的线性微分方程是一致的。下面通过例子演示非线性方程的解析解求解方法，同时给出不能求解的例子。 
2.6.1可解的非线性微分方程
一般情况下无法预先判断一个给定的非线性常微分方程是否解析可解，可尝试将微分方程以规范的格式输入 MATLAB工作空间，并利用 dsolve()函数求解。非线性方程的可解性只能通过尝试求解判定。下面通过例子演示几个非线性微分方程的求解方法。
　　例 2-37前面介绍的低阶方程很多是非线性方程。例如，例 2-4中给出了一阶非线性微分方程，试求解该微分方程。 
dy(t) 
+8y(t)+ y 2(t)= .15,y(0) = 0
dt 
　　解可以用前面介绍的常规方法描述并求解微分方程，得出的解为 Y = .(15e2t . 15)/(5e2t . 3)，与前面给出的是等效的。可以看出，这里介绍的方法更直接，无须任何技巧，只要规范地将微分方程用指定的格式描述出来，再调用求解函数，即可得出原始非线性微分方程的解。 
>> 	syms t y(t) 
Y=dsolve(diff(y)+8*y+y^2==-15, y(0)==0); 
simplify(Y) %进行求解与化简

　　例 2-38试求出一阶非线性微分方程 x ′ (t)= x(t)(1 . x2(t))的解析解。　　解这类一阶非线性方程可以考虑用 dsolve()函数直接求解。 
>> syms t x(t); x=dsolve(diff(x)==x*(1-x^2)) %非线性方程的直接求解该微分方程的解析解为 x(t)= √ .1/ (eC.2t . 1)
此外，常数 ±1与 0均为方程的解。
　　例 2-39试求解下面的二阶非线性微分方程 [12]。 
(d2y(x) )2 ()d2y(x) (dy(x))2 dy(x)
+	2ay(x)+ bx + c. a . b + k =0 
dx2 	dx2 dx dx 
·40·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
　　解这个非线性微分方程看起来比较复杂，对非专业人士而言，很难手工求解。可以考虑用 MATLAB将其表示出来，再通过求解语句求解。 
>> 	syms a b c k x y(x) 
d1y=diff(y); d2y=diff(y,2); 
y=dsolve(d2y^2+(2*a*y+b*x+c)*d2y-a*d1y^2-b*d1y+k==0)

得出问题的解如下所示。可以看出，这里得到的两个解的分支都是正确的。 
C1x2 1 
+ C2x . (C12 + cC1 . aC22 . bC2 + k) 
y(x)= C32 e.叫 ax . 1 (2cC+ 1abx)+ 1 e叫.ax(b2 +4ak)2a 16C3a3 
. .. .. 
　　例 2-40试求解下面给出的三阶非线性微分方程在 t → (0.2,丘)时的解析解 [12]，并绘制解的曲线。 
5 
xy ′′′ (x) = 2(xy ′ (x) . 2y(x)),y(1) = 1,y ′ (1) = 0.5,y ′′ (1) = .1
　　解与一阶、二阶非线性微分方程相比，存在解析解的三阶微分方程更少。实际上，无论哪类微分方程，都可以尝试使用常规方法求解。对这里给出的微分方程与初值条件而言，可直接定义中间变量并求解微分方程。如果能得到方程的解，则可以绘制出方程解的曲线，如图 2-13所示。 
>> 	syms x y(x) 

d1y=diff(y); d2y=diff(d1y); d3y=diff(d2y); 

y=dsolve(x^5*d3y==2*(x*d1y-2*y),... 

y(1)==1, d1y(1)==0.5, d2y(1)==-1) fplot(y,[0.2,pi]) %求解微分方程并绘图得出方程的解析解为 
2 .∨ 叫 22e.叫 ∨ 2e.叫 22e叫
3	2e3
2/x 	2/x
y(x)= xx + x 
8	8 

图 2-13方程解曲线

2.6.2解析解不存在的非线性微分方程
从前面的例子可以看出， MATLAB的微分方程求解函数 dsolve()特别强大，只需描述微分方程就可以直接求解。事实上并非如此，前面的例子中特意选择了有解析解的非线性微分方程，而非任意的微分方程，绝大多数非线性微分方程是没有解析解的，下面给出几个例子。
　　例 2-41试求出一阶非线性微分方程 x ′ (t)= x(t)(1 . x2(t)) + 1的解析解。　　解该微分方程只是对例 2-38作了微小变化，即等号右侧加上 1。可尝试使用下面的语句求解该方程。但是执行不成功，并得到警告信息 “Warning: Unable to find explicit solution. Returning implicit solution instead”（警告：无法找到显式解，返回隐式解），说明该方程的解析解是不存在的。 
>> 	syms t x(t); 
x=dsolve(diff(x)==x*(1-x^2)+1) %方程是无解的

　　例 2-42荷兰物理学家 Balthasar van der Pol（1889.1959）与其同事在 1927年报道了一个振荡电路 [15]，该电路对应的微分方程通常被称为 van der Pol方程，其数学表达式如下： d2y(t) []dy(t) 
dt2 + μy 2(t) . 1dt + y(t)=0（2-6-1）
试用 dsolve()函数求解该方程，看看能得出什么结论。　　解对该 van der Pol方程，尝试如下的 MATLAB命令，返回警告信息 “Unable to find explicit solution”（无法找到显式解），表明求解不成功。 
>> 	syms t y(t) mu; 
y=dsolve(diff(y,2)+mu*(y^2-1)*diff(y)+y==0) %直接求解

可见，微分方程解析解求解函数 dsolve()并不能直接应用于一般非线性方程解析解的求解。因此，非线性微分方程只能使用数值解法求解。即使看起来很简单的非线性微分方程也可能是没有解析解的，只有极特殊的非线性微分方程解析可解。后续各章内容将集中介绍各类非线性微分方程的数值解方法。
本章习题 
2.1	试求出下面一阶微分方程的解 [12]。
√ 
x(x + 1) 
（1）y ′ (x)= ∨ x + ∨ 1+ x
()
t . sin 2t 
（2）y ′ (t)= t2 + cos 2ty(t),y(0) = 4

（3）y ′ (t)+ (1+ t23/4)3/2 y(t)= (1+ 2tt 2/4)2 ,y(0) = 8。


·42·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 
2.
2试证明超几何函数满足下面的条件 [16]。

（1）1F1(a; a; x)= ex

（2）2F1(a, 1; 1; z)= (1 . 1 z)a

（3）2F1(1, 1; 2; z)= z 1 ln(z + 1)

（4）2F1(1/2, 1; 3/2; z 2)= 1 lnz +1 

2z 1 . z

（5）2F1(1/2, 1/2; 3/2; z 2)= 1 arcsinz 



2.3试利用 
Laplace变换求出下分方程的解析解，且已知初值 x(0) = 2，y(0) = 1， x ′ (0) = y ′ (0) = .1。 
面微z

{ 
x ′′ (t)+ y ′′ (t)+ x(t)+ y(t)=0 ′′ (t) . y
2x ′′ (t) . x(t)+ y(t)= sin t 
2.4试求出下面线性微分方程的通解。 
y(5) + 13y(4)+ 64y ′′′ (t) + 152y ′′ (t) + 176y ′ (t) + 80y(t) 
= e.2t [sin (2t +丘 ) + cos 3t]3 假设上述微分方程满足已知条件 y(0) = 1，y(1) = 3，y(丘)=2，y ′ (0) = 1， 
y ′ (1) = 2，试求出满足该条件的微分方程的解析解，并绘制解的函数曲线。 
2.5试求出下面微分方程的通解。
{ 
x ′′ (t)+ y ′ (t)+ x(t) . 3y(t)=0 
（1） { 4y ′′′′ (t() t) .. 22yx ′′ (t) . x ′ (t) . 2x(t)+5y(t)=0

（2） 2x ′′ (t)+2x ′ (t) . x(t)+3y ′′ (t)+ y ′ (t)+ y(t)=0 


x ′′ (t)+4x ′ (t) . x(t)+3y ′′ (t)+2y ′ (t) . y(t)=0 
2.6试求出下面的变系数微分方程的通解。 
4 32 
xy(4)(x) + 14xy ′′′ (x) + 55xy ′′ (x) + 65xy ′ (x) + 16y(x)=0假设 y(0)=y(丘)=1，y ′ (0)=y ′ (丘)=1，试得出微分方程的解，并绘制解曲线。 
2.7试求解下面微分方程的通解以及满足条件 x(0)=1，x(丘)=2，y(0)=0的解析解。 
{ 
x ′′ (t)+5x ′ (t)+4x(t)+3y(t)= e.6t sin 4t 2y ′ (t)+ y(t)+4x ′ (t)+6x(t)= e.6t cos 4t 
2.
8试求出下面的时变线性微分方程的解析解，并求出 n =2且 x(0) = 1，x ′ (0) = 0的解析解，绘制出微分方程解的曲线。

（1）Legendre微分方程：(1 . t2)x ′′ (t) . 2tx ′ (t)+ n(n + 1)x(t)=0

（2）Bessel微分方程：t2 x ′′ (t)+ tx ′ (t)+(t2 . n 2)x(t)=0 



2.
9试求出微分方程 y ′′ (x) . (2 . 1/x) y ′ (x) + (1 . 1/x) y(x)= x 2e.5x的解析解通解，并求出满足边界条件 y(1) =丘，y(丘)=1的解析解。 

2.
10试求解下面的微分方程边值问题，并绘制解的曲线，检验是否满足边值条件。 y ′′ (x)+ xy ′ (x)+ y(x)= cos x, y(0) = 0,y(2) = 1 

2.
11试求出下面微分方程的通解。

（1）x ′′ (t)+2tx ′ (t)+ t2x(t)= t +1

（2）y ′ (x)+2xy(x)= xe.x 2

（3）y ′′′ (t)+3y ′′ (t)+3y ′ (t)+ y(t)= e.t sin t 



2.12试利用 
MATLAB的微分方程解析解函数直接求解习题 2.3，并检验得出结果的正确性。 

2.
13试求解下面一阶微分方程的通解。如果 y(0) = 1，试求解并绘制该微分方程。 dy(x)2x 


2)3e.2x
. y(x) = (4+ x
dxx2 +4 
2.
14试求解下面的非线性微分方程的解析解。

（1）y ′ (x)= y4(x) cos x + y(x) tan x

（2）xy2(x)y ′ (x)= x2 + y2(x)

（3）xy ′ (x)+2y(x)+ x5y3(x)ex =0 



2.15已知初值为 
x(0) = x ′ (0) = 1，y(0) = y ′ (0) =0，试求出下面方程的解析解，并绘制出 (x, y)的轨迹曲线。 


{ 
′′ (t) . x
(2x ′ (t)+9x(t)) . (y 	′′ (t)+ y ′ (t)+3y(t)) = 0 ′′ (t) . y
(2x ′′ (t)+ x ′ (t)+7x(t)) . (y ′ (t)+5y(t)) = 0 
2.
16试求解下面的时变线性微分方程。

（1）(x 2 . 2x + 3)y ′′′ (x) . (x 2 + 1)y ′′ (x)+2xy ′ (x) . 2y(x)=0

（2）x 2 ln xy ′′ (x) . xy ′ (x)+ y(x)=0

（3）(et + 1)y ′′ (t) . 2y ′ (t) . et y(t)=0 



2.
17试求解下面的非线性微分方程。

（1）(y ′′ (t))2 =叫(ty ′ (t) . y(t)) + Δy ′ (t)+ Ψ

（2）y ′′ (t)=(叫eβyy(t)+ Δ)y ′ (t)

（3）t5 y ′′′ (t)= a(ty ′ (t) . 2y(t))

（4）y ′ (x)= y(x)(2x2 . xy(x)+ y2(x)) 


·44·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 

2.
18试求出下面非线性微分方程的解析解。 


2x . y(x) 
)2
2 ′′ (x) .2
xy(x)y 2xy ′ (x)+ axy(x)y ′ (x)+ ay 2(x)=0 
]
e.t + cos t。 
.3 .20 .2 10 
.. 
3203 31 
x ′ (t)=x(t)+u(t)
22 .12 21 .2 .20 .3 10 
..
(
[.. 
T2.19sin试求解下面的线性微分方程。其中， (0)[0121]()()tttx= u= u=，，,,, 12..
.
. 
微分方程的初值问题  
第 3章 
MATLAB微分方程求解 


前面介绍了微分方程的解析解方法，同时指出很多非线性微分方程是不存在解析解的，需要使用数值解法对其进行研究。本章开始探讨各类微分方程的数值求解方法。 3.1节首先给出一阶显式微分方程和初值问题的一般数学模型。 3.2节介绍定步长 Runge–Kutta算法，并介绍定步长多步算法，给出这些算法的 MATLAB底层实现，通过例子比较各种数值算法的精度。在实际应用中，一般不使用定步长方法，因为该方法的精度、效率等指标都不理想，通常采用带有精度监控的自适应变步长方法。 3.3节介绍 MATLAB提供的变步长微分方程求解算法与函数调用格式，并给出一般微分方程的求解步骤，通过例子演示各种微分方程的数值求解与分析过程。 3.4节介绍微分方程数值求解的一个关键步骤解的验证，并通过例子演示微分方程数值解的验证过程，确保得到的解的正确性；此外，还介绍了更高精度的数值求解方法、变步长方法的步长显示等问题。
本章前一部分主要介绍微分方程数值求解的一些思路和早期算法，旨在帮助读者全面了解微分方程的数值解与算法实现。如果读者只关心如何利用 MATLAB求取微分方程初值问题的数值解，可从 3.3节开始阅读。 


3.1一阶显式微分方程组的初值问题
一般微分方程的数值解法很大一类是关于微分方程初值问题的数值解法。这类问题需要使用一阶显式微分方程组描述整个微分方程，然后介绍通用的求解方法，得出方程的数值解。本节首先给出一阶显式微分方程组初值问题的数学形式，然后给出解的存在性与唯一性定理。 
3.1.1初值问题的数学形式
本节先给出微分方程组的一般显式微分方程组数学描述，并为以后微分方程数值求解奠定必要的基础。
·46·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
　　定义 3-1一阶显式微分方程组的数学模型为 
x ′ (t)= ft, x(t)（3-1-1）
其中，未知变量向量 x(t)= [x1(t),x2(t)(, ··· ,x)n(t)]T称为状态向量，函数向量 
[]T
f(·)= f1(·),f2(·), ··· ,fn(·)可以是任意非线性函数。
[]T
　　定义3-2如果已知状态变量向量的初值 x0(t0)=x1(t0),x2(t0), ···,xn(t0)，且已知定义 3-1描述的一阶显式微分方程组，则该问题称为微分方程的初值问题。
初值问题的数值解法是，用数值方法推演出微分方程在某个时间区间 [t0,tn]各个时刻状态变量 x(t)的数值，这里 tn又称为终止时间。微分方程数值解需要研究的是，如何采用相应的算法与工具，可靠地得出未知状态变量在 t → [t0,tn]时的演化过程。
从给出的初值问题数学模型看，只要微分方程可以表示成一阶显式微分方程组，就可以直接采用数值方法得出问题的解。与解析解不同，方程的数值解与方程的结构本身关系不大，线性与非线性微分方程的求解方法没有显著差异。后面章节还将介绍，如果微分方程不是以标准型形式给出的，则可将其变换成相应的标准型，再进行数值求解。 

3.1.2初值问题解的存在性与唯一性
在微分方程理论中有两个重要的定理，分别描述初值问题解的存在性与唯一性，这里列出这两个定理 [17]，并作必要的解释。
()
　　定理 3-1（存在性定理）：假设 ft,x(t)在矩形区域 a<t<b，c<x(t) <d

()
内连续，且初值 t0,x(t0)也定义在这个区域内，则存在一个 .> 0，使得在区间 t0 . .<t<t0 + .上定义的函数 x(t)满足式（ 3-1-1）。

()
也就是说，如果 ft, (t)是连续的函数，则微分方程初值问题的解是存在的。但是，这并不意味着不连续的函数不存在解。

()
　　定理 3-2（唯一性定理）：假设 ft,x(t)在矩形区域 a<t<b，c<x(t) <d

()
内连续，初值 t0,x(t0)也定义在这个区域内，且 df/dx也是连续函数，则存在一个 .> 0，使得在时间区间 t0 ..<t<t0 + .内，微分方程初值问题的两个解函数 x1(t)和 x2(t)满足 x1(t)= x2(t)。也就是说，解是唯一的。 



3.2定步长数值算法与实现
传统的微分方程数值解研究起源于定步长求解方法，即：在初始时刻 t = t0，从给定的初始状态 x(t0)出发，由数值解方法求出一步近似 x(t0 + h)，再以此为起点求出 x(t0 +2h)，用这样的方法一步一步推演下去，则可找出感兴趣区间内状态向量的数值解。
本节首先介绍一些常用的定步长方法，例如最原始的 Euler算法、二阶与四阶 Runge–Kutta算法、高阶 Runge–Kutta算法等，并介绍 Adams–Bashforth算法、 Adams–Mouton算法等线性多步算法，试图通过例子演示计算误差与误差传播方面的内容，为进一步学习更好的数值解方法打下一定的基础。 
3.2.1 Euler算法
对于多元非线性常微分方程初值问题， Euler算法是最直观的一类求解算法。 Euler算法是以瑞士数学家 Leonhard Euler（1707.1783）命名的。该方法是 Euler在其 1768年的著作 “Institutionum Calculi Integralis”（积分学基础）中给出的，是最早的常微分方程数值解方法。
虽然该算法比较简单，但理解该算法对学习其他复杂的微分方程算法很有帮助，故这里以 Euler算法为例介绍微分方程初值问题的数值算法与 MATLAB实现。
假设已知 t0时刻系统状态向量的初值为 x(t0)，选择一个很小的计算步长 h，由导数定义可以将式（ 3-1-1）直接写成 
lim x(t0 + h) . x(t0)= f(t, x(t))（3-2-1）
h 0 t0 + h . t0 当然，从数值求解角度看，不能真正使 h ～ 0，而应该取一个比较小的步长 h。在这种情况下，式（ 3-2-1）左侧可以人为地去掉极限运算，函数 x(t)在 t0 + h时刻的
取值可以近似地写成 x.(t0 + h) . x(t0)+ hf(t0, x(t0))（3-2-2）
更严格地，这样的近似解存在误差，所以可以写出 t0 + h时刻系统状态向量的
值为 x(t0 + h)= x.(t0 + h)+ R0 = x(t0)+ hf(t, x(t0)) + R0（3-2-3）
利用 Taylor幂级数展开去逼近 x(t0 + h)，则 
x(t0 + h)= x(t0)+ hx ′ (t0)+ h2 2 x ′′ (t0)+ o(h3)（3-2-4）
可见，比较两个表达式得出的累积误差为 h2x ′′ (t0)/2，在数学上称该算法是 o(h)精度的。在实际求解过程中，为简单起见，通常舍弃式（ 3-2-2）中的 “.”记号，而将数值解直接记为 x1。此外，x(t0 + kh)一般简记为 xk。
　　定理 3-3一般地，假设已知 tk时刻系统的状态向量为 xk，则 tk + h时刻 Euler
算法的数值解可以写成 xk+1 = xk + hf(t, xk)（3-2-5）
·48·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
用递推的方法可以由给定的初值问题逐步求出所选择的时间段 t → [0,T ]内各个时刻 t0 + h, t0 +2h, ···处的原问题的数值解。
根据给出的 Euler算法，可以编写出如下的 MATLAB函数来解一阶显式微分方程的数值解。 
function [t,x]=ode_euler(f,vec,x0) 

if length(vec)==3, h=vec(2); t=[vec(1):h:vec(3)]'; 

else, h=vec(2)-vec(1); t=vec(:); end 

n=length(t); x0=x0(:); x=[x0'; zeros(n-1,length(x0))]; 

for i=1:n-1 

x1=f(t(i),x0); x1=x0+h*x1(:); %求解一步 x(i+1,:)=x1'; x0=x1; %保存并更新得出的解 end注意，程序中 (:)的作用是强制将向量转换成列向量，增加容错功能。即使将微分方程的函数不小心写成行向量，该程序仍然能正确执行。该函数的调用格式为 [x,t]=ode_euler(f,vec,x0)，其中，微分方程用函数句柄 f表示，后面将通过例子演示其定义方法；x0为状态变量初值向量；vec有两种定义方法，其一是给出实际的定步长向量 t，另一种方法是给出参数向量 [t0,h,tn]，其中 h为计算步长， (t0,tn)为感兴趣的时间区间。返回的 t为时间向量， x为 n列矩阵，每一列对应一个状态变量的数值解， n为状态变量的个数。下面通过例子演示 Euler算法的使用方法和效果。
　　例 3-1考虑例 2-19给出的二阶微分方程，如果 λ =2，a =1，b =0，可以重新写出该微分方程模型，为方便起见，将 yi(t)改写成 xi(t)。 
{ 
′ 
x1(t)= x2(t) 
′ 
x2(t)= .2x1(t),x1(0) = 1,x2(0) = 0
试得出该微分方程的解析解。另外，选择不同的计算步长，利用 Euler算法获得该微分方程的数值解，评价计算步长对精度的影响。　　解采用例 2-19中给出的方法，可通过下面的命令得出方程的解析解。 
>> syms t x1(t) x2(t); [x10,x20]=dsolve(diff(x1)==x2,diff(x2)==-2*x1,... x1(0)==1,x2(0)==0) %求解析解得出方程的解析解为 x1(t)= cos ∨ 2 t，x2(t)= .∨ 2 sin ∨ 2 t，后面将利用该解评价数值解的精度。
考虑采用数值方法求解上述微分方程。首先，将该方程以 MATLAB能够接受的格式输入 MATLAB环境。若给出了一阶显式微分方程组的模型，则应该编写一个 MATLAB函数或匿名函数，由当前时刻 t和 t时刻状态变量向量的值，计算出 t时刻的
向量 x ′ (t)，标准的匿名函数可以如下给出： >> f=@(t,x)[x(2); -2*x(1)]; %用匿名函数定义微分方程
这样，就可以通过下面的语句求取不同计算步长 h =0.001、0.01和 0.02时的数值解，同时计算出最大误差分别为 e1 =0.0086，e2 =0.0874，e3 =0.1789。解析解与各个数值解的曲线如图 3-1所示。显然，当步长稍大时，得出的解与理论值相差较大，可以得出结论，Euler算法的精度不是很理想。 
>> 	t0=0; tn=5; sq=sqrt(2); fplot([x10,x20],[0,5]) 
H=[0.001,0.01,0.02,0.05,0.1,0.2]; E=[]; 
for h=H %尝试不同步长时微分方程的数值解 

[t1,x1]=ode_euler(f,[t0,h,tn],[1;0]); sol=[cos(sq*t1), -sq*sin(sq*t1)]; %求解析解 e1=norm(sol-x1,inf); E=[E,e1]; line(t1,x1) 
end 


图 3-1不同计算步长的数值解比较
在图 3-1中，未标注的实线表示解析解与 h =0.001时的数值解，二者接近重合；虚线表示 h =0.01时的数值解，点线表示 h =0.02时的数值解，点画线表示 h =0.05时的数值解，可以看出，这时误差比较大。即使 h =0.001时数值解和解析解在图中也看不出区别，由于最大误差为 e =0.0086，而 h增大，则计算误差明显增大，得出的结果偏离理论值曲线。即使 h =0.001时的误差也远远大于一般意义下双精度数值解的误差，所以需要更好的求解算法。
前面介绍过，微分方程有两种描述方法，刚才演示了匿名函数的描述方法，另一种方法需要建立一个实体的 M函数文件，例如可以建立如下的文件。 
function dx=c3mtest(t,x) 

dx=[x(2); -2*x(1)]; %用 MATLAB函数描述微分方程
这样，就可以用下面的语句求解微分方程，得出的结果与前面的完全一致。 
>> 	h=0.001; [t1,x1]=ode_euler(@c3mtest,[t0,h,tn],[1;0]); 

·50·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 

3.2.2二阶 Runge–Kutta算法
由于 Euler算法计算精度比较低，所以需要精度更高的数值算法，Runge–Kutta方法是一族应用较广的数值算法。 Runge–Kutta算法是以两位德国数学家 Carl David Tolmé Runge（1856.1927）与 Martin Wilhelm Kutta（1867.1944）命名的算法。其中， Runge在 1895年提出了在一步之内求解多步导数的思想， Kutta给出了四阶及以下阶次的一系列计算公式。 
Euler算法也属于特殊的 Runge–Kutta算法，是一阶算法。这里介绍另一种简单的 Runge–Kutta算法二阶 Runge–Kutta算法。
　　定理 3-4已知 tk时刻系统的状态向量为 xk，则 tk + h时刻二阶 Runge–Kutta
算法的数值解可以写成 xk+1 = xk + h 2(k1 + k2)（3-2-6）
其中，中间变量为 {
k1 = f(tk, xk) （3-2-7）k2 = f(tk+1, xk + hk1)
可以证明，该算法的累积误差为 o(h2)。 
o(h)与 o(h2)的误差在实际应用中有什么区别呢？可以由下面的例子描述累积误差：如果计算步长选择为 h =0.01，则 o(h)的累积误差是 0.01级别的，而 o(h2)累积误差是 0.012，即 0.0001级别的。显然，对较小的计算步长 h而言， o(h2)算法的累积误差明显小于 o(h)算法。由于计算步长 h的取值一般很小，所以在算法选择时，应该优先选择 p值较大的 o(hp)算法。
采用与前面 Euler算法完全一致的函数框架，可以编写出二阶定步长 Runge– Kutta算法的数值求解函数，该函数的调用格式与 ode_euler()函数完全一致。 function [t,x]=ode_rk2(f,vec,x0) if length(vec)==3, h=vec(2); t=[vec(1):h:vec(3)]'; else, h=vec(2)-vec(1); t=vec(:); end n=length(t); x0=x0(:); x=[x0'; zeros(n-1,length(x0))]; for i=1:n-1 k1=f(t(i),x0); k2=f(t(i+1),x0+h*k1); x1=x0+h*(k1(:)+k2(:))/2; x(i+1,:)=x1'; x0=x1; end
　　例 3-2试用二阶 Runge–Kutta法重新求解例 3-1中的问题，并评价精度。　　解可以用下面的语句直接求解该微分方程问题，尝试不同的计算步长，得出的误差分别为 e1 =4.0246×10.6，e2 =4.0299×10.4，e3 =0.0016，明显低于 Euler算法的误差。也可以选择更大的计算步长 h =0.1，这时误差为 e4 =0.0407。得出的数值解曲线
比较如图 3-2所示。可见，除了步长 h =0.1时的结果可以分辨出来，其他结果非常相近。即使选择大步长 h =0.2，得出的结果也明显优于前面介绍的 Euler算法。 
>> 	f=@(t,x)[x(2); -2*x(1)]; t0=0; tn=5; sq=sqrt(2); H=[0.001,0.01,0.02,0.05,0.1,0.2]; E=[]; for h=H %尝试不同的步长求解微分方程 
[t1,x1]=ode_rk2(f,[t0,h,tn],[1;0]); sol=[cos(sq*t1), -sq*sin(sq*t1)]; %求解析解 e1=norm(sol-x1,inf); E=[E,e1]; line(t1,x1) 
end


)
(
图 3-2不同计算步长的数值解比较
可以看出，与前面介绍的 Euler算法相比，同等计算步长时的误差显著减小了，但对稍大的计算步长而言，计算误差仍然很大。 

3.2.3四阶 Runge–Kutta算法
四阶定步长的 Runge–Kutta算法是传统数值分析课程和系统仿真课程中经常介绍的算法，被认为是求解微分方程的一种有效的方法，该算法结构很简单，易于计算机实现。
　　定理 3-5四阶 Runge–Kutta算法可以递推求解出状态变量值为 
xk+1 = xk + h k1 +2k2 +2k3 + k4（3-2-8）
6
其中，四个中间向量可以如下计算： 
. 

k1 
= f(tk, xk) hh
)
)
（3-2-9）
.
..

k2 = f
tk + , xk + k1
22
(
(
hh 

k3 = f
tk + , xk + k2
22 
.
.. 

k4 = f(tk + h, xk + hk3)
·52·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
其中，h为计算步长，在实际应用中该步长是一个常数。
可以证明，四阶 Runge–Kutta算法的累积误差为 o(h4)。仍以 h =0.01为例，四阶 Runge–Kutta算法的累积误差可能达到 0.014，即 10.8级别，故该算法的精度远远高于 Euler算法与二阶 Runge–Kutta算法。
因此，采用递推的方法可由给定的初值问题逐步求出所选择的时间段 t→[t0,tn]内各个时刻 t0 + h, t0 +2h, ···的原问题的数值解。有了上面的数学算法，可以编写出四阶定步长 Runge–Kutta算法的 MATLAB语言函数如下所示。 
function [t,x]=ode_rk4(f,vec,x0) 
if length(vec)==3, h=vec(2); t=[vec(1):h:vec(3)]'; 
else, h=vec(2)-vec(1); t=vec(:); end 
n=length(t); x0=x0(:); x=[x0'; zeros(n-1,length(x0))]; 
for i=1:n-1 %用循环结构求取微分方程的数值解 

k1=f(t(i),x0); k2=f(t(i)+h/2,x0+h*k1/2); 
k3=f(t(i)+h/2,x0+h*k2/2); k4=f(t(i+1),x0+h*k3); 
x1=k1+2*k2+2*k3+k4; x1=x0+x1(:)*h/6; 
x(i+1,:)=x1'; x0=x1; %更新初值，并保存本步得出的数值解 

end

　　例 3-3试用四阶 Runge–Kutta算法重新求解例 3-1中的问题，并评价精度。　　解使用四阶 Runge–Kutta算法，选取计算步长 h =0.001，h =0.01，h =0.1和更大的计算步长 h =0.2，得出的最大误差分别为 e1 =4.0712×10.13，e2 =4.0307×10.9， e3 =4.0756×10.5，e4 =0.00065。可以看出，选择四阶 Runge–Kutta算法，计算精度较其他方法有了显著的提高。即使计算步长取到 h =0.2，该算法仍然可用，得出的误差小于 Euler算法取 h =0.001时的误差。因此，该算法的效率明显优于 Euler算法与二阶 Runge–Kutta算法。 
>> 	f=@(t,x)[x(2); -2*x(1)]; t0=0; tn=5; sq=sqrt(2); H=[0.001,0.01,0.02,0.05,0.1,0.2]; E=[]; for h=H %尝试不同的步长求解微分方程 
[t1,x1]=ode_rk4(f,[t0,h,tn],[1;0]); sol=[cos(sq*t1), -sq*sin(sq*t1)]; %求解析解 e1=norm(sol-x1,inf); E=[E,e1]; %计算微分方程解的最大误差 
end 


3.2.4 Gill算法 
Gill算法是四阶 Runge–Kutta算法的改进形式。该算法对其中的参考点位置与加权系数改进后的结果，该算法如下给出。　　定理 3-6 Gill算法可以递推求解出状态变量值为 
h[(∨ )(∨ )]
xk+1 = xk + k1 +2 . 2k2 +2+ 2k3 + k4（3-2-10）
6 
其中，四个中间向量可以如下计算： 
. 
k1 = ftk, xk
. 

hh 
k2 = f tk + , xk + k1
.

. 22
(	)
. 
k3 = f ((tk + ) h, xk + h ( ∨) 2 . 1)k1 + h (2 .∨ 2)k2（3-2-11）
22 2

(	)
∨∨
hh 
.

. 	k4 = f tk + h, xk . 2k2 + (2+ 2)k3
22 
其中，h为计算步长。Gill算法的累积误差也是 o(h4)。
根据上述算法，可以写出下面的 MATLAB实现。 
function [t,x]=ode_gill(f,vec,x0) 
if length(vec)==3, h=vec(2); t=[vec(1):h:vec(3)]'; 
else, h=vec(2)-vec(1); t=vec(:); end 
n=length(t); x0=x0(:); x=[x0'; zeros(n-1,length(x0))]; 
for i=1:n-1 

k1=f(t(i),x0); k2=f(t(i)+h/2,x0+k1*h/2); 
k3=f(t(i)+h/2,x0+(sqrt(2)-1)*h*k1/2+(2-sqrt(2))*h*k2/2); 
k4=f(t(i+1),x0-sqrt(2)*h*k2/2+(2+sqrt(2))*h*k3/2); 
x1=k1+(2-sqrt(2))*k2+(2+sqrt(2))*k3+k4; 
x1=x0+x1(:)*h/6; x(i+1,:)=x1'; x0=x1; 

end

　　例 3-4试用 Gill算法程序重新求解例 3-1中的问题，并评估算法的精度。
　　解仿照前面的求解语句，可以直接给出下面的求解命令，得出各种计算步长时的
求解结果，并得出残差向量的范数。 
>> 	f=@(t,x)[x(2); -2*x(1)]; t0=0; tn=5; sq=sqrt(2); 
H=[0.001,0.01,0.02,0.05,0.1,0.2]; E=[]; 
for h=H 

[t1,x1]=ode_gill(f,[t0,h,tn],[1;0]); 
sol=[cos(sq*t1), -sq*sin(sq*t1)]; 
e1=norm(sol-x1,inf); E=[E,e1]; 

end 



3.2.5 m阶Runge–Kutta算法
前面介绍了若干数值求解算法，其中，Euler算法可以认为是 Runge–Kutta
算法的一个特例，该算法是一阶 Runge–Kutta算法，另外还给出了二阶与四阶
·54·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 
Runge–Kutta算法。从算法的结构看，它们之间是有相似之处的，其精度也分别为 o(h2)与 o(h4)。更一般地，还可以给出一族 m阶 Runge–Kutta算法。 Runge–Kutta算法采用插值的方法，由已知初值通过定积分计算下一个点的
函数值。 x(b)= x(a)+ ∫b f(t, x(t))dt（3-2-12）
a 其中， a = tk，b = tk+1，计算步长 h = tk+1 . tk。可以在区间 (tk,tk+1)设置一些插
值点： tij = tk +叫jh, j =1, 2, ··· ,m（3-2-13）
由选定的插值点可以得出如下的 m阶 Runge–Kutta算法。
　　定理 3-7 m阶 Runge–Kutta算法的通用数学公式为 
m
∑ 
xk+1 = xk + hΨiki（3-2-14）i=1 
其中，中间变量可以如下计算： 
  
i.1
∑ 
ki = f.tk +叫ih, xk + hΔijkj. ,i =1, 2, ··· ,m（3-2-15）j=1 
可以看出，前面介绍的二阶 Runge–Kutta算法、四阶 Runge–Kutta算法及 Gill算法都是 m阶 Runge–Kutta算法的特例，且 m阶 Runge–Kutta算法的累积误差为 o(hm)。
　　定理3-8 Runge–Kutta算法中，各系数之间的关系可以由 Butcher表（Butcher 
tableau）描述 [6]，如下所示：
叫1 =0叫2 叫3 
.
.
. 
叫m 
其中，系数满足如下关系： 
Δ21 
Δ31 ...  Δ32 ...  ... （3-2-16） 
Δm1  Δm2  · · ·  Δm,m.1  
Ψ1  Ψ2  · · ·  Ψm.1  Ψm 

i.1∑ Δij =叫j,  j = 2, 3, · · · , m （3-2-17） 
j=1  
且  m.1∑ i=1 Ψi叫γ.1 i = 1 λ,  λ = 1, 2, · · · , m . 1 （3-2-18）  
∑ ΨiΔij叫j = 1 6 ,  ∑ Ψi叫iΔij叫j = 1 8 , · · · （3-2-19） 

Δ21 =叫2,Δ31 + Δ32 =叫3 
Ψ3Δ32叫2 = 16
a2=1/2; a3=1; 

Butcher发明了形象的树状记号来简记系数项的组合形式。例如，使用实心圆	3-2-18点和记号、分别表示式（ ）中和的项。|∞λ =1λ =2λ =3.，..3-5Runge–Kutta　　例 试列写出三阶算法的待定系数方程，并写出该算法的公式。　　解选择则可以建立起如下六个方程。=3m ，.2.
.1122ΨΨΨ=1Ψ叫Ψ叫Ψ叫Ψ叫++++== 123 , 2233 , 232323.
. 从上述方程可知，有八个未知变量需要求解，所以需要人为地指定其中两个变量。MATLAB为简单起见，一般给出的值，例如选择。利用符号运算可叫叫 =1/2叫 =1，i23 以直接求解相应的代数方程组。 symsb21b31b32g1g2g3>> [b21b31b32g1g2g3]=solve(b21==a2,b31+b32==a3,... g1+g2+g3==1,g2*a2+g3*a3==1/2,g2*a2^2+g3*a3^2==1/3,... RungeKuttag3*b32*a2==1/6)%求算法的系数.Runge–Kutta由此可以写出三阶算法的公式为 
得出方程的解为 
1 	121 
Δ21 = ,Δ31 = .1,Δ32 =2,Ψ1 = ,Ψ2 = ,Ψ3 = 
636
h 
xk+1 = xk +(k1 +4k2 + k3)
6 
. 
(
(()
h 
其中， 
k1 = f
k2 = f
tk, xktk + 
.
. .
. 
k1
)
h 
, xk + 
2 2 
)
重新求解相应的
k3 = f
tk + h, xk + h(.k1 +2k2)
当然，叫2和叫3
也可以选择为其他数值。例如，叫2 = 1/2，叫3 = 3/4，
代数方程组。 
>> 	syms b21 b31 b32 g1 g2 g3 
a2=1/2; a3=3/4; 
[b21 b31 b32 g1 g2 g3]=solve(b21==a2,b31+b32==a3,... 

g1+g2+g3==1,g2*a2+g3*a3==1/2, g2*a2^2+g3*a3^2==1/3,... 
g3*b32*a2==1/6)

得出方程的解为 
1 	3214 
Δ21 = ,Δ31 =0,Δ32 = ,Ψ1 = ,Ψ2 = ,Ψ3 = 
2 	4939
由此可以写出另一个三阶 Runge–Kutta算法的公式（又称为 Ralston公式）为 
h 
xk+1 = xk + (2k1 +3k2 +4k3)
9 
·56·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
其中， 	. k1 = f((tk, xk)).

. h	h 
k2 = f tk + , xk + k1
22
(	)
2h 3h 
.

. k3 = f tk + , xk + k2
34 
这两个不同算法的 Butcher表分别描述如下： 
0  0  
1/2 1  1/2 .1  2  1/2 2/3  1/2 0  3/4  
1/6  2/3  1/6  2/9  1/3  4/9 

四阶或更高阶 Runge–Kutta算法公式的推导比较烦琐，可参考文献 [18, 19]。其中，四阶 Runge–Kutta算法公式的推导可以建立起 11个代数方程，含有 13个未知数。例如，需要人为指定叫2和叫3的值，再去求解其他待定系数。 

3.2.6定步长多步算法与实现
前面介绍的 Runge–Kutta算法属于单步方法，即根据 xk求解 xk+1的方法。在求解过程中需要得出 ki，这些中间值可以认为是一个步长 h内的插值。在实际应用中，如果已知一些初始点 x0，x1，···的函数值，则可以考虑引入多步算法。
线性多步法的主要思路是，将微分方程的解写成线性差分方程的形式，从而使用递推方法得出微分方程的近似解。线性多步法的一般公式为 
xk+m = .amxk+m.1 .· ··. a1xk 
[	]（3-2-20）
+hbmf(tn+m, xk+m)+ ··· + b1f(tk, xk)
线性多步法需要确定常数 ai、bi，i =1, 2, ··· ,m，用递推方法得出微分方程的数值解。比较常用的线性多步方法是四阶 Adams–Bashforth算法与 Adams– Mouton算法。
　　定理 3-9取 am =1，ai =0，i =1, 2, ··· ,m . 1，且计算出合适的 bi，则可以由下面的四阶 Adams–Bashforth公式计算微分方程的数值解。 
h 
xk+1 = xk + 55f(tk, xk) . 59f(tk.1, xk.1)24 []（3-2-21）+37f(tk.2, xk.2) . 9f(tk.3, xk.3)
其中，k =3, 4, ···，该算法的精度为 o(h4) [18]。
　　定理 3-10基于四阶 Adams–Mouton公式的微分方程的数值解为 
h 
xk+1 = xk +9f(tk+1, xk+1) + 19f(tk, xk)24 []（3-2-22）+37f(tk.1, xk.1) . 9f(tk.2, xk.2)
其中，k =2, 3, ···，该算法的精度为 o(h4) [18]。
Adams–Mouton算法中，计算 xk+1时等号右边也含有 xk+1，这就给计算带来了困难，通常需要引入预估 –校正算法来实现，此处不介绍该算法的 MATLAB实现，只探讨 Adams–Bashforth算法的 MATLAB实现。
由式（ 3-2-21）可知，需要知道 x3（或 Adams–Mouton公式中的 x2）的值，而一般情况下的微分方程已知初值 x0，其他初值例如 x3不是已知的，需要采用其他方法得出。传统的线性多步方法的做法是，采用 Euler算法计算 x1，由 x0与得出的 x1，采用低精度的梯形近似法求 x2，再用三步法求 x3。显然，这样的方法在初始点生成时可能伴随有巨大的误差，严重影响后续计算的精度，往往不可取。所以，应该采用高精度算法重新构造这些初始点。例如，可以将四阶 Runge–Kutta算法求解程序嵌入多步法程序求解 x1、x2和 x3。 
function [t,x]=ode_adams(f,vec,x0) 

if length(vec)==3, h=vec(2); t=[vec(1):h:vec(3)]'; 

else, h=vec(2)-vec(1); t=vec(:); end 

[t0,x]=ode_rk4(f,[vec(1),h,vec(1)+3*h],x0) 

n=length(t); x=[x; zeros(n-4,length(x0))]; 

for i=4:n-1 %利用多步方法求解微分方程的数值解 
k1=f(t(i),x(i,:)); k2=f(t(i-1),x(i-1,:)); 

k3=f(t(i-2),x(i-2,:)); k4=f(t(i-3),x(i-3,:)); 

x1=55*k1-59*k2+37*k3-9*k4; x1=x(i,:).'+h*x1(:)/24; 

x1=x1(:); x(i+1,:)=x1'; 
end

　　例 3-6试用 Adams–Bashforth公式重新求解例 3-1中的问题，并评价精度。　　解仿照前例，可以直接使用下面的语句，在不同步长下求出微分方程的数值解，并得出精度。 
>> f=@(t,x)[x(2); -2*x(1)]; t0=0; tn=5; sq=sqrt(2); 

H=[0.001,0.01,0.02,0.05,0.1,0.2]; E=[]; 

for h=H %尝试不同的步长求解微分方程 
[t1,x1]=ode_adams(f,[t0,h,tn],[1;0]); %求解微分方程 
sol=[cos(sq*t1), -sq*sin(sq*t1)];  %求解析解  
e1=norm(sol-x1,inf); E=[E,e1];  %求最大误差  
end 

对相同的微分方程，采用了不同的求解算法，不同算法的最大误差如表 3-1所示。可以看出，四阶 Runge–Kutta算法的精度与 Gill算法几乎完全一致。此外，四阶 Runge–Kutta算法的精度（ o(h4)精度）明显高于二阶 Runge–Kutta算法（ o(h2)精度），更高于 Euler算法（o(h)精度）。即使使用 h =0.2这样的大步长，由 o(h4)算法得到的结果也是具有比较高的精度的。虽然 Adams–Bashford线性多步算法也被称
·58·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
为 o(h4)精度，但相比之下，对上述具体例子而言，求解精度比 Runge–Kutta算法差很多。
表 3-1不同底层算法的最大误差比较
算法 h =0.001 h =0.01 h =0.02 h =0.05 h =0.1 h =0.2 
Euler算法 0.00856 0.08741 0.1789 0.4792 1.069 2.554二阶 RK算法 4.025×10.6 0.000403 0.00161 0.01013 0.040735 0.16457四阶 RK算法 4.071×10.13 4.03×10.9 6.46×10.8 2.535×10.6 4.076×10.5 0.00065 
Gill算法 4.070×10.13 4.03×10.9 6.46×10.8 2.535×10.6 4.076×10.5 0.00065线性多步法 1.68×10.11 1.68×10.7 2.68×10.6 0.0001 0.0016 0.0229
与 Runge–Kutta算法这类单步方法相比，多步算法的劣势在于需要用其他方法生成足够的初始点，然后才能由这些初始点搜索后续的微分方程的解。这些初始点将影响多步法的最终求解精度，而它们需要由 Runge–Kutta这类算法计算。从另一个角度看，这里介绍的多步算法只能用于定步长算法，否则不能直接使用以往的信息，多步法也就无法实现微分方程的求解了。变步长多步法的结构与实现更复杂 [5]，这里不再讨论。 


3.3变步长数值算法与实现
前面介绍的数值算法都属于定步长算法，即选择一个计算步长 h，通过递推方法从给定的初值一步一步递推出微分方程的数值解。如果选定了算法的精度阶次，则提高数值解精度的一种方法是减小步长 h的值。然而，在实际应用中并不能无限制地减小 h的值，这主要有两个原因：
（1）减慢计算速度。对于确定的求解时间，减小步长意味着增加在这个时间段内的计算点数目，故计算速度减慢。

（2）增加累积误差。不论选择多小的步长，所得到的数值解都将有一个舍入误差，减小计算步长将增加计算的次数，从而使得整个计算过程的舍入误差的叠加和传递次数增多，会产生较大的累积误差。舍入误差、累积误差和总误差关系的示意图如图 3-3所示。


本节主要探讨变步长方法，介绍 MATLAB中微分方程的变步长求解函数，并利用该函数直接求解微分方程，探讨该函数的各种调用格式与使用方法，同时给出更高效的求解函数。 
3.3.1提高求解效率的措施
虽然数值分析类课程主要介绍定步长算法，但是定步长算法存在较多局限，在实际应用中几乎不使用定步长算法来求解微分方程。要想提高求解问题的效率，则
计算
误差


舍入误差总误差
累积误差步长 h
图 3-3计算误差示意图
在对微分方程的求解过程中，应采取下列措施：
（1）选择适当的步长。采用像 Euler算法这样简单的算法时，应选择适当的步长，既不能太大，又不能太小。

（2）改进近似算法精度。由于 Euler算法只是将原积分问题进行梯形近似，其近似精度很低，因而不能有效地逼近原始问题。可以使用更精确的插值方法取代 Euler算法，从而改进运算精度。比较成功的是 Runge–Kutta算法、Adams算法等。

（3）采用变步长方法。前面提及“适当”地选择步长，这本身就是个模糊的概念，如何适当地选择步长取决于经验。事实上，很多种方法都允许变步长求解。当误差较小时，可自动增大步长，而误差较大时再自动减小步长，从而精确、有效地求解给出的常微分方程初值问题。 



3.3.2变步长方法简介
一般变步长算法的原理如图 3-4所示。已知 tk时刻的状态变量 xk，首先在某步长 h下计算出 tk + h时刻的状态变量 x.k+1。然后，将步长变成原来步长的一半，分两步从 xk计算出 tk + h时刻的状态变量 x.k+1。如果两种运算步长下的误差 . = ||x.k+1 . x.k+1||小于给定的误差限，则可以采用该步长或适当增大步长；如果误差 
一步求解结果 x.k+1
xk+1 

图 3-4变步长算法示意图
·60·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
大于误差限，则进一步减小步长。自适应变步长算法可以较好地解决计算速度问题，并且能够保证计算的精度。
另一种常用的变步长方法是，在一步计算中用两种方法计算出 xk+1和 x.k+1的值，如果 ||xk+1 . x.k+1|| . φ，则可以接受当前的步长，或增大步长以保证求解速度；如果该条件不满足，则减小步长。
定步长算法相当于用开环控制的思想求解微分方程，选定步长后将不考虑计算是不是有误差，也不考虑误差是否能被接受，一直采用单一步长计算全程。而变步长算法则采用由误差作为监测指标的闭环控制思想，在计算过程中不断修正步长，使得预期的计算精度得以保证。同时，由于采用变步长思想，大部分问题的求解时间将明显缩短，因为计算过程中在保证计算精度的前提下也允许加大步长。 
3.3.3四级五阶 Runge–Kutta变步长算法
美国国家宇航局（NASA）的研究员 Erwin Felhberg对传统的 Runge–Kutta方法进行了改进 [19]，在每一个计算步长内对函数 fi(·)进行六次求值，以保证更高的精度和数值稳定性，该算法效率很高，精度也比较理想，所以实际应用中经常使用。
　　定理 3-11假设当前的步长为 hk，则可以定义下面的六个 ki变量： 
i.1
∑ 
ki = f .tk +叫ihk, xk + Δijkj. ,i =1, 2, ··· , 6（3-3-1）j=1 
式中，tk为当前计算时刻，中间参数叫i、Δij及其他参数由表 3-2给出。其中，叫i和 Δij参数对又称为 Dormand–Prince对。这时，下一步的状态向量可以由下式求出。 
6
∑ 
xk+1 = xk + Ψihkki（3-3-2）i=1 
当然，直接采用这一方法仍为定步长的方法。在实际问题中，往往希望在一些情况（例如解变化很快时）下采用较小的步长，而在另一些情况（例如解的变化很缓
表 3-2四级五阶 RKF算法系数的 Butcher表 
0  
1/4 3/8 12/13 1 1/2  1/4 3/32 1932/2197 439/216 .8/27  9/32 .7200/2197 .8 2  7296/2197 3680/513 .3544/2565  .845/4104 1859/4104  .11/40  
ìi ì i  16/135 25/216  0 0  6656/12825 1408/2565  28561/56430 2197/4104  .9/50 .1/5  2/55 0  

慢时）下采用较大的步长，这样做既可以保证较高的精度，又可以保证较高的计算速度。在此算法中，定义误差向量 
6
∑ 
= 讨 （3-3-3）
.k (Ψi . Ψi )hki
i=1 
可以根据其大小改变计算步长，这种能自动变换步长的方法又称为自适应变步长算法。在该算法中，计算用 Ψi误差级别为 o(h5)，而检验部分的 Ψ讨i误差级别为 o(h4)，所以该算法又称为四级五阶 Runge–Kutta–Felhberg（RKF）算法。也就是说，这里的“级”就是检验的阶次，而“阶”是指实际计算的阶次。四级五阶 RKF算法系数的 Butcher表如表 3-2所示。 

3.3.4基于 MATLAB的微分方程求解函数 
MATLAB中求解一阶微分方程组初值问题数值解最常用的是 ode45()函数，该函数实现了变步长四级五阶 Runge–Kutta–Felhberg算法。该函数的调用格式为 
[t,x]=ode45(Fun,[t0,tn],x0), %直接求解 
[t,x]=ode45(Fun,[t0,tn],x0,options), %带有控制选项 
[t,x]=ode45(Fun,[t0,tn],x0,options,p1,p2, ··· ,pm), %带有附加参数其中，微分方程应该用 MATLAB函数 Fun或匿名函数按指定的格式描述，有关该函数的编写方法后面将通过例子专门介绍； tspan是描述求解区间的变元，通常情况下由 tspan=[t0,tn]给出，用于描述微分方程求解的区间，如果只给出 tn，则表示初始时刻为 t0 =0。为求解微分方程，还应该已知初值问题的初始状态变量 x0。
另外，该函数允许 tn <t0，即可以认为 t0为终止时刻， tn为初始时刻， x0表示状态变量的终值，故该函数也可以直接求解终值问题。
除了函数 ode45()外，类似的函数还有 ode113()，ode15s()，ode23()等，它们有相同的调用格式，但算法是不同的。
求解一阶显式微分方程组的关键是用 MATLAB语言编写一个函数描述需要求解的微分方程组。该函数的入口应该为 
function xd=funname(t,x), %无附加参数的格式 
function xd=funname(t,x,p1,p2,··· ,pm), %可以使用附加参数其中， t是时间变量或自变量，即使需要求解的微分方程不是时变的，也需要给出 t占位，否则 MATLAB在变量传递过程中将出现问题。变量 x为状态向量，返回的 xd为状态向量的导数。该函数允许求解带有附加参数 p1,p2, ··· ,pm的微分方程，附加参数必须与微分方程求解程序中的完全对应。下面将通过例子详细介绍相关的调用格式。
·62·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
　　例 3-7假设著名的 . Lorenz模型的状态方程表示为 
′
.x1(t)= .Δx1(t)+ x2(t)x3(t)
. 
x (t)= .λx2(t)+ λx3(t)
2′
.
. 
x (t)= .x1(t)x2(t)+ νx2(t) . x3(t)
3′
其中， Δ = 8/3，λ = 10，ν = 28。若其初值为 x1(0) = x2(0) = 0，x3(0) = .，.为机器上可以识别的小常数，例如取一个很小的正数 . = 10.10，试求解该微分方程组。　　解该方程是非线性微分方程，所以不存在解析解，只能用数值解法求解。若想用
()
数值算法求解该微分方程，需要将该微分方程变换为标准型 x ′ (t)= ft, x(t)的形式。构造状态变量向量 x(t)=[x1(t),x2(t),x3(t)]T，则原始的微分方程可以写成 
. 
.Δx1(t)+ x2(t)x3(t)
()()
x ′ (t)= ft, x(t),其中 ft, x(t)= . .λx2(t)+ λx3(t) ..x1(t)x2(t)+ νx2(t) . x3(t) 在实际求解过程中，只需将 f(t, x(t))用向量形式描述出来。例如，可以编写一个匿名函数来描述系统的动态模型。 
>> f=@(t,x)[-8/3*x(1)+x(2)*x(3); 
-10*x(2)+10*x(3); 
-x(1)*x(2)+28*x(2)-x(3)];

这时，可以调用微分方程数值解 ode45()函数对匿名函数 f描述的系统进行数值求解，并将结果进行图形显示。 
>> tn=100; x0=[0;0;1e-10]; 

[t,x]=ode45(f,[0,tn],x0); plot(t,x) %求方程的数值解其中， tn为设定的仿真终止时间， x0为初始状态。绘图命令绘制出系统的各个状态和时间的关系的二维曲线图，如图 3-5所示。

图 3-5 Lorenz方程的仿真结果图示
如果想绘制出三个状态变量的相空间曲线，则可以使用下面的语句，得出的曲线如图 3-6所示。 
>> plot3(x(:,1),x(:,2),x(:,3)); grid %相空间轨迹三维图

图 3-6 Lorenz方程的相空间轨迹三维图
看似很复杂的三元一阶非线性常微分方程组的数值解问题通过几条简单直观的 MATLAB语句就求解出来了。此外，还可以使用 MATLAB语言直观地将结果用可视化的方式显示出来，这也是本书选择 MATLAB语言作为主要语言的原因之一。
观察相空间轨迹的最好方法是采用 comet3()函数绘制动画式的轨迹，故可将上述最后一个语句改成 comet3(x(:,1),x(:,2),x(:,3))。 
Lorenz方程是美国气象学家 Edward Norton Lorenz（1919.2008）于 1963年提出的，是用于描述大气对流（atmospheric convection）的数学模型 [20]。其原始模型与这里给出的状态变量顺序稍有不同，此处采用的状态变量次序是为了更好地展示相空间曲线的“蝴蝶”形状。Lorenz方程是混沌研究的开端，有关混沌的内容详见第 7章。
下面简单总结微分方程的数值求解步骤。
（1）写出微分方程的数学形式标准型。x ′ (t)= f(t, x(t))，且已知 x(t0)。
（2）用MATLAB描述微分方程。可以用匿名函数或 MATLAB函数描述微分方程的标准型。注意，对于函数输入变元 t，即使原始微分方程不显含 t，也应该保留该变元占位，不能舍弃。
（3）微分方程的求解。可调用 ode45()或其他函数直接求解。

（4）解的检验。解的检验是微分方程求解的一个重要环节，后面将专门介绍。


有了规范的微分方程数值解方法与步骤，就可以通过例子进一步学习与掌握一般微分方程的直接求解方法与技巧。
　　例 3-8例 1-5曾给出 Lotka–Volterra捕食者与猎物模型方程，现在重新给出如下。试求解该微分方程，并绘制响应曲线。 
{ 
x ′ (t)=4x(t) .叫x (t)y(t) y ′ (t)= Δx(t)y(t) . 3y(t)
其中，叫 =2，Δ =1，且初值为 x(0) = 2，y(0) = 3。
·64·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
　　解本例要求解的函数为 x(t)和 y(t)，而标准型中要求的向量 x(t)，所以需要引入
状态变量。例如，令 x1(t)= x(t)，x2(t)= y(t)，就可以写出如下的标准型模型。 
[ 	][]
x ′ (t)= f(t, x(t)),其中 f(t, x(t))= 4x1(t) . 2x1(t)x2(t) , x(0)= 2 
x1(t)x2(t) . 3x2(t)3 
有了标准型模型，可用匿名函数直接描述该微分方程，并进行求解，最后绘制出 
x(t)与 y(t)之间的关系曲线（又称为相平面曲线），如图 3-7所示。另外，可以测出该求解
语句耗时 0.00648s，计算的点数为 153。 
>> 	x0=[2; 3]; tn=10; 
f=@(t,x)[4*x(1)-2*x(1)*x(2); x(1)*x(2)-3*x(2)]; 
tic, [t,x]=ode45(f,[0,tn],x0); toc %求数值解并计时 
length(t), plot(x(:,1),x(:,2)) %绘制解的相平面曲线


图 3-7微分方程的相平面曲线（比较粗糙）
从上述微分方程求解实例看，如果能建立一个描述微分方程组的 M函数或匿名函数，则调用 ode45()可以立即解出方程的数值解。因此，编写 MATLAB函数描述微分方程是常微分方程初值问题数值求解的关键。
下面再看一个更复杂的例子及其求解方法。
　　例 3-9试求解下面给出的三天体运动模型 [6]，这里略去了 (t)。 
. 
′ 
y = y4
1 
′ 
y = y5
2 
′ 
y = y6
3 
.
μ(y1 + μ . 1) (1 . μ)(y1 + μ)
. ′ 
y4 =2y5 + y1 . √ 	. √ 
22 	22
(y2 + y+(y1 +μ.1)2)3 (y2 + y+(y1 +μ)2)3 
3	3 
μy2 	(1 . μ)y2
′ 
y = .2y4 + y2 . √ 	. √
5 
22 	22
(y2 + y+(y1 +μ.1)2)3 (y2 + y+(y1 +μ)2)3 
3	3 
μy3 	(1 . μ)y3
′ 
.
y = .√ . √ 
22 	22
. 	6 (y+(y1 + μ . 1)2)3 (y+(y1 + μ)2)3
2 + y3	2 + y3 
其中，两个较大的天体是地球、月球，最小的天体是卫星。y1、y2和 y3为卫星的空间
坐标， y4、y5和 y6为卫星在三个坐标系下的速度分量。假设卫星状态变量的初始值为 y0 = [0.994, 0, 0, 0, .2.0015851063790825224, 0]T，常数 μ = 1/81.45。试绘制卫星运动的空间轨迹。　　解观察给出的一阶显式微分方程组，可以用匿名函数直接描述，并调用 ode45()函
√ 
数进行求解。从数学模型表示上，可以看到很多项都有公共项 (y22 + y32 +(y1 + μ)2)3
√ 
与 (y22 + y32 +(y1 + μ . 1)2)3，如果逐项描述是比较烦琐的，所以最好将其作为两个中间变量预先计算出来。由于匿名函数不支持中间变量的使用，所以这里只能用 MATLAB函数来描述微分方程。 
function dy=threebody(t,y) 
mu0=1/81.45; 
D1=sqrt((y(2)^2+y(3)^2+(y(1)+mu0-1)^2)^3); 
D2=sqrt((y(2)^2+y(3)^2+(y(1)+mu0)^2)^3); 
dy=[y(4:6); 

2*y(5)+y(1)-mu0*(y(1)+mu0-1)/D1-(1-mu0)*(y(1)+mu0)/D2; 
-2*y(4)+y(2)-mu0*y(2)/D1-(1-mu0)*y(2)/D2; 
-mu0*y(3)/D1-(1-mu0)*y(3)/D2];

令 t0 =0，tn = 40，则可以调用下面的语句求解该微分方程，并绘制出卫星的轨迹曲线，如图 3-8所示。然而，该结果是错误的，这也带来了一个新的问题：如何判定结果的正确性。 
>> 	y0=[0.994,0,0,0,-2.0015851063790825224,0]'; 
[t,y]=ode45(@threebody,[0,40],y0); %求微分方程数值解 
plot3(y(:,1),y(:,2),y(:,3)) %绘制卫星的相空间轨迹


图 3-8卫星的三维空间轨迹（结果是错误的） 


3.3.5基于 MATLAB的带有附加参数的微分方程求解
在基于 MATLAB的微分方程求解中，引入附加参数的目的是，当微分方程的某些参数可以选择不同的值时，可以避免对不同值进行求解时对模型文件的修改。
·66·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解
例如，例 3-7中给出的 Lorenz方程中， Δ、λ和 ν可以看成附加参数，当表示它们的参数变化时，无须修改描述原始微分方程的 MATLAB函数，而只需在调用求解器时从外部进行修改。
　　例 3-10试编写带有附加参数的 MATLAB函数描述例 3-7中的 Lorenz方程，并利用该函数分别求解在例 3-7给定参数下和一组新参数 Δ =2，λ =5，ν = 20时 Lorenz方程的数值解。　　解选定附加参数为 Δ、λ和 ν，编写如下的 MATLAB函数来描述给定的常微分方程组。 
function dx=lorenz1(t,x,b,r,s) %带有附加参数的微分方程描述函数 dx=[-b*x(1)+x(2)*x(3); 
-r*x(2)+r*x(3); -x(1)*x(2)+s*x(2)-x(3)]; %描述微分方程然后求取方程的数值解。调用函数时无须使用和函数本身完全一致的变量名，只要对应关系正确就可以了。在 ode45()函数中，空矩阵[]表示使用默认的控制模板。 
>> 	b1=8/3; r1=10; s1=28; tn=100; %设置附加参数的值 
x0=[0;0;1e-10]; %数值微分方程初值 
[t,x]=ode45(@lorenz1,[0,tn],x0,[],b1,r1,s1); 
plot(t,x) %微分方程的数值求解，并绘制时间响应曲线

可以使用下面的语句绘制三维相空间曲线，得出的结果与前面是完全一致的，说明
两种方程描述方法是等效的。 >> plot3(x(:,1),x(:,2),x(:,3)); %打开新图形窗口，绘制三维相空间轨迹有了带有附加参数的微分方程 M函数后，就可以在不同参数值下求解微分方程，
而无须改变 lorenz1.m文件。例如，选择 Δ =2,λ =5,ν = 20，可以用下面的语句直接求出数值解，并分别得出如图 3-9和图 3-10所示的二维和三维图形。 
>> 	tn=100; x0=[0;0;1e-10]; b2=2; r2=5; s2=20; %另一组附加参数 [t2,x2]=ode45(@lorenz1,[0,tn],x0,[],b2,r2,s2); %求解方程 plot(t2,x2), figure; %绘制时间响应曲线

图 3-9新参数下 Lorenz方程的时间响应曲线
plot3(x2(:,1),x2(:,2),x2(:,3)); grid; %绘制三维相空间轨迹

图 3-10新参数下 Lorenz方程的三维相空间轨迹
由结果可知，如果混沌系统的参数发生了变化，得到的结果性质也可能完全不同。
如果微分方程比较简单，则可以采用匿名函数的形式描述该方程，无须使用附加参数的形式，因为匿名函数可以直接使用 MATLAB工作空间中的变量。这时，求解语句变成如下的形式，得到的结果也是完全一致的。 
>> 	b=8/3; r=10; s=28; %使用匿名函数则可以不用附加参数 
f=@(t,x)[-b*x(1)+x(2)*x(3); -r*x(2)+r*x(3); 

-x(1)*x(2)+s*x(2)-x(3)]; [t,x]=ode45(f,[0,tn],x0); %重新求解微分方程 plot(t2,x2), figure; %时域响应曲线的绘制 plot3(x2(:,1),x2(:,2),x2(:,3)); grid; %绘制三维相空间轨迹 

3.3.6避免附加参数的方法
在实际应用中，采用附加参数的方式并不是很方便，因为调用函数时，也需要对应地给出这些附加参数。对于简单问题，可以采用匿名函数直接描述。下面将介绍一种避免附加参数的方法。
考虑使用匿名函数为 MATLAB函数作一个接口函数，附加参数可以由 MAT-LAB工作空间直接传入 MATLAB函数，无须采用附加参数传递这些变量。下面举例说明。
　　例 3-11考虑例 3-10中的微分方程模型，该例中采用了附加参数描述来函数微分方程，构建了文件 lorenz1.m，试建立无须附加参数的模型描述与求解方法，重新求解原始的微分方程。　　解为避开附加参数，可以将这些参数的值先输入 MATLAB工作空间，然后利用匿名函数为该函数编写一个接口来描述原始的微分方程。 
>> 	b=8/3; r=10; s=28; %将参数输入 MATLAB工作空间
·68·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 
f=@(t,x)lorenz1(t,x,b,r,s) %为模型写一个匿名接口函数这样构造的匿名函数可以完全描述原始的微分方程模型，然后调用 ode45()函数直接求解，得出的结果与前面例子完全一致。 
>> 	tn=100; x0=[0;0;1e-10]; 
[t,x]=ode45(f,[0,tn],x0); plot(t,x) %求解并绘图 



3.4微分方程数值解的验证
前面通过例子说明，仿真算法和控制参数（例如相对误差限）选择不当，可能得到不可信甚至错误的结果。因此，在方程求解结束之后，应该对仿真结果进行检验。
在实际应用中，所求解的问题往往不存在解析解，怎么检验所得结果的正确性呢？下面介绍几种可以参考的验证思路。
（1）修改仿真控制参数，如可接受的误差限。将精度控制选项设置成一个更小的值，观察所得的结果是否和上次得出的结果完全一致，如果存在不能接受的差异，则应该考虑进一步减小误差限。为稳妥起见，也可以将这些控制选项直接设置为机器精度 eps，得出双精度数据结构下最可信的解。

（2）同样的问题选择不同的微分方程数值求解算法，可以交叉检验所得结果的正确性。 


3.4.1计算结果的验证
求解微分方程时，有时需要对求解算法及控制条件进行进一步设置，这可以通过 options变量进行修改。初始 options变量可以通过 odeset()函数获取，该变量是一个结构体变量，其中有众多成员变量，表 3-3中列出了常用的一些成员变量。修改该变量有两种方式，一是用 odeset()函数设置；二是直接修改 options的成员变量。例如，若想将相对误差设置成较小的 10.7，则需要给出 
options=odeset('RelTol',1e-7); %或更直观地采用下面的格式 
options=odeset; options.RelTol=1e-7;
下面通过例子演示微分方程数值解的验证方法。
　　例 3-12由图 3-7给出的曲线可以看出，结果比较粗糙，似乎原始求解结果的精度有问题，因此，应该考虑设置更严格的精度控制选项，得出精确的最终结果。　　解将相对误差限与绝对误差限设置为非常小的数值，例如 10.10，重新求解微分方程，得出的相平面曲线如图 3-11所示。可见，得出的曲线是光滑的，说明解的精度是比较高的，进一步降低误差限也无法得到更精确的结果。 
>> 	x0=[2; 3]; tn=10; 
f=@(t,x)[4*x(1)-2*x(1)*x(2); x(1)*x(2)-3*x(2)]; 

表 3-3常微分方程求解函数的控制参数表
参数名 参数说明  
RelTol 相对误差容许上限，默认值为 0.001（即 0.1％的相对误差）。在一些特殊的微分方程 
求解中，为了保证较高的精度，还应当再适当减小该值  
AbsTol 一个向量，其分量表示每个状态变量允许的绝对误差，其默认值为 10.6。可以自由 
设置其值，以改变求解精度  
MaxStep 求解方程时允许的最大步长  
Mass 微分代数方程中的质量矩阵，用于描述微分代数方程  
Jacobian 描述 Jacobi矩阵函数 df/dx的函数名。已知该 Jacobi矩阵，能加速仿真过程，得出 
所期望方程的解  
Events 事件响应属性，可以将其设置为事件的函数句柄  
OutputFcn 每步成功计算后自动调用自定义函数  

ff=odeset; ff.RelTol=1e-10; ff.AbsTol=1e-10; tic, [t,x]=ode45(f,[0,tn],x0,ff); toc %解方程 length(t), plot(x(:,1),x(:,2)) %绘制相平面曲线

图 3-11微分方程的相平面光滑曲线
本例的计算精度设置提升了 107倍，是否会因此显著增大计算量呢？上述代码执行耗时 0.02007s，计算点数增加到 3593，从耗时角度看并未显著增加运算量，所以值得用这样的方法得出问题的精确结果。
如果将误差限设置为机器精度 eps（实际上，所求解函数允许的最小误差限为 100eps），运算耗时仅增加到 0.0487s，计算点数增加到 19329。 >> ff.RelTol=eps; ff.AbsTol=eps; %设置误差限 tic, [t,x]=ode45(f,[0,tn],x0,ff); toc, length(t)
　　例 3-13试利用 ode45()函数重新求解例 3-9中的问题。　　解由于原始问题不存在解析解，所以检验得到的解并不是一件简单的事情。一种比较可行的方法是设置不同的精度控制参数，对方程进行求解，观察在不同的控制参数下得到的结果是否完全一致。在求解控制参数中，绝对误差限 AbsTol与相对误差限
·70·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 
RelTol是两个重要的控制参数，必要时可以将它们都设置成最苛刻的值 eps。如果在这样的控制精度下仍然不能得出方程的精确解，则说明 ode45()函数不适合求解该问题，只能寻求更有效的算法了。
可以将两个误差限设置成最苛刻的 eps，则可以由下面的语句重新求解方程，得出方程的数值解。 
>> 	y0=[0.994,0,0,0,-2.0015851063790825224,0]'; 
ff=odeset; ff.AbsTol=eps; ff.RelTol=eps; 
tic, [t,y]=ode45(@threebody,[0,40],y0,ff); toc 
plot3(y(:,1),y(:,2),y(:,3)), length(t) %绘制相空间轨迹

执行上面的代码，将得到警告信息“ Warning: RelTol has been increased to 
2.22045e-14.”（RelTol已经增加到 2.22045×10.14），说明实际使用的相对误差限没有达到预期的 eps，这样得出的结果已经是双精度数据结构下能得出的最精确的解了。绘制的相空间轨迹如图 3-12所示。仔细观察图中的曲线，发现在图中标记处曲线开始发散，若增大仿真终止时间，则最终的轨迹将发散，导致错误的结果。由 length()函数还可以测出求解的计算点数为 83485。 

图 3-12卫星的相空间轨迹（精确结果） 

3.4.2中间计算结果的动态处理
正常情况下，可以先用匿名函数或 MATLAB函数描述微分方程，然后调用 ode45()这类求解函数求解微分方程，一次性得出微分方程的数值解。如果求解过程中想显示中间结果，则可以通过设置控制选项中的 OutputFcn参数来实现。可以编写一个处理函数，以便在每一步成功计算后， MATLAB运行机制自动调用该函数处理中间结果。 
MATLAB还提供了 4个现成的函数用作常规处理，这 4个函数分别是
（1）@odeplot 自动绘制每个状态与时间的动态关系曲线； 
（2）@odephas2 自动绘制前两个状态变量的相平面动态曲线； 
（3）@odephas3 自动绘制前三个状态变量的相空间动态曲线； 

（4）@odeprint用数字的方式显示当前时刻和状态变量的值。
下面通过例子演示中间结果的处理方法。

　　例 3-14考虑例 3-7中的 Lorenz方程模型，试动态显示状态的相空间轨迹。　　解与例 3-7相比，只需修改控制选项中的 OutputFcn参数，将其设定为 @odephas3即可。另外，为保证计算精度，这里使用更严格的精度控制选项。运行以下程序，可以直接显示状态变量的动态相空间轨迹。 
>> f=@(t,x)[-8/3*x(1)+x(2)*x(3); 

-10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3)]; 

tn=100; x0=[0;0;1e-10]; 

ff=odeset('RelTol',100*eps,'AbsTol',100*eps,... 

'OutputFcn',@odephas3); %一次性修改多个控制选项 
[t,x]=ode45(f,[0,tn],x0,ff); %解方程并绘制动态相空间轨迹
为获得醒目的动态轨迹效果， MATLAB设置的默认轨迹线宽过大，是由 'Marker'选项产生的。如果不满意该显示效果，可以考虑自己编写中间曲线绘制的函数或直接修改底层的 odephas3.m文件。具体的方法是，在源程序编辑界面中打开该文件，找到下面的语句： 
ud.line = plot3(y(1),y(2),y(3),'-o','Parent',TARGET_AXIS);将'-o'选项改写成 '-'即可。如果需要，也可以相应地修改其他函数源程序，例如 odeplot.m、odephas2.m等。 

3.4.3更高精度的数值计算函数
从前面的例子可知，即使设置了很苛刻的精度控制选项，得到的仿真结果仍有可能存在较大误差，已经超出 ode45()这类函数的能力范围，所以考虑引入更精确的求解函数。
俄罗斯学者 Vasiliy Govorukhin开发了基于八阶七级的 Runge–Kutta变步长算法的求解函数 ode87() [21]，该算法的理论精度可达 o(h8)，每步求解过程调用 13次模型函数，所以求解效率远远高于 ode45()等函数，在某种意义上可以得到某些微分方程更高精度的解。该函数的调用格式与 ode45()等函数完全一致，注意不能采用过小的误差限，否则可能得到错误信息。
　　例 3-15试利用高精度求解函数重新求解例 3-9中的卫星轨迹绘制问题。　　解可以采用上述高精度的求解函数重新求解原始微分方程，然而得出的精度与前面完全一致，可能出现的现象也是一致的。因为在双精度数据结构下， ode45()函数已经达到了精度极限，即使使用高精度算法也不能得出更精确的结果。不同的是，本例采用的方法的计算点数大大减少，减少到 2511，只有 ode45()函数计算点数的 3％左右，故该算法的效率明显高于 ode45()。
·72·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 
>> 	ff=odeset; ff.AbsTol=1e-15; ff.RelTol=1e-15; 
y0=[0.994,0,0,0,-2.0015851063790825224,0]'; 
[t,y]=ode87(@threebody,[0,40],y0,ff); %求解方程 
plot3(y(:,1),y(:,2),y(:,3)), length(t) %绘制相空间轨迹

函数 ode87()是在双精度框架下编写的，不能脱离双精度框架来执行，即使采用了高精度的运算公式，也可能得不出更精确的结果，需要寻求更精确的方法。此外，由于这里采用了两个不同的求解函数，得出的结果是完全一致的，也可以通过这种方法验证结果的正确性。 

3.4.4计算步长与定步长显示
与前面编写的 ode_rk2()等函数不同， MATLAB自带的微分方程求解函数 ode45()与其他第三方函数（例如 ode87()等）均采用了变步长算法。这里将通过例子演示如何显示求解过程的步长情况，并介绍变步长算法中如何返回定步长输出的方法。
　　例 3-16试求解例 3-8中的微分方程，并观察各个时刻的计算步长。　　解函数 ode45()可以返回时间向量 t，对向量 t作差分运算（后项减前项）则可以得到每一步的步长值。差分运算可以通过 diff(t)实现，该函数得到的差分向量的长度比向量 t少一。所以，如果将误差限设置成 100倍的机器精度 eps，则可以由下面语句求解微分方程，并绘制出步长曲线，如图 3-13所示。可以看出，计算步长一直在变化，确保满足误差限的要求。 
>> 	x0=[2; 3]; tn=10; 
f=@(t,x)[4*x(1)-2*x(1)*x(2); x(1)*x(2)-3*x(2)]; 
ff=odeset; ff.RelTol=100*eps; ff.AbsTol=100*eps; 
[t,x]=ode45(f,[0,tn],x0,ff); %求解方程 
plot(t(1:end-1),diff(t)), min(diff(t)) %绘制步长曲线


图 3-13各个时刻的步长
从结果看，初始时刻的步长取得很小，最小的步长值为 1.859×10.4。
在介绍 ode45()函数调用格式时，求解范围用 tspan变元描述。如果 tspan为给定的时间向量 t，则 ode45()函数会采用变步长方法求出向量 t的每个子区间的微分方程的解，并返回向量各点处方程的解。即使 t为等间距的向量，也并不意味采用定步长算法求解微分方程。
　　例 3-17例 3-1给出了一个微分方程，并在后续例子中使用各种定步长算法求解该微分方程，得出的结论是，当步长 h比较大时，得到的方程的解误差比较大，甚至可能超出可以接受的范围。试利用 ode45()函数得出步长 h =0.2时各点处微分方程的解。　　解将 tspan设置成步距为 h =0.2的等间距向量 t，并设置苛刻的误差限，就可以得出 t向量各点处方程的解，在图中用 “o”记号标记出来，如图 3-14所示。可见，等间距时间点处的函数值精确地位于解析解曲线上，得出的最大误差为 2.6423×10.14。 
>> 	f=@(t,x)[x(2); -2*x(1)]; t0=0; tn=5; sq=sqrt(2); 
h=0.2; t=t0:h:tn; 
ff=odeset; ff.RelTol=100*eps; ff.AbsTol=100*eps; 
[t1,x1]=ode45(f,t,[1;0],ff); %求数值解 
sol=[cos(sq*t1), -sq*sin(sq*t1)]; norm(sol-x1,inf) 
plot(t,sol,t1,x1,'o') %绘制数值解与解析解曲线


图 3-14等间距时刻的状态响应图
注意，这里虽然采用了等间距的返回向量，但实际计算过程中采用的是变步长机制，所以得出的结果是精确的。
计算结果的定步长显示是有意义的，下面举例说明定步长显示的应用。
　　例 3-18考虑例 {2-19给出的二阶微分方程，试求出不同 λ下的微分方程解 x2(t)。 
′ 
x1(t)= x2(t) 
′ 
x2(t)= .λx1(t),x1(0) = 1,x2(0) = 0
　　解若直接使用 ode45()函数的调用格式，对每一个 λ求解微分方程，得到的解的

·74·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 
t向量可能各不相同，甚至长度也可能不同，利用这样的解很难绘制出不同 λ的曲面。考虑使用定步长的显示方式，选择不同的 λ值，可以绘制出 x2(t)的曲面，如图 3-15所示。 
>> t=0:0.2:10; Lam=0.2:0.2:3; H=[]; 

for lam=Lam %在不同参数下求解微分方程 
f=@(t,x)[x(2); -lam*x(1)]; 
[t1,x]=ode45(f,t,[1;0]); H=[H x(:,2)]; 

end 
surf(Lam,t,H) %绘制解的三维曲面 

t轴
α轴 
图 3-15 x2
............ 

3.4.5高阶非线性微分方程的求解实例
参考文献[22]提供的网站给出了很多微分方程初值问题数值解的测试实例，其中很多实例都是有应用背景的。这里用其中最简单直观的例子演示微分方程的求解。有兴趣的读者可以尝试该测试问题集中更复杂的例子。
　　例 3-19考虑下面的一阶显式微分方程组 [22]。 
.1.71y1(t)+0.43y2(t)+8.32y3(t)+0.0007 1.71y1(t) . 8.75y2(t) .10.03y3(t)+0.43y4(t)+0.035y5(t) 8.32y2(t)+1.71y3(t) . 1.12y4(t) 
y ′ (t)= 
.1.745y5(t)+0.43y6(t)+0.43y7(t) 
.280y6(t)y8(t)+0.69y4(t)+1.71y5(t) . 0.43y6(t)+0.69y7(t) 
280y6(t)y8(t) . 1.81y7(t) 
.280y6(t)y8(t)+1.81y7(t) 

且初值为 y0 = [1, 0, 0, 0, 0, 0, 0, 0.0057]T。试求出 0 . t . 321.8122时该微分方程的数值解。　　解采用 MATLAB求解，可以使用下面的命令，将微分方程用匿名函数描述出来，然后调用 ode45()函数求解。求解过程很快， 0.15s就可以得出方程的解，在比较严苛的
在不同 λ下的时域响应解 
. ............
精度要求下，总共计算出 41521个点。图 3-16给出了这些状态变量的时域响应曲线，与文献 [22]给出的结果完全一致。 
>> f=@(t,y)[-1.71*y(1)+0.43*y(2)+8.32*y(3)+0.0007; 
1.71*y(1)-8.75*y(2); 
-10.03*y(3)+0.43*y(4)+0.035*y(5); 
8.32*y(2)+1.71*y(3)-1.12*y(4); 
-1.745*y(5)+0.43*y(6)+0.43*y(7); 
-280*y(6)*y(8)+0.69*y(4)+1.71*y(5)-0.43*y(6)+0.69*y(7); 
280*y(6)*y(8)-1.81*y(7) 
-280*y(6)*y(8)+1.81*y(7)]; 

y0=[1;0;0;0;0;0;0;0.0057]; %微分方程初值 
ff=odeset; ff.RelTol=100*eps; ff.AbsTol=100*eps; 
tic, [t,x]=ode45(f,[0,321.8122],y0); toc 
plot(t,x), xlim([0,5]) %微分方程的时域响应 


图 3-16测试例子的时域响应解（0 . t . 5局部放大图）还可以由下面的语句绘制出步长的曲线，如图 3-17所示。在计算初始阶段，最小的步长为 2.9379×10.5，而 t比较大时可能采用较大计算步长，在保证精度的前提下加快计算速度，所以这里给出的求解算法与工具是有效的。
在原始文献中，该方程被认定为刚性微分方程（有关刚性微分方程的内容将在第 5章进一步介绍），难以用普通方法求解。而这里采用 MATLAB直接求解，过程是直观的，在极短的时间内就可以得出问题在严苛条件下的精确解。由此可见， MATLAB在求解微分方程时功能非常强大。


本章习题 
3.1试利用代数方程求解的方法验证四阶 Runge–Kutta算法与 Gill算法。
·76·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解

图 3-17求解过程的计算步长 
3.2试使用 MATLAB函数编写 5阶和 6阶 Adams–Bashforth算法，并评价其有效性。
（1）5阶 Adams–Bashforth算法： 
[
h 
xk+1 = xk + 1901f(tk, xk) . 2774f(tk.1, xk.1) + 2616f(tk.2, xk.2)
720 
]
.1274f(tk.3, xk.3) + 251f(tk.4, xk.4)
（2）6阶 Adams–Bashforth算法： 
[
h 
xk+1 = xk + 4277f(tk, xk) . 7923f(tk.1, xk.1) + 9982f(tk.2, xk.2)
1440 
]
.7298f(tk.3, xk.3) + 2877f(tk.4, xk.4) . 475f(tk.5, xk.5)
3.3极限环是非线性常微分方程中一种常见的现象，对某些非线性微分方程来说，不论初始状态为何值，微分方程的相轨迹都将稳定在一条封闭的曲线上，该曲线称为微分方程的极限环。试绘制出下面微分方程的极限环，并对不同初值验证微分方程的相平面曲线确实收敛于极限环。 
{ 
x ′ (t)= y(t)+ x(t)(1 . x2(t) . y2(t)) y ′ (t)= .x(t)+ y(t)(1 . x2(t) . y2(t)) 
3.4考虑下面给出的非线性微分方程。文献 [23]指出，该方程具有多个极限环， r = 1/(n丘)，n =1, 2, 3, ···。试用数值方法求解此方程并观察多极限环的情况，其中， f(r)= r2 sin(1/r). 。 (√ )
. x ′ (t)= .y(t)+ x(t)fx2(t)+ y2(t)
(√ )
. y ′ (t)= x(t)+ y(t)fx2(t)+ y2(t)
3.5考虑如下的著名的 R.ssler.化学反应方程组。 
.x ′ (t)= .y(t) . z(t)
. 
y ′ (t)= x(t)+ ay(t)
.
. z ′ (t)= b +(x(t) . c)z(t)选定 a = b =0.2，c =5.7，且 x(0) = y(0) = z(0)，绘制仿真结果的三维相空间
轨迹，并得出其在 xy平面上的投影。建议将 a、b和 c作为附加参数，若设 a =0.2， b =0.5，c = 10时，绘制出状态变量的二维图和三维图。 
3.6考虑下面的微分方程组 [24]。 
y ′ (t)= tan .(t) gsin .(t)Ψv2(t) 
v ′ (t)= .
v(t) cos .(t) 
. .. .. Chua其中，为电路的二极管分段线性特性，即f()x
. ′ (t)= .g/v2(t)
其中，g =0.032，Ψ =0.02，初值为 y(0) =0，v(0) = 0.5，.(0)分别取 0.3782和 
9.7456，试求微分方程的数值解。 
3.7 Chua电路方程是混沌理论中经常提到的微分方程 [25]。 
x ′ (t)=叫[ y(t) . x(t) . f(x(t))] y ′ (t)= x(t) . y(t)+ z(t) z ′ (t)= .Δy(t) . Ψz(t)
{
..
..
. {
1 
f(x)= bx +(a . b)(|x +1|.|x . 1|)
2 
且 a<b< 0。试编写出 MATLAB函数描述该微分方程，并绘制出叫 = 15， Δ = 20，Ψ =0.5，a = .120/7，b = .75/7，且初值为 x(0) = .2.121304， y(0) = .0.066170，z(0) = 2.881090时的相空间曲线。 
3.8试求解下面微分方程的数值解 [19]。 
y ′ (x)= .2xy(x) ln z(x) 
z ′ (x)=2xz(x) ln y(x)
已知初值为 y(0) = e，z(0) = 1，且该微分方程的解析解为 y(x)= ecos x 2，z(x)= esin x 2。试对比各种算法在求解该问题时的速度与精度。 
3.9考虑下面的微分方程 [26]： 
w ′ (t)=叫Kw + Δy(t) . Ψx(t)w(t) .叫w(t)z(t) x ′ (t)= .x(t)+ Δy(t) . Ψx(t)w(t) . y(t)z(t)/Ka y ′ (t)= x(t) . Δy(t)+ Ψx(t)w(t) . y(t)z(t)/Ka z ′ (t)=叫Kw + x(t)+叫w(t)z(t) . y(t)z(t)/Ka
其中，Ka = 31 × 10.8，Kw =3.25 × 10.18，叫 = 55.5 × 106，Δ = 10 × 10.6， Ψ =1.11×106。若已知各个状态变量的初值均为 0，试求解该微分方程。 
3.10 Lotka–Volterra捕食模型方程如下，试求解该微分方程，并绘制响应曲线，其中， x(0) = 2，y(0) = 3。 x ′ (t)=4x(t) . 2x(t)y(t) 
y ′ (t)= x(t)y(t) . 3y(t) 
. .
. .
. 
·78·薛定宇教授大讲堂（卷 V）：MATLAB微分方程求解 
3.11 Van der Pol方程的数学形式如下，试绘制出不同 μ时，状态变量 x1(t)的时域响应曲面，其中，x1(0) = .0.2，x2(0) = .0.7。 
′ 
x1(t)= x2(t) 
′ 2
x2(t)= .μ(x1(t) . 1)x2(t) . x1(t) 
3.12化学反应动力学方程如下 [27]： 
{
. .. 
′ (t)= s
′ (t)= h
′ (t)= x(t) . z(t)
[
[]
qy(t) . x(t)y(t)+ x(t)(1 . x(t)). qy(t) . x(t)y(t) . pz(t)
x 
]
y
.. 
z 
其中， h 
20，且已知初值为 x(0)
8.333，p =0.3，q =0.01，s 
= 
1，y(0) = 2，
= = 
z(0) = 0.6，若 tn = 60，试求解该微分方程。 
3.13 FitzHugh–Nagumo模型 [27]。 
[
v ′ (t)= a′ (t)= .
{
]
]
[
v(t)+ w(t) . v3(t) . 0.5v(t) . b + cw(t)
/a
w 
其中， a =3， 
试用三种方法
b =0.7，c =0.8，且已知初值 v(0) = 0.5，w(0) = 0.2，
{
. .. .. 
求解该微分方程。 
{
3.14考虑下面的微分方程模型，试选择合适的终止时间求解微分方程，以便能观察到相平面轨迹的全过程。 
′ 2	2
x	1(t)= .x2(t) . 10x1(t)+5x1(t)x2(t)+ x2(t) ′ 2
x2(t)= x1(t)+ x1(t) . 25x1(t)x2(t)
初值为 x1(0) = .0.0914，x2(0) = .0.1075。 

3.15试求解下面的微分方程模型 [5]。 
′ 
I1y1(t)=(I2 . I3)y2(t)y3(t) 
I2y (t)=(I3 . I1)y3(t)y1(t)
2′
′ 
I3y3(t)=(I1 . I2)y1(t)y2(t)+ f(t)
其中，f(x)为分段函数 
f(t)=0.25 sin2 t, 3丘 . t . 4丘 
0,其他

已知常数 I1 =0.5，I2 =2，I3 =3，y1(0) = 1，y2(0) = 0，y3(0) = 0.9。 

3.16 Hodgin–Huxley神经模型 [27]。 
i = CV ′ (t)+ iNa(t)+ iK (t)+ iL(t)
. ..
))
(
v . VNav . VK
其中， 
iNa(t)= gNam
(
3(t)h(t)4(t)
..(
iK (t)= gK niL(t)= gLv . VL
) 
已知常数 i = 10，C =1，VNa = 115，VK = 12，VL = 10.599，gNa = 120，gK = 36， gL =0.3。另外，已知微分方程 
m ′ (t)= am(1 . v) . bmm(t) n ′ (t)= an(1 . v) . bnn(t) h ′ (t)= ah(1 . v) . bhh(t)
其中， am =0.1(25 . v) ,bm =4e.v/18 
e0.1(25.v) . 1 
0.01(10 . v) 

an = ,bn =0.125e.v/80 
e0.01(10.v) . 1 
ah =0.07e.v/20,bh =1 
e0.1(30.v) +1
且已知初值 V (0) = 0，m(0) = n(0) = h(0) = 0.5，若 v在 (.80, 80)取值，试求解不同 v值下的微分方程，其中 t → (0, 50)。 
3.17参考文献 [22]提供的网站给出了很多微分方程数值解的基准测试问题，例 3-19只
. .. .. 
演示了其中最简单直观的微分方程的求解。读者可以访问该网站并下载相关的问题，探讨在 MATLAB下如何求解，效果如何，求解效率如何。