利用 MATLAB 求解常微分方程数值解
目录
1. 内容简介 .................................................................................................................................... 1
2.
Euler Method(欧拉法)求解 ....................................................................................................... 1
2.1. 显式 Euler 法和隐式 Euler 法 ....................................................................................... 2
2.2. 梯形公式和改进 Euler 法 ..............................................................................................3
2.3.
Euler 法实用性 ...............................................................................................................5
3.
Runge-Kutta Method(龙格库塔法)求解 ................................................................................... 6
Runge-Kutta 基本原理 ...................................................................................................6
3.1.
3.2. MATLAB 中使用 Runge-Kutta 法的函数 ....................................................................... 8
4. 使用 MATLAB 求解常微分方程................................................................................................ 8
4.1. 使用 ode45 函数求解非刚性常微分方程 ................................................................... 8
4.2. 刚性常微分方程 ............................................................................................................9
5. 总结 .......................................................................................................................................... 10
参考文献 .......................................................................................................................................... 11
附录 .................................................................................................................................................. 12
1.
2.
3.
4.
显式 Euler 法数值求解................................................................................................12
改进 Euler 法数值求解................................................................................................12
四阶四级 Runge-Kutta 法数值求解............................................................................13
使用 ode45 求解 ..........................................................................................................14
1. 内容简介
把《高等工程数学》看了一遍,增加对数学内容的了解,对其中数值解法比较感兴趣,
这大概是因为在其它各方面的学习和研究中经常会遇到数值解法的问题。理解模型然后列出
微分方程,却对着方程无从下手,无法得出精确结果实在是让人难受的一件事情。
实际问题中更多遇到的是利用数值法求解偏微分方程问题,但考虑到先从常微分方程下
手更为简单有效率,所以本文只研究常微分方程的数值解法。把一个工程实际问题弄出精确
结果远比弄清楚各种细枝末节更有意思,因此文章中不追求非常严格地证明,而是偏向如何
利用工具实际求解出常微分方程的数值解,力求将课程上所学的知识真正地运用到实际方程
的求解中去,在以后遇到微分方程的时候能够熟练运用 MATLAB 得到能够在工程上运用的结
果。
文中求解过程中用到 MATLAB 进行数值求解,主要目的是弄清楚各个函数本质上是如何
对常微分方程进行求解的,对各种方法进行 MATLAB 编程求解,并将求得的数值解与精确解
对比,其中源程序在附录中。最后考察 MATLAB 中各个函数的适用范围 ,当遇到实际工程
问题时能够正确地得到问题的数值解。
2. Euler Method(欧拉法)求解
Euler 法求解常微分方程主要包括 3 种形式,即显式 Euler 法、隐式 Euler 法、梯形公式
法,本节内容分别介绍这 3 种方法的具体内容,并在最后对 3 种方法精度进行对比,讨论
Euler 法的实用性。
本节考虑实际初值问题
1
使用解析法,对方程两边同乘以 得到下式
两边同时求积分并采用分部积分得到解析解:
本节后面将对此方程进行求解,并与精确解进行对比,分析 Euler 的可行性。
2.1. 显式 Euler 法和隐式 Euler 法
显式和隐式 Euler 法都属于一阶方法,显式 Euler 法的迭代公式简单,如下所示:
对过上述公式对式 进行迭代,其中步长
,计算
之间的数值,迭代求解
的 MATLAB 程序见附录 ,能够得出精确解和数值解的图像,如图 所示。
图 2.1 显式 Euler 法精确解和数值解图像
2
从图 2.1 中可以看出,显式 Euler 法在斜率很大的时候存在非常大的误差。本质上是 Euler
法只计算了每一步差值中的一阶部分,由 Taylor 级数可知:
当公式中的二阶导数较大时就会产生明显的偏差,同时迭代过程中由于使用到上一部的
结果,误差会在迭代中传播,因此这种 Euler 法在实际中是无法使用的,但是却给求解微分
方程数值解提供了好的开始。
另外一种 Euler 法是隐式 Euler 法,其迭代公式是
,它并
没有解决上面所说的问题,同时它的计算更加繁琐,对于无法化简成显示迭代的公式时还需
要用迭代法求解非线性方程。
为了解决上面的方法,就需要提高迭代公式中计算差值的阶数,下面介绍了梯形法和改
进 Euler 法,它们都是二阶方法。
2.2. 梯形公式和改进 Euler 法
梯形公式以及改进 Euler 法都属于二阶方法,下面证明它是二阶方法,使用两次 Taylor
公式,将
和
展开:
将
得到
从上式可以看出,梯形法的局部截断误差的主要部分是
,是关于步长的三次
3
式,这说明了梯形法取到了差值中的二次项,因此梯形法是二阶方法。
从上面可以得到梯形的迭代公式:
但是上式并不容易计算,因为上式中的 为带求量,当
无法化成显式形
式时,需要对上式进行迭代求解。因此梯形公式不易通过计算机编程求解,实际上改进的
Euler 法更容易求解。
改进 Euler 法迭代公式先通过显式 Euler 法求出一个估计值 ,通过这个估计值来计算
,然后方程 就变成了显式方程,从而可以得到修正值 ,改进 Euler 法更
适合计算机编写程序,同样解决初值问题
,详细 MATLAB 程序见附录 2,得到的对比图
像如图 2.2 所示。
图 2.2 改进 Euler 法精确解与数值解对比
4
由于改进 Euler 法用来求解
的并不是精确解,所以得到的导数会有一定误
差,因此改进 Euler 法的实际局部截断误差不仅仅是
。
2.3. Euler 法实用性
从图 2.1 可以看出来一阶方法精确度非常差,基本上是无法用到实际工程中的,因此显
式和隐式 Euler 法只是提供一种对微分方程求解的思想。
从图 2.2 中得到的数值解相对图 2.1 已经有了明显的改善,但是对于精度要求较高的工
程问题,梯形法和改进 Euler 法这样的二阶方法同样是不满足要求的。但通过这两个方法,
了 解 到 提 高 方 法 取 差 值 中 的 更 高 阶 数 项 能 够 达 到 提 高 精 度 的 目 的 , 后 面 内 容 中 的
Runge-Kutta 法就是由此思想而来。
将细节部分放大更能够看出两种方法的精度,如图 2.3 所示(左图是 Euler 法,右图是改
进 Euler 法)。
图 2.3 两种方法细节部分
5
3. Runge-Kutta Method(龙格库塔法)求解
这一节的目的是弄清楚 Runge-Kutta 法(后面简称 R-K 法)的基本原理,并弄清楚 MATLAB
中 ode 系列适合解决哪一种工程问题,希望达到的目的对于任何一个工程上提出的常微分方
程问题都能够利用 MATLAB 求得误差在能够接受范围内的数值解。
3.1. Runge-Kutta 基本原理
MATLAB 中使用 R-K 法的函数有 ode23 和 ode45,其中 ode 是 ordinary differential equation
的缩写,其中 ode23 表示使用的是 2 阶 3 级 R-K 法,ode45 使用的是 4 阶 5 级 R-K 法。
将局部截断误差表示如下式:
将上式中的方法称作 阶 级 Runge-Kutta 法。其中:
下面将构造出四阶四级经典 R-K 法,并通过 MATLAB 编程实现对初值问题
求解。
四级四阶 R-K 法的本质是用四项级数来取得 Taylor 展开的前面四阶项,局部截断误差首
项是步长的五次方,即:
6