(完整版)实验七用matlab求解常微分方程

1、实验七用matlab求解常微分方程一、实验目的：1、熟悉常微分方程的求解方法，了解状态方程的概念；2、能熟练使用dsolve函数求常微分方程(组)的解析解；3、能熟练应用ode45ode15s函数分别求常微分方程的非刚性、刚性的数值解；4、掌握绘制相图的方法二、预备知识：1 .微分方程的概念未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为 微分方程。如果未知函数是一元函数，称为常微分方程。常微分方程的一般形式为F(t,y,y,y, ,y(n) 0如果未知函数是多元函数，成为偏微分方程。联系一些未知函数的一组微分方程组称为 微分方程组。微分方程中出现的未知函数的导数的最

2、高阶解数称为微分方程的阶。若方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，一般表示为y ai(t)y(n1) an i(t)y an(t)y b(t)若上式中的系数ai (t)，i1，2，n均与t无关，称之为常系数。2 .常微分方程的解析解dy .y 1有些微分方程可直接通过积分求解.例如，一解常系数常彳分方程 出 可化为dydtty 1，两边积分可得通解为 y ce 1.其中c为任意常数.有些常微分方程可用一些技巧,如分离变量法，积分因子法，常数变异法，降阶法等可化为可积分的方程而求得解析解.线性常微分方程的解满足叠加原理，从而他们的求解可归结为求一个特解和相应齐次微 分方程的通

3、解.一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分 方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。一阶常微分方程与高阶微分方程可以互化，已给一个 n阶方程(n)(n 1)、y f(t, y ,y , , y )(n 1)设y1 y，y2 y，yny ,可将上式化为一阶方程组yy2y2, y3 yn 1 Vnyn f (t,y1, y2, yn)反过来，在许多情况下，一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高 阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根法求解。3 .微分方程的数值解法除常系数线性微分方程可

4、用特征根法求解，少数特殊方程可用初等积分法求解外，大 部分微分方程无限世界，应用中主要依靠数值解法。考虑一阶常微分方程初值问题y(t)f(t,y(t),t0 t tfy(to) v。其中 y( yl, y2,ym ) , f ( f1,f2, fm ) ,y0(y10,y20,Vm。).所谓数值解法，就是寻求y在一系列离散节点tot1 tn tf上的近似值yk，k。，1，n称hktk 1tk为步长，通常取为常量h。最简单的数值解法是Euler法。Euler法的思路极其简单：在节点出用差商近似代替导数y(tk)y(tki) y(tk)h这样导出计算公式(称为Euler格式)ykiyk hf(tk

5、,yk),k 0,1,2,他能求解各种形式的微分方程。Euler法也称折线法。Euler方法只有一阶精度，改进方法有二阶Runge-Kutta法、四阶Runge-Kutta法、五阶Runge-Kutta-Felhberg法和先行多步法等，这些方法可用于解高阶常微分方程(组)初值问题。边值问题采用不同方法，如差分法、 有限元法等。数值算法的主要缺点是它缺乏物理理 解。4 .解微分方程的 MATLAB命令MATLAB中主要用dsoke求符号解析解，ode45,ode23,ode15s求数值解。s=dsolve(方程1方程2工，初始条件1,初始条件2，自变量)用字符串方程表示，自变量缺省值为t。导数

6、用D表示，2阶导数用D2表示，以此类推。S返回解析解。在方程组情形， s为一个符号结构。tout,yout=ode45( yprime ,t0tffyOji 步长四阶 Runge-Kutta 法和五阶 Runge-Kutta-Felhberg法求数值解，yprime是用以表示 f(t,y)的M文件名，t0表示自变 量的初始值，tf表示自变量的终值， y0表示初始向量值。输出向量tout表示节点(to,t1,tn)T,输出矩阵yout表示数值解，每一列对应 y的一个分量。若无输出参数，则 自动作出图形。ode45是最常用的求解微分方程数值解的命令，对于刚性方程组不宜采用。 ode23与ode45

7、类似，只是精度低一些。ode12s用来求解刚性方程组，是用格式同ode45。可以用helpdsolve, help ode45查阅有关这些命令的详细信息例1求下列微分方程的解析解(1) y ay b y sin(2x) y,y(0) 0,y(0) 1 f f g,g g f,f(0) 1,g(0) 1方程(1)求解的MATLAB 代码为：clear;s=dsolve(Dy=a*y+b)结果为5 =-b/a+exp(a*t)*C1方程(2)求解的MATLAB 代码为：clear;s=dsolve(D2y=sin(2*x)-y,y(0)=0,Dy(0)=1,x)simplify(s) %以最简形式

8、显示 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.f) %s是一一个结构simplify(s.g)结果为ans =exp(t)*cos(t)+exp(t)*sin(t)ans =-exp(t)*sin(t)+exp(t)*cos(t)例2求解微分方程y y t 1

9、, y(0) 1,先求解析解，再求数值解，并进行比较。由clear;s=dsolve(Dy=-y+t+1,y(0)=1,t)simplify(s)可得解析解为y t e t O下面再求其数值解，先编写 M文件fun8.m%MK数 fun8.mfunction f=fun8(t,y)f=-y+t+1;再用命令clear; close; t=0:0.1:1;y=t+exp(-t); plot(t,y); %化解析解的图形hold on; %保留已经画好的图形，如果下面再画图，两个图形和并在一起t,y=ode45(fun8,0,1,1);plot(t,y,ro); %画数值解图形，用红色小圈画xla

10、bel(t),ylabel(y)结果见图7.11.4_tIIk,1.35.1.3.1.25.y 1.2.r1.15 一/-1.1-.1.05-.1 1c00.20.40.60.81t图16.6.1解析解与数值解由图16.6.1可见，解析解和数值解吻合得很好。例3求方程ml mg sin , (0)0, (0) 0的数值解.不妨取l1,g 9.8, (0) 15 .则上面方程可化为9.8sin , (0) 15, (0) 0先看看有没有解析解.运行MATLAB代码clear;s=dsolve(D2y=9.8*sin(y),y(0)=15,Dy(0)=0,t)simplify(s)知原方程没有解析

11、解.下面求数值解.令y1,y2可将原方程化为如下方程组yV2V29.8sin(y1)y(0)15(0)0建立 M文彳fun9.m如下%MC件 fun9.mfunction f=fun9(t,y)f=y(2), 9.8*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)结果见图7.20123456t78910图7.2数值解由图7.2可见， 随时间t周期变化。以上这些都是常微分方程的精确解法，也称为常

12、微分方程的符号解。但是，我们知道，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面， MATLAB具有丰富的函数，我们将其统称为 solver,其一般格式为：T,Y=solver（odefun,tspan,y0）该函数表示在区间tspan=t0,tf上，用初始条件y0求解显式常微分方程y f(t,y)solver 为命令 ode45, ode23, ode113, ode15s ode23s, ode23t, ode23tb 之,这些命令各有特点。我们列表说明如下：求解 器ODE类型特点说明ode45非刚性算法，

13、4,5 阶 Runge-Kutta方法累积截断误差（x）大部分场合的首选算 法ode23非刚性算法，2,3 阶 Runge-Kutta方法累积截断误差（x）使用于精度较低的情 形ode113非刚性多步法，Adams算法，高低精度均可达到10 3 10 6计算时间比ode45短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法，Gear反向 数值积分，精度中等若ode45失效时， 可尝试使用ode23s刚性2 阶 Rosebrock算法，低精度。当精度较低时，计算时间比ode15s短odefun为显式常微分方程 寸f&y)中的f (t, y)tspan为求解区间，要获得问题在其他

