实验七用matlab求解常微分方程
matlab-常微分方程
Events
含义 为‘on’时，控制解向量 有效值： 范数的相对误差，使每 on、off 步计算中，满足： 缺省值： norm(e)<=max(RelTol*n off orm(y),AbsTol) 有效值： 为‘on’时，返回相应的 on、off 事件记录
取值
参数设置
属性名 含义 若无输出参量，则solver 将执行下面操作之一： 有效值： 画出解向量中各元素随 odeplot、 时间的变化； odephas2、画出解向量中前两个分 odephas3、量构成的相平面图； odeprint 画出解向量中前三个分 缺省值： 量构成的三维相空间图 odeplot ； 随计算过程，显示解向 量 取值
使用于精度较低 的情形
OD 求解器 E类 Solver 型 非 ode113 刚 性 适 度 ode23t 刚 性 刚 ode15s 性
特点
说明
多步法；Adams算 计算时间比ode45 法；高低精度均可 短 -3～10-6 到10 采用梯形算法 适度刚性情形
多步法；Gear’s反 若ode45失效时， 向数值微分；精度 可尝试使用 中等
• (6)若没有给定输出参量，则在命令窗口显 示解列表。若该命令找不到解析解，则返 回一警告信息，同时返回一空的sym对象。 这时，用户可以用命令ode23或ode45求解 方程组的数值解。
y ′′ = −a y ′ y ( 0) = 1 y ′(π / a ) = 0
2
例1
例2
• >> [u,v] = dsolve('Du=v,Dv=u') u= C1*exp(-t)+C2*exp(t) V= -C1*exp(-t)+C2*exp(t)
matlab常微分方程和常微分方程组的求解
下载的，感觉不错，共享一下常微分方程和常微分方程组的求解一、实验目的：熟悉Matlab 软件中关于求解常微分方程和常微分方程组的各种命令，掌握利用Matlab 软件进行常微分方程和常微分方程组的求解。
二、相关知识在MATLAB 中，由函数dsolve()解决常微分方程（组）的求解问题，其具体格式如下：X=dsolve(‘eqn1’,’eqn2’,…)函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。
例1：求解常微分方程1dy dx x y =+的MATLAB 程序为：dsolve('Dy=1/(x+y)','x')，注意，系统缺省的自变量为t ，因此这里要把自变量写明。
结果为：-lambertw(-C1*exp(-x-1))-x-1其中：Y=lambertw(X)表示函数关系Y*exp(Y)=X 。
例2：求解常微分方程2'''0yy y -=的MATLAB 程序为：Y2=dsolve('y*D2y-Dy^2=0’,’x’) 结果为：Y2 =[ exp((x+C2)/C1)][ C2]我们看到有两个解，其中一个是常数。
例3：求常微分方程组253tt dx x y e dt dy x y e dt ⎧++=⎪⎪⎨⎪--=⎪⎩通解的MATLAB 程序为：[X,Y]=dsolve('Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t)','t')例4：求常微分方程组020210cos ,224,0t t t dx dy x t x dt dt dx dy y e y dt dt =-=⎧+-==⎪⎪⎨⎪++==⎪⎩通解的MATLAB 程序为：[X,Y]=dsolve('Dx+2*x-Dy=10*cos(t),Dx+Dy+2*y=4*exp(-2*t)','x(0)=2','y(0)=0')以上这些都是常微分方程的精确解法，也称为常微分方程的符号解。
matlab常微分方程的数值解法实验报告
实验四常微分方程的数值解法 指令：[t,y]=ode23(‘fun ’,tspan,yo) 2/3阶龙格库塔方法 [t,y]=ode45(‘fun ’,tspan,yo) 4/5阶龙格库塔方法 [t,y]=ode113(‘fun ’,tspan,yo) 高阶微分方程数值方法其中fun 是定义函数的文件名。
该函数fun 必须以为dx 输出量，以t,y 为输入量。
tspan=[t0 tfina]表示积分的起始值和终止值。
yo 是初始状态列向量。
考虑到初始条件有00d , (0)0,d d , (0)0.d SSI S S tI SI I I I tββμ⎧=-=>⎪⎪⎨⎪=-=≥⎪⎩ (5.24) 这就是Kermack 与McKendrick 的SIR 仓室模型. 方程(5.24)无法求出()S t 和()I t 的解析解.我们先做数值计算。
Matlab 代码为：function dy=rigid(t,y) dy=zeros(2,1); a=1; b=0.3;dy(1)=a*y(1).*y(2)-b*y(1); dy(2)=-a*y(1).*y(2);ts=0:.5:50; x0=[0.02,0.98];[T,Y]=ode45('rigid',ts,x0); %plot(T,Y(:,1),'-',T,Y(:,2),'*') plot(Y(:,2),Y(:,1),'b--') xlabel('s') ylabel('i')任务：1 画出i （t ），2分析各参数的影响例57：求解两点边值问题：0)5(,0)1(,32==='-''y y x y y x 。
（注意：相应的数值解法比较复杂）。
y=dsolve('x*D2y-3*Dy=x^2','y(1)=0,y(5)=0','x') ↙ y =-1/3*x^3+125/468+31/468*x^4例：用数值积分的方法求解下列微分方程 π21''2t y y -=+设初始时间t0=0；终止时间tf=3*pi ；初始条件0|',0|00====x x y y 。
MATLAB常微分方程的数值解法
MATLAB常微分⽅程的数值解法MATLAB常微分⽅程的数值解法⼀、实验⽬的科学技术中常常要求解常微分⽅程的定解问题，所谓数值解法就是求未知函数在⼀系列离散点处的近似值。
⼆、实验原理三、实验程序1. 尤拉公式程序四、实验内容选⼀可求解的常微分⽅程的定解问题，分别⽤以上1, 4两种⽅法求出未知函数在节点处的近似值，并对所求结果与分析解的（数值或图形）结果进⾏⽐较。
五、解答1. 程序求解初值问题取n=10源程序：euler23.m:function [A1,A2,B1,B2,C1,C2]=euler23(a,b,n,y0)%欧拉法解⼀阶常微分⽅程%初始条件y0h = (b-a)/n; %步长h%区域的左边界a%区域的右边界bx = a:h:b;m=length(x);%前向欧拉法y = y0;for i=2:my(i)=y(i-1)+h*oula(x(i-1),y(i-1));A1(i)=x(i);A2(i)=y(i);endplot(x,y,'r-');hold on;%改进欧拉法y = y0;for i=2:my(i)=y(i-1)+h/2*( oula(x(i-1),y(i-1))+oula(x(i),y(i-1))+h*(oula(x(i-1),x(i-1))));B1(i)=x(i);B2(i)=y(i);endplot(x,y,'m-');hold on;%欧拉两步公式y=y0;y(2)=y(1)+h*oula(x(1),y(1));for i=2:m-1y(i+1)=y(i-1)+2*h*oula(x(i),y(i));C1(i)=x(i);C2(i)=y(i);endplot(x,y,'b-');hold on;%精确解⽤作图xx = x;f = dsolve('Dy=-3*y+8*x-7','y(0)=1','x');%求出解析解y = subs(f,xx); %将xx代⼊解析解，得到解析解对应的数值plot(xx,y,'k--');legend('前向欧拉法','改进欧拉法','欧拉两步法','解析解');oula.m:function f=oula(x,y)f=-3*y+8*x-7;2. 运算结果A1,A2为前向欧拉法在节点处的近似值，B1,B2为改进的欧拉法在节点处的近似值，C1,C2为欧拉公式法在节点处的近似值。
matlab_常微分方程数值解法
dt 2
简朴问题可以求得解析解，多数实际问题靠数值求解 。
第4页
一阶常微分方程(ODE )初值问题 ： ODE ：Ordinary Differential Equation
dy
f
(x,
y)
dx
x0 x xn
y(x0 ) y0
数值解法就是求y(x)在某些分立旳节点 xn 上旳近似值 yn，用以近似y(xn)
x0
y0
x1 f y(x), x dx
x0
x2 f y(x), x dx
x1
y(x1) f y(x1), x1 h
第17页
同样，在[x0,xn+1] ，积分采用矩形近似，得：
y(xn1) y0
f xn1
x0
y(x), x dx
y(xn ) f y(xn ), xn h
yn y(xn )
第5页
2、欧拉近似办法
2.1 简朴欧拉(L.Euler, 1707-1783)办法。
dy
dx
f
(y, x)
y(x0 ) y0
欧拉数值算法就是由初值通过递推求解，递推求解
就是从初值开始，后一种函数值由前一种函数值得到。核 心是构造递推公式。
y0 y1 y2 yn
第6页
i 1,2,...
第36页
没有一种算法可以有效地解决所有旳 ODE 问题，因此 MATLAB 提供了多种ODE函数。
函数 ODE类
特点
阐明
型
ode45
非刚性 单步法；4，5 阶 R-K 措施；合计 大部分场合旳首选措施
截断误差为 (△x)3
ode23
非刚性 单步法；2，3 阶 R-K 措施；合计 使用于精度较低旳情形
MATLAB求解常微分方程数值解
利用MATLAB求解常微分方程数值解目录1.容简介12.Euler Method(欧拉法)求解12.1.显式Euler法和隐式Euler法22.2.梯形公式和改进Euler法32.3.Euler法实用性53.Runge-Kutta Method(龙格库塔法)求解53.1.Runge-Kutta基本原理63.2.MATLAB中使用Runge-Kutta法的函数74.使用MATLAB求解常微分方程84.1.使用ode45函数求解非刚性常微分方程84.2.刚性常微分方程95.总结9参考文献11附录121.显式Euler法数值求解122.改进Euler法数值求解123.四阶四级Runge-Kutta法数值求解134.使用ode45求解141.容简介把《高等工程数学》看了一遍，增加对数学容的了解，对其中数值解法比较感兴趣，这大概是因为在其它各方面的学习和研究中经常会遇到数值解法的问题。
理解模型然后列出微分方程，却对着方程无从下手，无法得出精确结果实在是让人难受的一件事情。
实际问题中更多遇到的是利用数值法求解偏微分方程问题，但考虑到先从常微分方程下手更为简单有效率，所以本文只研究常微分方程的数值解法。
把一个工程实际问题弄出精确结果远比弄清楚各种细枝末节更有意思，因此文章中不追求非常严格地证明，而是偏向如何利用工具实际求解出常微分方程的数值解，力求将课程上所学的知识真正地运用到实际方程的求解中去，在以后遇到微分方程的时候能够熟练运用MATLAB得到能够在工程上运用的结果。
文中求解过程中用到MATLAB进行数值求解，主要目的是弄清楚各个函数本质上是如何对常微分方程进行求解的，对各种方法进行MATLAB编程求解，并将求得的数值解与精确解对比，其中源程序在附录中。
最后考察MATLAB中各个函数的适用围，当遇到实际工程问题时能够正确地得到问题的数值解。
2.Euler Method(欧拉法)求解Euler法求解常微分方程主要包括3种形式，即显式Euler法、隐式Euler法、梯形公式法，本节容分别介绍这3种方法的具体容，并在最后对3种方法精度进行对比，讨论Euler法的实用性。
matlab数值求解常微分方程快速方法
MATLAB是一种用于科学计算和工程应用的高级编程语言和交互式环境。
它在数学建模、模拟和分析等方面有着广泛的应用。
在MATLAB 中，常微分方程的数值求解是一个常见的应用场景。
在实际工程问题中，通常需要对常微分方程进行数值求解来模拟系统的动态行为。
本文将介绍MATLAB中对常微分方程进行数值求解的快速方法。
1. 基本概念在MATLAB中，可以使用ode45函数来对常微分方程进行数值求解。
ode45是一种常用的Runge-Kutta法，它可以自适应地选取步长，并且具有较高的数值精度。
使用ode45函数可以方便地对各种类型的常微分方程进行求解，包括一阶、高阶、常系数和变系数的微分方程。
2. 函数调用要使用ode45函数进行常微分方程的数值求解，需要按照以下格式进行函数调用：[t, y] = ode45(odefun, tspan, y0)其中，odefun表示用于描述微分方程的函数，tspan表示求解的时间跨度，y0表示初值条件，t和y分别表示求解得到的时间序列和对应的解向量。
3. 示例演示为了更好地理解如何使用ode45函数进行常微分方程的数值求解，下面我们以一个具体的例子来进行演示。
考虑如下的一阶常微分方程：dy/dt = -2*y其中，y(0) = 1。
我们可以编写一个描述微分方程的函数odefun：function dydt = odefun(t, y)dydt = -2*y;按照上述的函数调用格式，使用ode45函数进行求解：tspan = [0 10];y0 = 1;[t, y] = ode45(odefun, tspan, y0);绘制出解曲线：plot(t, y);4. 高级用法除了基本的函数调用方式外，MATLAB中还提供了更多高级的方法来对常微分方程进行数值求解。
可以通过设定选项参数来控制数值求解的精度和稳定性，并且还可以对刚性微分方程进行求解。
5. 性能优化在实际工程应用中，常常需要对大规模的常微分方程进行数值求解。
实验报告七常微分方程初值问题的数值解法
浙江大学城市学院实验报告课程名称数值计算方法实验项目名称常微分方程初值问题的数值解法 实验成绩指导老师签名日期2015/12/16 一.实验目的和要求1． 用Matlab 软件掌握求微分方程数值解的欧拉方法和龙格－库塔方法； 2． 通过实例学习用微分方程模型解决简化的实际问题;二.实验内容和原理编程题2-1要求写出Matlab 源程序m 文件,并有适当的注释语句；分析应用题2-2,2-3,2-4,2-5要求将问题的分析过程、Matlab 源程序和运行结果和结果的解释、算法的分析写在实验报告上; 2-1 编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序,问题如下：在区间[],a b 内(1)N +个等距点处,逼近下列初值问题的解,并对程序的每一句添上注释语句; Euler 法y=eulera,b,n,y0,f,f1,b1改进Euler 法y=eulerproa,b,n,y0,f,f1,b1 2-2 分析应用题假设等分区间数100n =,用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度; 2-3 分析应用题用以下三种不同的方法求下述微分方程的数值解,取10h = 画出解的图形,与精确值比较并进行分析; 1欧拉法； 2改进欧拉法； 3龙格－库塔方法；2-4 分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型;假设在时刻t 单位为年,社会上有人口()x t 人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人;而固定比例为r 的所有其他的后代也是与众不同的人;如果对所有人来说出生率假定为常数b ,又如果普通的人和与众不同的人的婚配是任意的,则此问题可以用微分方程表示为：其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例,()i x t 表示在时刻t 人口中与众不同的人的数量;1假定(0)0.01,0.02p b ==和0.1r =,当步长为1h =年时,求从0t =到50t =解()p t 的近似值,并作出近似解的曲线图形;2精确求出微分方程的解()p t ,并将你当50t =时在分题b 中得到的结果与此时的精确值进行比较; MATLAB 相关函数求微分方程的解析解及其数值的代入dsolve‘egn1’,‘egn2’,‘x ’ subsexpr,{x,y,…},{x1,y1,…}其中‘egn i ’表示第i 个方程,‘x ’表示微分方程中的自变量,默认时自变量为t ; subs 命令中的expr 、x 、y 为符合型表达式,x 、y 分别用数值x1、x2代入; >>symsxyz>>subs'x+y+z',{x,y,z},{1,2,3} ans= 6>>symsx>>subs'x^2',x,2 ans= 4>>s=dsolve‘12Dy y ∧=+’,‘(0)1y =’,‘x ’ ans= >>symsx >>subss,x,2 ans=右端函数(,)f x y 的自动生成f=inline ‘expr ’,’var1’,‘var2’,……其中’expr ’表示函数的表达式,’var1’,‘var2’表示函数表达式中的变量,运行该函数,生成一个新的函数表达式为fvar1,var2,……; >>f=inline'x+3y','x','y' f=Inlinefunction: fx,y=x+3y >>f2,3 ans= 114,5阶龙格－库塔方法求解微分方程数值解t,x=ode45f,ts,x0,options其中f 是由待解方程写成的m 文件名；x0为函数的初值；t,x 分别为输出的自变量和函数值列向量,t的步长是程序根据误差限自动选定的;若ts=t0,t1,t2,…,tf,则输出在自变量指定值,等步长时用ts=t0:k:tf,输出在等分点；options 用于设定误差限可以缺省,缺省时设定为相对误差310-,绝对误差610-,程序为：options=odeset ‘reltol ’,rt,’abstol ’,at,这里rt,at 分别为设定的相对误差和绝对误差;常用选项见下表;选项名 功能 可选值 省缺值 AbsTol 设定绝对误差正数 RelTol 设定相对误差 正数InitialStep 设定初始步长 正数 自动 MaxStep设定步长上界正数MaxOrder 设定ode15s 的最高阶数 1,2,3,4,5 5 Stats 显示计算成本统计 on,off off BDF 设定ode15s 是否用反向差分on,offoff例：在命令窗口执行>>odefun =inline ‘2*y t y -’,‘t ’,‘y ’；>>[],45(,[0,4],1)t y ode odefun =；ans=>>t y ‘o-’,%解函数图形表示>>45(,[0,4],1)ode odefun %不用输出变量,则直接输出图形 >>[],45(,0:4,1)t y ode odefun =；[],t yans=三.操作方法与实验步骤包括实验数据记录和处理2-1编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序,问题如下：在区间[],a b 内(1)N +个等距点处,逼近下列初值问题的解,并对程序的每一句添上注释语句; Euler 法y=eulera,b,n,y0,f,f1,b1改进Euler 法y=eulerproa,b,n,y0,f,f1,b1Euler 法y=eulera,b,n,y0,f,f1,b1 y=zeros1,n+1; y1=y0; h=b-a/n; x=a:h:b; fori=1:n; yi+1=yi+hfxi,yi; end plotx,y holdon%求微分方程的精确解 x1=linspacea,b,100; '精确解为' s=dsolvef1,b1,'x' symsxy1=zeros1,100; for i=1:100y1i=subss,x,x1i; endplotx1,y1,'r'title'红色代表精确解'改进Euler 法y=eulerproa,b,n,y0,f,f1,b1 %求微分方程的数值解 y=zeros1,n+1; y1=y0; h=b-a/n; x=a:h:b; fori=1:n; T1=fxi,yi; T2=fxi+1,yi+hT1; yi+1=yi+h/2T1+T2; end plotx,y holdon%求微分方程的精确解 x1=linspacea,b,100; '精确解为' s=dsolvef1,b1,'x' symsxy1=zeros1,100; fori=1:100 y1i=subss,x,x1i; endplotx1,y1,'r'title'红色代表精确解' 2-2分析应用题假设等分区间数100n =,用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度;1向前欧拉法>>euler0,10,100,10,inline'y-20','x','y','Dy=y-20','y0=10' ans= 精确解为 s= 20-10expx ans= +005Columns1through8（2）改进欧拉法>>eulerpro0,10,100,10,inline'y-20','x','y','Dy=y-20','y0=10' ans= 精确解为 s= 20-10expx ans= +005Columns1through8改进欧拉法的精度比向前欧拉法更高; 2-3分析应用题用以下三种不同的方法求下述微分方程的数值解,取10h = 画出解的图形,与精确值比较并进行分析; 1欧拉法； 2改进欧拉法；2-4分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型;假设在时刻t 单位为年,社会上有人口()x t 人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人;而固定比例为r 的所有其他的后代也是与众不同的人;如果对所有人来说出生率假定为常数b ,又如果普通的人和与众不同的人的婚配是任意的,则此问题可以用微分方程表示为：其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例,()i x t 表示在时刻t 人口中与众不同的人的数量;1假定(0)0.01,0.02p b ==和0.1r =,当步长为1h =年时,求从0t =到50t =解()p t 的近似值,并作出近似解的曲线图形;2精确求出微分方程的解()p t ,并将你当50t =时在分题b 中得到的结果与此时的精确值进行比较;1>>euler0,50,50,,inline'','t','p','Dp=','p0= 1' ans= 精确解为 s=1-99/100expx/500 ans=Columns1through82>>dsolve'Dp=','p0=','t' ans=1-99/100expt/500 >>1-99/100exp ans=与欧拉法求得的精确值差0,0001四.实验结果与分析。
用MATLAB求解微分方程
1. 微分方程的解析解
求微分方程（组）的解析解命令:
dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’)
结 果：u = tan(t-c)
解 输入命令：dsolve('Du=1+u^2','t')
STEP2
STEP1
解 输入命令: y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')
导弹追踪问题
设位于坐标原点的甲舰向位于x轴上点A(1, 0)处的乙舰发射导弹，导弹头始终对准乙舰.如果乙舰以最大的速度v0(是常数)沿平行于y轴的直线行驶，导弹的速度是5v0，求导弹运行的曲线方程.又乙舰行驶多远时，导弹将它击中？ 解法一（解析法）
由(1),(2)消去t整理得模型:
解法二(数值解)
结 果 为：x = (c1-c2+c3+c2e -3t-c3e-3t)e2t y = -c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t z = (-c1e-4t+c2e-4t+c1-c2+c3)e2t
2、取t0=0，tf=12，输入命令： [T,Y]=ode45('rigid',[0 12],[0 1 1]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')
3、结果如图
图中，y1的图形为实线，y2的图形为“*”线，y3的图形为“+”线.
用MATLAB计算常微分方程
dsolve('Dy=f(x,y)','y(0) =a','x')
dsolve('D2y=f(x,y,Dy)',' y(0)=a','Dy(0)=b','x')
二、例题解析
例5、求以下微分方程和微分方程组： (1)求 y x 的通解； (2)求 y y 1 的通解； (3)求 y y 1 满足初始条件y(0)=1，Dy(0)=0
的特解。
解 输入及结果如下：
(1) >>syms x y >>y=dsolve('Dy=x','x') y =1/2*x^2+C1 (2) >>dsolve (‘D2y= Dy+1’,’x’) ans =-x+C1+C2*exp(x) (3) >>dsolve('D2y= Dy+1','y(0)=1','Dy(0)=0','x') ans =-x+exp(x)
输 入 格 式 dsolve('Dy=f(x,y)','x') 含 义
求一阶微分方程 y ' f ( x, y ) 的通解 求一阶微分方程 y' f ( x, y), y(0) a 的特解 求二阶微分方程 y ' ' f ( x, y, y ' ), y (0) a, y ' (0) b 的特解
dy y y tan 的通解 例2 求微分方程 dx x x 解：>> y=dsolve(‘Dy=y/x+tan(y/x)’,’x’) y= asin(x*C1)*x dy 2 xy 2 的通解和当 y (0) 3 的特解 例3 求 dx 解：>> y=dsolve(‘Dy=2*x*y^2’,’x’) y= -1/(x^2-C1)
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验七 用matlab 求解常微分方程一、实验目的：1、熟悉常微分方程的求解方法，了解状态方程的概念；2、能熟练使用dsolve 函数求常微分方程（组）的解析解；3、能熟练应用ode45\ode15s 函数分别求常微分方程的非刚性、刚性的数值解；4、掌握绘制相图的方法 二、预备知识： 1．微分方程的概念未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。
如果未知函数是一元函数，称为常微分方程。
常微分方程的一般形式为0),,",',,()(=n y y y y t F如果未知函数是多元函数，成为偏微分方程。
联系一些未知函数的一组微分方程组称为微分方程组。
微分方程中出现的未知函数的导数的最高阶解数称为微分方程的阶。
若方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，一般表示为)()(')()(1)1(1)(t b y t a y t a y t a y n n n n =++++--若上式中的系数n i t a i ,,2,1),( =均与t 无关，称之为常系数。
2．常微分方程的解析解有些微分方程可直接通过积分求解.例如,一解常系数常微分方程1+=y dt dy可化为dt y dy=+1,两边积分可得通解为1-=tce y .其中c 为任意常数.有些常微分方程可用一些技巧,如分离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解.线性常微分方程的解满足叠加原理,从而他们的求解可归结为求一个特解和相应齐次微分方程的通解.一阶变系数线性微分方程总可用这一思路求得显式解。
高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。
一阶常微分方程与高阶微分方程可以互化，已给一个n 阶方程),,",',()1()(-=n n y y y t f y设)1(21,,',-===n n y y y y y y ，可将上式化为一阶方程组⎪⎪⎪⎩⎪⎪⎪⎨⎧====-),,,,(''''2113221n n nn y y y t f y yy y y y y反过来，在许多情况下，一阶微分方程组也可化为高阶方程。
所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根法求解。
3．微分方程的数值解法除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程无限世界，应用中主要依靠数值解法。
考虑一阶常微分方程初值问题⎩⎨⎧=<<=000)()),(,()('y t y t t t t y t f t y f其中)'.,,,(,)',,,(,)',,,(020*******m m m y y y y f f f f y y y y ===所谓数值解法，就是寻求)(t y 在一系列离散节点f n t t t t ≤<<< 10上的近似值n k y k ,,1,0, =称k k k t t h -=+1为步长，通常取为常量h 。
最简单的数值解法是Euler 法。
Euler 法的思路极其简单：在节点出用差商近似代替导数h t y t y t y k k k )()()('1-≈+这样导出计算公式（称为Euler 格式）,2,1,0),,(1=+=+k y t hf y y k k k k他能求解各种形式的微分方程。
Euler 法也称折线法。
Euler 方法只有一阶精度，改进方法有二阶Runge-Kutta 法、四阶Runge-Kutta 法、五阶Runge-Kutta-Felhberg 法和先行多步法等，这些方法可用于解高阶常微分方程（组）初值问题。
边值问题采用不同方法，如差分法、有限元法等。
数值算法的主要缺点是它缺乏物理理解。
4．解微分方程的MATLAB 命令MATLAB 中主要用dsolve 求符号解析解，ode45,ode23,ode15s 求数值解。
ode45类似，只是精度低一些。
ode12s 用来求解刚性方程组，是用格式同ode45。
可以用help dsolve, help ode45查阅有关这些命令的详细信息. 例1 求下列微分方程的解析解（1）b ay y +='（2）1)0(',0)0(,)2sin(''==-=y y y x y（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)=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.f) %s是一个结构>>simplify(s.g)结果为ans =exp(t)*cos(t)+exp(t)*sin(t)ans =-exp(t)*sin(t)+exp(t)*cos(t)例２求解微分方程,1)0(,1'=++-=ytyy先求解析解，再求数值解，并进行比较。
由>>clear;>>s=dsolve('Dy=-y+t+1','y(0)=1','t')>>simplify(s)可得解析解为tety-+=。
下面再求其数值解，先编写M文件fun8.m%M函数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'); %画数值解图形，用红色小圈画>>xlabel('t'),ylabel('y')结果见图7.1图16.6.1 解析解与数值解由图16.6.1可见，解析解和数值解吻合得很好。
例3 求方程0)0(',)0(,sin "0===θθθθθmg ml的数值解.不妨取15)0(,8.9,1===θg l .则上面方程可化为0)0(',15)0(,sin 8.9"===θθθθ先看看有没有解析解.运行MATLAB 代码>>clear;>>s=dsolve('D2y=9.8*sin(y)','y(0)=15','Dy(0)=0','t') >>simplify(s)知原方程没有解析解.下面求数值解.令',21θθ==y y 可将原方程化为如下方程组⎪⎩⎪⎨⎧====0)0(,15)0()sin(8.9''211221y y y y y y建立M 文件fun9.m 如下%M 文件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.2图7.2 数值解由图7.2可见，θ随时间t周期变化。
以上这些都是常微分方程的精确解法，也称为常微分方程的符号解。
但是，我们知道，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB具有丰富的函数，我们将其统称为solver，其一般格式为：[T,Y]=solver(odefun,tspan,y0)该函数表示在区间tspan=[t0,tf]上，用初始条件y0求解显式常微分方程y f t y='(,)solver为命令ode45，ode23，ode113，ode15s，ode23s，ode23t，ode23tb之一，这些命令各有特点。
我们列表说明如下：odefun 为显式常微分方程'(,)y f t y =中的(,)f t ytspan 为求解区间，要获得问题在其他指定点012,,,t t t 上的解，则令012[,,,,]f tspan t t t t = （要求i t 单调），y0初始条件。
例5：求解常微分方程2'222y y x x =-++，00.5x ≤≤，(0)1y =的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.5000 y =1.0000,0.9247,0.8434,0.7754,0.7199,0.6764,0.6440,0.6222,0.6105,0.6084,0.6154,0.6179例6：求解常微分方程222(1)0,(0)1,'(0)0d y dy y y y y dt dt μ--+===的解，并画出解的图形。