MATLAB求解常微分方程数值解

利用MATLAB求解常微分方程数值解

目录

1. 内容简介 (1)

2. Euler Method(欧拉法)求解 (1)

2.1. 显式Euler法和隐式Euler法 (2)

2.2. 梯形公式和改进Euler法 (3)

2.3. Euler法实用性 (4)

3. Runge-Kutta Method(龙格库塔法)求解 (5)

3.1. Runge-Kutta基本原理 (5)

3.2. MATLAB中使用Runge-Kutta法的函数 (7)

4. 使用MATLAB求解常微分方程 (7)

4.1. 使用ode45函数求解非刚性常微分方程 (8)

4.2. 刚性常微分方程 (9)

5. 总结 (9)

参考文献 (11)

附录 (12)

1. 显式Euler法数值求解 (12)

2. 改进Euler法数值求解 (12)

3. 四阶四级Runge-Kutta法数值求解 (13)

4.使用ode45求解 (14)

1.内容简介

把《高等工程数学》看了一遍，增加对数学内容的了解，对其中数值解法比较感兴趣，这大概是因为在其它各方面的学习和研究中经常会遇到数值解法的问题。理解模型然后列出微分方程，却对着方程无从下手，无法得出精确结果实在是让人难受的一件事情。

实际问题中更多遇到的是利用数值法求解偏微分方程问题，但考虑到先从常微分方程下手更为简单有效率，所以本文只研究常微分方程的数值解法。把一个工程实际问题弄出精确结果远比弄清楚各种细枝末节更有意思，因此文章中不追求非常严格地证明，而是偏向如何利用工具实际求解出常微分方程的数值解，力求将课程上所学的知识真正地运用到实际方程的求解中去，在以后遇到微分方程的时候能够熟练运用MATLAB得到能够在工程上运用的结果。

文中求解过程中用到MATLAB进行数值求解，主要目的是弄清楚各个函数本质上是如何对常微分方程进行求解的，对各种方法进行MATLAB编程求解，并将求得的数值解与精确解对比，其中源程序在附录中。最后考察MATLAB中各个函数的适用范围，当遇到实际工程问题时能够正确地得到问题的数值解。

2.Euler Method(欧拉法)求解

Euler法求解常微分方程主要包括3种形式，即显式Euler法、隐式Euler法、梯形公式法，本节内容分别介绍这3种方法的具体内容，并在最后对3种方法精度进行对比，讨论Euler法的实用性。

本节考虑实际初值问题

使用解析法，对方程两边同乘以得到下式

两边同时求积分并采用分部积分得到解析解：

本节后面将对此方程进行求解，并与精确解进行对比，分析Euler的可行性。

2.1.显式Euler法和隐式Euler法

显式和隐式Euler法都属于一阶方法，显式Euler法的迭代公式简单，如下所示：

对过上述公式对式进行迭代，其中步长，计算之间的数值，迭代求解

的MATLAB程序见附录，能够得出精确解和数值解的图像，如图所示。

图2.1 显式Euler法精确解和数值解图像

从图2.1中可以看出，显式Euler法在斜率很大的时候存在非常大的误差。本质上是Euler 法只计算了每一步差值中的一阶部分，由Taylor级数可知：

当公式中的二阶导数较大时就会产生明显的偏差，同时迭代过程中由于使用到上一部的结果，误差会在迭代中传播，因此这种Euler法在实际中是无法使用的，但是却给求解微分方程数值解提供了好的开始。

另外一种Euler法是隐式Euler法，其迭代公式是，它并没有解决上面所说的问题，同时它的计算更加繁琐，对于无法化简成显示迭代的公式时还需要用迭代法求解非线性方程。

为了解决上面的方法，就需要提高迭代公式中计算差值的阶数，下面介绍了梯形法和改进Euler法，它们都是二阶方法。

2.2.梯形公式和改进Euler法

梯形公式以及改进Euler法都属于二阶方法，下面证明它是二阶方法，使用两次Taylor 公式，将和展开：

将得到

从上式可以看出，梯形法的局部截断误差的主要部分是，是关于步长的三次式，这说明了梯形法取到了差值中的二次项，因此梯形法是二阶方法。

从上面可以得到梯形的迭代公式：

但是上式并不容易计算，因为上式中的为带求量，当无法化成显式形

式时，需要对上式进行迭代求解。因此梯形公式不易通过计算机编程求解，实际上改进的Euler法更容易求解。

改进Euler法迭代公式先通过显式Euler法求出一个估计值，通过这个估计值来计算

，然后方程就变成了显式方程，从而可以得到修正值，改进Euler法更

适合计算机编写程序，同样解决初值问题，详细MATLAB程序见附录2，得到的对比图

像如图2.2所示。

图2.2 改进Euler法精确解与数值解对比

由于改进Euler法用来求解的并不是精确解，所以得到的导数会有一定误

差，因此改进Euler法的实际局部截断误差不仅仅是。

2.3.Euler法实用性

从图2.1可以看出来一阶方法精确度非常差，基本上是无法用到实际工程中的，因此显式和隐式Euler法只是提供一种对微分方程求解的思想。

从图2.2中得到的数值解相对图2.1已经有了明显的改善，但是对于精度要求较高的工程问题，梯形法和改进Euler法这样的二阶方法同样是不满足要求的。但通过这两个方法，了解到提高方法取差值中的更高阶数项能够达到提高精度的目的，后面内容中的Runge-Kutta法就是由此思想而来。

将细节部分放大更能够看出两种方法的精度，如图2.3所示(左图是Euler法，右图是改进Euler法)。

图2.3 两种方法细节部分

3.Runge-Kutta Method(龙格库塔法)求解

这一节的目的是弄清楚Runge-Kutta法(后面简称R-K法)的基本原理，并弄清楚MATLAB 中ode系列适合解决哪一种工程问题，希望达到的目的对于任何一个工程上提出的常微分方程问题都能够利用MATLAB求得误差在能够接受范围内的数值解。

3.1.Runge-Kutta基本原理