MATLAB求解常微分方程数值解文档格式.docx

MATLAB求解常微分方程数值解文档格式.docx

《MATLAB求解常微分方程数值解文档格式.docx》由会员分享，可在线阅读，更多相关《MATLAB求解常微分方程数值解文档格式.docx（13页珍藏版）》请在冰豆网上搜索。

对过上述公式对式

进行迭代，其中步长

，计算

之间的数值，迭代求解的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-KuttaMethod（龙格库塔法）求解

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

3.1.Runge-Kutta基本原理

MATLAB中使用R-K法的函数有ode23和ode45，其中ode是ordinarydifferentialequation的缩写，其中ode23表示使用的是2阶3级R-K法，ode45使用的是4阶5级R-K法。

将局部截断误差表示如下式：

将上式中的方法称作

阶

级Runge-Kutta法。

其中：

下面将构造出四阶四级经典R-K法，并通过MATLAB编程实现对初值问题

求解。

四级四阶R-K法的本质是用四项级数来取得Taylor展开的前面四阶项，局部截断误差首项是步长的五次方，即：

从上述表达式中可以看出，实际上R-K法是希望通过线性组合来近似Taylor公式中的

阶导数组合，这样能够避免求解函数的

阶导数，从而大量减小了计算量。

实际计算中将

用

来表示，

则利用二元Taylor公式展开，舍去

，并将同次数项系数相等，得到经典四阶四级R-K法的公式：

上面从一阶到二阶的数值求法让误差有很大的改观，利用上述公式写出四阶四级的MATLABA程序（见附录3），得到的图像如图3.1所示，图中解析解和数值解基本上重叠了，并且在

处的误差已经小于

了，而显式Euler法的误差为

数量级，二阶改进Euler法误差在

数量级。

图3.1四阶四级R-K法数值解和精确解对比

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

MATLAB中利用R-K法求解常微分方程的函数主要是ode45，它用来求解非刚性常微分方程，是工程计算中首选函数。

另外的ode23是2阶3级R-K法函数，它可以用来求解轻微刚性方程，其精度相对较低，但是效率较高。

它们都是显式的R-K法，另外一个函数式ode23tb，它是隐式的R-K法可以用来求解非刚性的常微分方程。

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

MATLAB中求解常微分方程共有7个函数，表4-1列举出它们的使用围以及简介。

解算器

问题类型

精确度

适用条件

ode45

非刚性

中等

大多数情况下，最先尝试的求解器

ode23

低

低精度，稍微刚性问题

ode113

低-高

适用于对精度要求高的问题

ode15s

刚性

低-中等

精度一般，求解刚性问题

ode23s

低精度，刚性问题

ode23t

轻微刚性

ode23tb

表4-1MATLAB中求解常微分方程函数

通过表4-1了解到，当遇到一个常微分方程问题，首先使用ode45对其进行求解，如果求解非常慢或者无法求解，那么问题很可能是刚性问题，需要使用ode15s进行求解。

对于精度要求非常高的问题可使用ode113尝试进行求解，对于其他的刚性问题，需要具体研究使用哪一个函数。

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

ode45函数实现了变步长四阶五级Runge-Kutta-Felhberg算法（简称RKF法）。

调用格式如下所示。

公式中

是函数的句柄，

是求解区间的起点和终点，

是初始值，后面的

都是各种控制求解细节的选项。

其中的

可以是匿名函数、子函数、嵌套函数、单独M文件等形式。

实际上ode45能够求解常微分方程组，其中的

使用状态变量表示方法，左侧是各个状态变量的导数，右侧是方程。

ode45可以解决高阶常微分方程的数值求解问题，基本思想是将高阶通过变量替换变成低阶，然后列成下面所示的微分方程组，然后进行求解。

对于处置初值

问题，使用MATLAB求解（代码见附录4）如图3.2所示，其精度与上面自己编写的程序差不多，在

处的误差小于

图3.2使用ode45求解细节

4.2.刚性常微分方程

如上所述，ode45函数能够求解的是非刚性问题，而刚性（stiff）问题使用ode45求解会产生求解缓慢，或者完全无法求解的情况，这种情况可以判断这个方程是刚性的。

可以考虑使用ode23或者ode23s函数进行求解。

刚性常微分方程指的是其特征根相差很大，但似乎没有非常有效判断一个方程是刚性还是非刚性的办法，需要实际求解的时候通过不同的函数进行求解比较。

5.总结

文本主要研究了如何求常微分方程的数值解，从开始的一阶显式Euler法，到四阶R-K法，精度提高了很多。

同时通过编写MATLAB代码来对具体初值问题进行求解，对课本上将的机制有了更加深入的了解，对具体算法过程有了深入的认识，从中学到了很多。

另外文中总结了MATLAB中求解常微分方程数值解函数的用法，在面对不同的工程问题的时候，怎样去挑选正确的函数来求解具体问题，更多的是学会了如何去解决实际问题。

实际问题中更多是由偏微分方程进行描述的，由于时间和个人数学水平原因，文中没有对偏微分方程数值求解进行研究，这是文章的不足，有待在后续的学习当中对其进行解决。

参考文献

[1]于寅.高等工程数学[M].:

华中科技大学.2012

[2]红.数值分析（第2版）[M].:

华中科技大学.2010

[3]胡建伟,汤怀民.微分方程数值方法（第二版）[M].:

科学.2007

[4]薛定宇,.高等应用数学问题的MATLAB求解（第二版）[M].:

清华大学.2008

附录

1.显式Euler法数值求解

clear,clc%y1中存储0~5中的数值解y2中存储-5~0中的数值解

dyx=（x,y）（-2*y-4*x）;

%方程右边函数

fx=（x）（exp（-2*x）-2*x+1）;

%解析解

y1

（1）=2;

y2

（1）=2;

h=0.1;

i=2;

forx=0:

h:

5-h

y1（i）=y1（i-1）+h*dyx（x,y1（i-1））;

i=i+1;

end

-h:

-5+h

y2（i）=y2（i-1）-h*dyx（x,y2（i-1））;

i=i+1;

x=-5:

5;

y=fx（x）;

ye=[fliplr（y2（2:

end））2y1（2:

end）];

plot（x,y,x,ye,'

--'

）;

legend（'

\fontsize{12}\fontname{楷体}精确解'

...

'

\fontsize{12}\fontname{楷体}数值解'

'

interpreter'

latex'

text（x（8）+0.1,y（8）,'

$\leftarrowe^{-2x}-2x+1$'

fontsize'

15）;

2.改进Euler法数值求解

clear,clc%y1中存储0~5