《用matlab求解常微分方程--修订编选》由会员分享，可在线阅读，更多相关《用matlab求解常微分方程--修订编选（7页珍藏版）》请在金锄头文库上搜索。
1、实验六用 matlab 求解常微分方程 1 微分方程的概念 未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为 微分方程。如果未知函数是一元函数，称为常微分方程。常微分方程的一般形式为 0), ,( )(n yyyytF 如果未知函数是多元函数，成为偏微分方程。 联系一些未知函数的一组微分方程组称为 微分方程组。 微分方程中出现的未知函数的导数的最高阶解数称为微分方程的阶。若方程中 未知函数及其各阶导数都是一次的，称为线性常微分方程，一般表示为 )()()()( 1 )1( 1 )( tbytaytaytay nn nn 若上式中的系数 nitai,2, 1),( 均与
2、 t无关，称之为常系数。 2 常微分方程的解析解 有些微分方程可直接通过积分求解. 例如 , 一解常系数常微分方程 1y dt dy 可化为 dt y dy 1 , 两边积分可得通解为 1 t cey . 其中 c为任意常数 . 有些常微分方程可用一些 技巧 , 如分离变量法, 积分因子法 , 常数变异法 , 降阶法等可化为可积分的方程而求得解析解. 线性常微分方程的解满足叠加原理, 从而他们的求解可归结为求一个特解和相应齐次微 分方程的通解 . 一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分 方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。 一阶常微分
3、方程与高阶微分方程可以互化，已给一个 n 阶方程 ), ,( )1()(nn yyytfy 设 )1( 21, , n nyyyyyy ，可将上式化为一阶方程组 ),( 21 1 32 21 nn nn yyytfy yy yy yy 反过来，在许多情况下，一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高 阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根 法求解。 3微分方程的数值解法 除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大 部分微分方程无限世界，应用中主要依靠数值解法。考虑一阶常微分方程初值问题 00 0 )( ),(
4、,()( yty ttttytfty f 其中 ).,(,),(,),(0201002121mmmyyyyffffyyyy 所谓数值解法，就 是 寻 求 )(ty 在 一 系 列 离 散 节 点 fn tttt 10 上 的 近 似 值 nkyk, 1 ,0, 称 kkk tth 1 为步长，通常取为常量 h。最简单的数值解法是 Euler 法。 Euler 法的思路极其简单：在节点出用差商近似代替导数 h tyty ty kk k )()( )( 1 这样导出计算公式（称为Euler 格式） ,2, 1 , 0),( 1 kythfyy kkkk 他能求解各种形式的微分方程。Euler法也称
5、折线法。 Euler 方法只有一阶精度，改进方法有二阶Runge-Kutta法、四阶 Runge-Kutta法、五 阶 Runge-Kutta-Felhberg法和先行多步法等，这些方法可用于解高阶常微分方程（组）初 值问题。边值问题采用不同方法，如差分法、有限元法等。数值算法的主要缺点是它缺乏物 理理解。 4 解微分方程的MATLAB 命令 MATLAB 中主要用dsolve 求符号解析解，ode45,ode23,ode15s求数值解。 s=dsolve( 方程1, 方程2, , 初始条件1，初始条件 2, 自变量) 用字符串方程表示，自变量缺省值为t 。导数用D 表示， 2 阶导数用 D2
6、表示，以此类推。S返回解析解。在方程组情形，s 为一个符号结构。 tout,yout=ode45(yprime ,t0,tf,y0) 采 用 变 步 长 四 阶 Runge-Kutta法和五阶Runge-Kutta-Felhberg法求数值解， yprime 是用以 表示 f(t,y)的 M文件名， t0 表示自变量的初始值，tf表示自变量的终值， y0 表示初始向量值。 输出向量 tout表示节点 (t0,t1, ,tn) T,输出矩阵 yout 表示数值解，每一列对应y 的一个分量。若无输出参数，则自动作出图形。 ode45 是最常用的求解微分方程数值解的命令，对于刚性方程组不宜采用。od
7、e23 与 ode45 类似，只是精度低一些。 ode12s 用来求解刚性方程组，是用格式同ode45。 可以用 help dsolve, help ode45查阅有关这些命令的详细信息. 例 1求下列微分方程的解析解 （1） bayy （2） 1)0( ,0)0(,)2sin( yyyxy （3） 1)0( , 1)0( ,gffgggff 方程（ 1）求解的MATLAB 代码为： clear; s=dsolve(Dy=a*y+b) 结果为 s =-b/a+exp(a*t)*C1 方程（ 2）求解的MATLAB 代码为： clear; s=dsolve(D2y=sin(2*x)-y,y(0)
8、=0,Dy(0)=1,x) simplify(s) %以最简形式显示s 结果为 s =(-1/6*cos(3*x)-1/2*cos(x)*sin(x)+(-1/2*sin(x)+1/6*sin(3*x)*cos(x)+5/3*sin(x) ans =-2/3*sin(x)*cos(x)+5/3*sin(x) 方程（ 3）求解的MATLAB 代码为： clear; s=dsolve(Df=f+g,Dg=g-f,f(0)=1,g(0)=1) simplify %s 是一个结构 simplify 结果为 ans =exp(t)*cos(t)+exp(t)*sin(t) ans =-exp(t)*si
9、n(t)+exp(t)*cos(t) 例求解微分方程 ,1)0(, 1ytyy 先求解析解，再求数值解，并进行比较。由 clear; s=dsolve(Dy=-y+t+1,y(0)=1,t) simplify(s) 可得解析解为 t ety 。下面再求其数值解，先编写M文件 %M 函数 function f=fun8(t,y) f=-y+t+1; 再用命令 clear; close; t=0:1; y=t+exp(-t); plot(t,y); %化解析解的图形 hold on; %保留已经画好的图形，如果下面再画图，两个图形和并在一起 t,y=ode45(fun8,0,1,1); plot(
10、t,y,ro); %画数值解图形，用红色小圈画 xlabel(t),ylabel(y) 结果见图 00.20.40.60.81 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 t y 图 16.6.1 解析解与数值解 由图 16.6.1可见，解析解和数值解吻合得很好。 例 3 求方程 0)0( ,)0(,sin 0 mgml 的数值解 .不妨取 15)0(, 8.9, 1 gl . 则上面方程可化为 0)0( ,15)0(,sin8.9 先看看有没有解析解. 运行 MATLAB 代码 clear; s=dsolve(D2y=*sin(y),y(0)=15,Dy(0)
11、=0,t) simplify(s) 知原方程没有解析解. 下面求数值解. 令 , 21 yy 可将原方程化为如下方程组 0)0(,15)0( )sin(8. 9 21 12 21 yy yy yy 建立 M文件如下 %M 文件 function f=fun9(t,y) f=y(2), *sin(y(1); %f 向量必须为一列向量 运行 MATLAB 代码 clear; close; t,y=ode45(fun9,0,10,15,0); plot(t,y(:,1); %画随时间变化图,y(:2)则表示 的值 xlabel(t),ylabel(y1) 结果见图 012345678910 15 15.5 16 16.5 t y 1 图 数值解 由图可见，随时间t周期变化。 习题 16-6 1 求下列微分方程的解析解 2 求方程 3)0( , 1)0(, 2)1 ( 2 yyxyyx 的解析解和数值解,并进行比较 3分别用ode45 和 ode15s 求解 Van-del-Pol方程 1)0, 0)0( 0)1(1000 2 2 2 xx x dt dx x dt xd 的数值解 ,并进行比较 .
