基于MATLAB的微分方程数值解法_冯元珍.docx

《基于MATLAB的微分方程数值解法_冯元珍.docx》由会员分享，可在线阅读，更多相关《基于MATLAB的微分方程数值解法_冯元珍.docx（3页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、科技信息 计算机与网络 基于 MATLAB的微分方程数值解法 南京人口管理干部学院基础部冯元珍南京医科大学屠小明 L摘要规实生活屮的许多问题都可以通过微分方程的形式进行表示，因此微分方程的求解具有很大的实际意义。本 文介绍了 MATLAB软件在微分方程数值求解屮的应用，并通过一个实际的例子对其具体应用进行了探讨。 咲键词 jMATLAB微分方程数值解法 1引言 微分方程是数学科学联系实际问题的主要桥梁之一，它是含有未知函数及其导数的方程。如果未知函数的自变量是一个， 称为常微分方程 （ Ordinary Differential Equations, ODE);自变量多于一个，称为偏微分方程

2、（ Partial Differential Equations, PDE)。 在科学研究和工程计算中碰到的许多微分方程，根本不存在解析解，或者求解析解的代价很大，求解过程 过于复杂，在这种情况下，我们只能借助于数值计算来求方程的数值解。 2常用的微分方程数值解法 常微分方程的求解问题大体上可以分为：初值问题和边值问题。相比较而言，求解边值问题的难度更大一些，必须根据 问题所在的背景和方程特征，具体问题具体解决。对于初值问题 ，常用的解法冇：欧拉方法，泰勒级数展开方法，龙格 -库塔 方法，外推方法等；边值问题的常用解法有：打靶法，有限差分方法，配置法等。对于偏微分方程，数值算法的实施也必须根

3、据方程的类型特征来进行选取，常用的方法有：半离散方法，全离散方法，有限差分方法等等 m。上述所有求解微分方程的数 值解法本质上都是通过一定的迭代过程来求解方程的可以接受的解。 对于工程中碰到的实际问题，直接采用上述方法求解的工作量通常都比较大，计算难度也很大。应用数学软件求解可以 在很大程度上提高计算效率 ,节约计算时间。关于微分方程的数值解法， 目前有许多比较成熟的数学软件，如 :MATLAB, 0DEPACK, NETLIB等。下面我们仅就 MALAB在微分方程数值解法方面的应用进行介绍。 3 MATLAB软件在微分方程数值求解中的应用 MATLAB是 MATRIX LABORATORY的

4、缩写，具有强大的科学计算、图形显示和程序设计功能。该软件提供了求解广义微分方 程初值和边值问题的完整的指令。 对于常微分方程初值问题， MATLAB提供了诸如 ode23, ode45, ode23t, odl5s, ode23tb等指令 t21。 其中 ode23和 ode45 分别采用普通 2-3阶 Runge-Kutta法和 4-5阶 Runge-Kutta法来求解常微分方程的初值问题，后面的几个指令分别用梯形法 求解适度刚性常微分方程，变阶法和低阶法求解刚性常微分方程。 de45因为计算精度相对较高而成为大多数场合的优先选择 算法。 MATLAB软件在求解常微分方程初值问题吋，对上述命

5、令采用了相同的格式。求解 ODE初值问题比较完整的调川格式如 下： /, Y, Te, Ye, Ie = solver ( fun, tspcin, Y 0, options, p, p2, ) (l) 其中 xo/ver是上述用来求解微分方程的指令，如 de23等， /im是需要求解的函数文件名 ， /5/7CW是方程求解对应的吋间 区间，F0以向量的形式存放问题的初始值， ;?/似用来设置算法的参数， pi,是向调用函数 /训传递的参数， /是 求得数值解对应自变量的数据列向量， 1T是以矩阵的形式对应存放每个计算时刻的数值解。另外三个输出宗量是在设置 op/bm的亊件属性以后对应的输出。

6、公式 （ 1)在使用时对应的参数相对较多，调用比较繁琐，在对算法的采 ;0和计钉的精 度不做特殊耍求时，通常会采爪下列简化公式： 413 科技信息 计算机与网络 /, Y = solver funjspan, 70) (2) 与式 （ 1)相比，式 （ 2)缺失的参数在实际计算时都采用系统默认的参数。 对于 ODE的边值问题，在假定解唯一的前提下， MATLAB提供了一个求解函数 bvp4c， 求得解的相对精度比较均匀。在应 用 MATLAB求解 ODE边值问题时，首先需要把待解的问题转化为标准的边值问题，即将边界条件写成形如： 0 的形式，然后再为期望解指定一个初始猜测。相关指令为 ： we