14、指定点 3t1,t2,L上的解，则令tspan t0，t1，t2，L，tf(要求 t 单调)，y0初始条件。2例5:求解常微分方程y2y 2x 2x , 0 x 0.5, y 1的MATLAB程序如下：fun=inline(-2*y+2*x*x+2*x) ； x,y=ode23(fun,0,0.5,1)结果为：x =0,0.0400,0.0900,0.1400,0.1900,0.2400,0.2900,0.3400,0.3900,0.4400,0.4900,0.5000y =1.0000,0.9247,0.8434,0.7754,0.7199,0.6764,0.6440,0.6222,0.61

15、05,0.6084,0.6154,0.6179例6:求解常微分方程d2y dt2(1 y* y 0,y(0) 1,y(0)0,、，一 的解，并回出解的图形。分析：这是一个二阶非线性方程，用现成的方法均不能求解，但我们可以通 过下面的变换，将二阶方程化为一阶方程组，即可求解。dyx2-7一.令：x1y, dt ,7,则得到:dx1x2, Xi(0)1dt2 7(1 x12)x2 x1, x2(0)0dt接着，编写vdp.m如下：function fy=vdp(t,x)fy=x(2);7*(1-x(1)A2)*x(2)-x(1);再编写m文件sy12_6.m如下：y0=1;0t,x=ode45(vdp,0,40,y0);y=x(:,1);dy=x(:,2);7 / 7plot(t,y,t,dy)三、实验内容:1 .求下列微分方程的解析解(1) 了rr +2了，-3尸=6一了一3y = 2/皿乂13) 了“十一尸=血工(口 0) y/-/2 -1 = 0了啾