7、sM = int(x, V, , X是边界区间 (/, 上 的初始网格，通过单调上升、等距的向量来进行表示， v是对应网格点上解的初始猜测值，是待定的未知参数。 初始网格确定好之后，就可以利用命令 solve - bvp4c( fun, bcfun, mesh, options, p, p2, ) (3) 来求解方程的近似解：其中厂 = &/丨 2(,_,厂训训祕 /， /)1， ._)， &； /是计算边值条件函数， 1， 76分別表示左 右边界的列向量，次序不可以进行交换： res是 边 界 条 件 满 足 时 的 残 差 形 成 的 列 向 量 。 最 后 采 用 命 令 va/we=6

8、v/?va$o/vgx)就可以求得给定积分区间内任意点 x的解值。如果求解的结果不太满意的话，可以修改式 （ 3) 中的选项来改进求解的质量，具体的修改过程在此不再详述，第四个及其后的参数在简单调用指令中可以略去不写。 关于 PDE的数值解法， MATLAB软件只针对含有一个时间变量和空间变量的抛物型 PDE和椭圆型 PDE给出了直接求解初边 值问题的系统指令，格式如下： solve = pdepem, fun, icfun, bcfun, xmesh, tspan, options, p, p2, ) (4) 其中 S 2)是将所求解的 PDE转化成形如 ( c x9tu, V dxj du

9、 =x dt、 dx Xm,f du X, ty W. dx w ) du l dx) (5) 时的数值， /c/iw和 分 别 是 对 初 始 条 件 和 边 界 条 件 计 算 的 函 数 。 xwes/7是需要求解的对应空间方向的网格点，其余 参数含义同上。 对于其他形式的 PDE,在进行数值求解的时候，可以根据问题所处的实际物理背景，分析函数特征，然后采用优化算法 , 使用 MATLAB中提供的相应优化命令来编程求解 4 实例 下面通过一个具体的实例来进一步形象地说明利用 MATLAB求解微分方程的问题。 对常微分方程 w” + c|w| = 0,求满足边界条件 M(0)=0,W(4)

10、 = _2的解。 为了求解上述方程，我们需要将已知方程转化为类似 =/ cj) 的 形 式 。 令 4 = w, _ y2 = z/ 则 有 ： _yl = 少 2,少 2 =z; =-c|z/|。 边界分别为 ： x = 0和 ;r = 4,边界条件对应的函数为：办 c(ya(l)， )(l)+2)=。 这样，通过编写下列 M文件就可以求解该方程了 fun. m以向量的形式书写需要求解的微分方程 function dydx=fun (x, y, c) dydx=y(2) -c*abs(y(l) ： be. m保存计算对应的边界条件 ， 函数输出是对应的满足边界时的残差向量 function

11、rcsi=bcfun (ya, yb, c) 414 科技信息 计算机与网络 对应的点就是所求得的数值点，对应的曲线是系统经过插值以后得到的近似曲线。 5结论 本文介绍了微分方程数值解法的理论及其在 MATLAB中的实现。对相关的函数调用和命令要求进行了详细的阐述，并给出具 体的求解实例，对 MATLAB软件的初接触者来讲具有很大的实用性和借鉴性。 参考文献 1 Michael T. Heath, Scientific Computing: an introductory survey 北京：淸华大学出版社， 2001 2 张志涌等编著，精通 MATLAB 6. 5版北京：北京航空航天大学出版

12、社， 2005 3 钱祥征著，常微分方程解题方法长沙：湖南铒科学技术出版社， 1987 resi = ya(l) yb(l)+2j,; program 调用指令求解给定边界条件的微分方程，这部分代码也可以写在命令行屮直接运行。 tmesh=bvpinit (linspace (0, 4, 5),1;0); c二 1: solve二 bvp4c(fun, bcfun, tmesh, , c); x=0:0.05:4; y二 bvpval (solve, x); plot(x， y(l, :)， k- ） hold on; plot (solve, x, solve, y (1, :) / ) xlabel ( x ）； ylabclC solution of ordinary differential equation*) legendC curve after interpolation ， solution point*,1) hold off; 将上述三个 M文件保存在相同的路径下面，运行 program, m后即可得到下面的图形： ualGnbsl-H-Bjjip-elllpjojcluollnlos 415