第五讲 Matlab求解微分方程
第五讲 Matlab求解微分方程
第五讲 Matlab求解微分方程教学目的：学会用MATLAB求简单微分方程的解析解、数值解，加深对微分方程概念和应用的理解；针对一些具体的问题，如追击问题，掌握利用软件求解微分方程的过程；了解微分方程模型解决问题思维方法及技巧；体会微分方程建摸的艺术性．教学重点：利用机理分析建模微分方程模型，掌握追击问题的建模方法，掌握利用MATLAB求解数值解．教学难点：利用机理分析建模微分方程模型，通过举例，结合图形以及与恰当的假设突破教学难点．1微分方程相关函数（命令）及简介因为没有一种算法可以有效地解决所有的ODE 问题，为此，Matlab 提供了多种求解器Solver，对于不同的ODE 问题，采用不同的Solver．阶常微分方程(组)的初值问题的解的 Matlab 的常用程序，其中：ode23 采用龙格-库塔2 阶算法，用3 阶公式作误差估计来调节步长，具有低等的精度．ode45 则采用龙格-库塔4 阶算法，用5 阶公式作误差估计来调节步长，具有中等的精度．2 求解微分方程的一些例子2.1 几个可以直接用 Matlab 求微分方程精确解的例子：例1：求解微分方程22x xe xy dxdy-=+，并加以验证． 求解本问题的Matlab 程序为：syms x y %line1 y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') %line2 diff(y,x)+2 *x*y-x*exp(-x^2) %line3 simplify(diff(y,x)+2*x*y-x*exp(-x^2)) %line4说明：(1) 行line1是用命令定义x,y 为符号变量．这里可以不写，但为确保正确性，建议写上；(2) 行line2是用命令求出的微分方程的解：1/2*exp(-x^2)*x^2+exp(-x^2)*C1(3) 行line3使用所求得的解．这里是将解代入原微分方程，结果应该为0，但这里给出：-x^3*exp(-x^2)-2*x*exp(-x^2)*C1+2*x*(1/2*exp(-x^2)*x^2+exp(-x^2)*C1)(4) 行line4 用 simplify() 函数对上式进行化简，结果为 0， 表明)(x y y =的确是微分方程的解．例2：求微分方程0'=-+x e y xy 在初始条件e y 2)1(=下的特解，并画出解函数的图形．求解本问题的 Matlab 程序为： syms x yy=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') ezplot(y)微分方程的特解为：y=1/x*exp(x)+1/x* exp (1) (Matlab 格式)，即xe e y x+=，此函数的图形如图 1：图1 y 关于x 的函数图象2.2 用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)．例３：求解微分方程初值问题⎪⎩⎪⎨⎧=++-=1)0(2222y xx y dx dy 的数值解，求解范围为区间[0, 0.5]．fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,0.5],1); x'; y'; plot(x,y,'o-') >> x' ans =0.0000 0.0400 0.0900 0.1400 0.1900 0.2400 0.2900 0.3400 0.3900 0.4400 0.4900 0.5000>> y' ans =1.0000 0.9247 0.8434 0.7754 0.7199 0.6764 0.6440 0.6222 0.6105 0.6084 0.6154 0.6179 图形结果为图 2．图2 y 关于x 的函数图像3 常微分在实际中的应用 3.1 导弹追踪问题1、符号说明，w ，乙舰的速率恒为v 0；设时刻t 乙舰的坐标为((),())X t Y t ，导弹的坐标为((),())x t y t ；当零时刻，((0),(0))(1,0)X Y =，((0),(0))(0,0)x y =．建立微分方程模型．202,01(0)0,'(0)0v d yk x dx w y y ⎧⎪⎪= =<<⎨⎪⎪==⎩由微分方程模型解得11(1)(1)11(),12(1)2(1)2(1)2(1)k k x x y x k k k k k +-+--=-+-≠+-+-++代入题设的数据1/5k =，得到导弹的运行轨迹为4655555(1)(1)81224y x x =--+-+当1=x 时245=y ，即当乙舰航行到点)245,1(处时被导弹击中. 被击中时间为:00245v v y t ==. 若v 0=1, 则在t =0.21处被击中. 利用MALAB 作图如图３．clear, x=0:0.01:1; y=-5*(1-x).^(4/5)/8+5*(1-x).^(6/5)/12+5/24; plot(x,y,'*')图３导弹运行轨迹（解析法） 图４两种方法对比的导弹运行轨迹2、数值方法求解．设导弹速率恒为w ，则得到参数方程为⎪⎪⎩⎪⎪⎨⎧--+-=--+-=)()()()()()(2222y Y y Y x X wdt dy x X y Y x X w dt dx因乙舰以速度0v 沿直线1x =运动，设01v =，5,1,w X Y t ===，因此导弹运动轨迹的参数方程为：⎪⎪⎪⎩⎪⎪⎪⎨⎧==--+-=--+-=0)0(,0)0()()()1(5)1()()1(52222y x y t y t x dtdyx y t x dt dxMATLAB 求解数值解程序如下，结果见图４ t0=0,tf=0.21;[t,y]=ode45('eq2',[t0 tf],[0 0]); X=1;Y=0:0.001:0.21;plot(X,Y,'-') plot(y(:,1),y(:,2),'*'),hold onx=0:0.01:1; y=-5*(1-x).^(4/5)/8+5*(1-x).^(6/5)/12+5/24; plot(x,y,'r')很明显数值计算的方法比较简单而且适用． 3.2 蚂蚁追击问题（1）建立平面直角坐标系．以正多边形的中心为原点, 设正多边形的一个顶点为起始点, 连接此点和原点作x 轴. 根据x 轴作出相应的y 轴; 选取足够小的t ∆进行采样．（2）在每一时刻t ，计算每只蚂蚁在下一时刻t t +∆时的坐标．不妨设甲追逐对象是乙，在时间t 时，甲的坐标为A 11(,)x y ，乙的坐标为B 22(,)x y ．甲在t t +∆时在'A 点(如图1), 其坐标为11(cos ,sin )x v t y v t θθ+∆+∆，其中2221212121cos ,sin ,()()x x y yd x x y y d dθθ--===-+-. 同理，依次计算下一只蚂蚁在t t +∆时的坐标．通过间隔t ∆进行采样，得到新一轮各个蚂蚁在一个新的正多边形位置坐标．（4）重复2）步，直到d 充分小为止．（5）连接每只蚂蚁在各时刻的位置，就得到所求的轨迹．用MALAB 求解并作图，函数zhuJi(x,y)在附录一定义，如图６ t=[1:8]; s=7*exp(t.*2*pi/length(t)*i); x=real(s); y=imag(s);zhuJi(x,y)图6当蚂蚁为7只时的图形习题1. 求微分方程0sin 2')1(2=-+-x xy y x 的通解．2. 求微分方程x e y y y x sin 5'2''=+-的通解．3. 求微分方程组'A 11,(,)A x y 22,(,)B x y dθ图1 在采样时间内，相连蚂蚁追击⎪⎪⎩⎪⎪⎨⎧=-+=++00y x dtdy y x dtdx在初始条件0|,1|00====t t y x 下的特解，并画出解函数()y f x =的图形．4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解)，求解区间为[0,2]t ∈．利用画图来比较两种求解器之间的差异．4 参考文献:[1] Mastering Matlab 6, D. Hanselman, B. Littlefield, 清华大学出版,2002[2] 赵静等编．数学建模与数学实验（第3版）．北京：高等教育出版社．2008. [3] 姜启源编. 数学模型（第二版）．北京：高等教育出版社.1993.[4] 石勇国. 蚂蚁追击问题与等角螺线. 宜宾学院学报. 2008，(6): 23-25. [5] 张伟年，杜正东，徐冰．常微分方程．北京：高等教育出版社．2006．5 附录附录一:zhuji(x,y)的M 文件 function zhuji(x,y) clf v=1; dt=0.05;x(length(x)+1)=x(1);y(length(y)+1)=y(1); plot(x,y,'*') holdfor it=1:1000for i=1:length(x)-1d=sqrt((x(i)-x(i+1))^2+(y(i)-y(i+1))^2); x(i)=x(i)+v*dt*(x(i+1)-x(i))/d; y(i)=y(i)+v*dt*(y(i+1)-y(i))/d; endplot(x,y,'.') hold onx(length(x))=x(1); y(length(y))=y(1); end。
如何使用MATLAB求解微分方程(组)ppt课件
差,输出参数,事件等，可缺省。 9
使用ODE?时如何编 写微分方程 ？
方式一：带额外参数，使用时需对参数进行赋值
function odefun(t,x,flag,R,L,C) %用flag说明R、L、C为变 量
xdot=zeros(2,1);
%表明其为列向量
xdot(1)=-R/L*x(1)-1/L*x(2)+1/L*f(t);
2
Where ?
工程控制
ODE
医学生理
航空航天
金融分析
3
Where ?
算法开发 数据分析
数值计算 MAT LAB
数据可视化
4
When ?
当对问题进行建模后，有常微分方程 需要求解时。
在生物建模中，经常需要求解常微分 方程。如药物动力学的房室模型的建模 仿真。
5
方法 ode23
ode45
数 ode113
当无法藉由微积分技巧求 得解析解时，这时便只能利 用数值分析的方式来求得其 数值解了。实际情况下，常 微分方程往往只能求解出其
数值解。
在数学中，刚性方程是指一 个微分方程，其数值分析的解 只有在时间间隔很小时才会稳 定，只要时间间隔略大，其解 就会不稳定。
目前很难去精确地去定义哪 些微分方程是刚性方程，但是 大体的想法是：这个方程的解
y(1)=x(2);
y1
y2
y(2)= -t*x(1)+exp(t)*x(2)+3*sin(2*t);
end
1000
求解程序ห้องสมุดไป่ตู้键步骤
[t,y]=ode45('odefun',[3.9 4.0],[2 8])
y
matlab_微分方程求解
解
du 的通解. = 1 + u 2 的通解 dt 输入命令： 输入命令：u=dsolve('Du=1+u^2','t')
结果：u = tan (t+c1)
例如求下例微分方程的特解。 例如求下例微分方程的特解。
dy = e x , y (0) = e dx
命令窗口中输入： 在Matlab命令窗口中输入： 命令窗口中输入 y=dsolve('Dy=exp(x)','y(0)=exp(1)','x') 输出结果为： 输出结果为： y =exp(x)-1+exp(1) 如想画出函数在自变量x在区间 如想画出函数在自变量 在区间 [-10,10]的函数图像，可在 的函数图像， 的函数图像 命令窗口中输入： 命令窗口中输入： ezplot(y,[-10,10])
例2
求微分方程的特解. 求微分方程的特解
d 2 y dy 2 + 4 + 29 y = 0 dx dx y (0) = 0, y ' (0) = 15
输入命令: 解 输入命令 y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x') 结 果 为 : y =3*exp(-2*x)*sin(5*x) 作图命令： 作图命令：ezplot(y,[1.0,4])
1、微分方程的解析解 、
求微分方程（ 求微分方程（组）的解析解命令: 的解析解命令 dsolve(‘方程 ‘方程 方程1’, 方程 方程2’,…‘方程 ‘初始条件’, ‘自变量’) 方程n’, 初始条件 初始条件’ 自变量 自变量’ 方程 方程
matlab 微分方程求解
例2
求微分方程的特解. 求微分方程的特解
d 2 y dy 2 + 4 + 29 y = 0 dx dx y (0) = 0, y ' (0) = 15
输入命令: 解 输入命令 y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x') 结 果 为 : y =3*exp(-2*x)*sin(5*x) 作图命令： 作图命令：ezplot(y,[1.0,4])
解析法假设导弹在t时刻的位置为pxtyt乙舰位于由于导弹头始终对准乙舰故此时直线pq就是导弹的轨迹曲线弧op处的切线即有又根据题意弧op的长度为aq3数学建模实例由12消去t整理得模型
《高等数学》 高等数学》
—上机教学（三） 上机教学（ 微分方程求解
上机目的
1、学会用 Matlab 求简单微分方程的解析解 、 求简单微分方程的解析解. 2、学会用 Matlab 求微分方程的数值解 、 求微分方程的数值解.
例3
求微分方程组的通解. 求微分方程组的通解 dx dt = 2 x − 3 y + 3 z dy = 4 x − 5 y + 3z dt dz = 4 x − 4 y &olve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z', 't') ; x=simple(x) % 将x化简 化简 y=simple(y) z=simple(z)
假设导弹在 t 时刻的位置为 P(x(t), y(t))，乙舰位于 Q(1, v0 t ) . y(t))
matlab解带参数的微分方程
matlab解带参数的微分方程要在MATLAB中解带参数的微分方程，你可以使用MATLAB的内置函数`dsolve`。
`dsolve`函数可以用于解析解或数值解微分方程。
首先，你需要定义微分方程，然后使用`dsolve`函数来解方程。
下面我将详细介绍一下这个过程。
首先，假设我们有一个带参数的一阶微分方程，例如：syms y(t) a.eqn = diff(y,t) == ay;这里的`y(t)`是未知函数，`a`是参数，`eqn`是微分方程。
接下来，我们可以使用`dsolve`函数来解这个微分方程。
如果我们要求解的是初值问题，可以通过指定初始条件来解微分方程。
例如，如果我们有初始条件`y(0) = 1`，我们可以这样使用`dsolve`函数：cond = y(0) == 1;ySol(t) = dsolve(eqn,cond);这将给出微分方程的解`ySol(t)`，其中包含参数`a`。
如果你想要数值解而不是解析解，你可以使用`ode45`或其他数值求解器。
例如，如果我们有一个带参数的二阶微分方程：syms y(t) a.eqn = diff(y,t,2) == -ay;我们可以使用`ode45`来求解微分方程：tspan = [0 10];y0 = 1;params = 2;[t,y] = ode45(@(t,y) -paramsy, tspan, y0);在这个例子中，`@(t,y) -paramsy`定义了微分方程的右侧。
参数`params`在这里是带参数微分方程中的参数。
总之，在MATLAB中解带参数的微分方程，你可以使用`dsolve`函数来获得解析解，或者使用数值求解器如`ode45`来获得数值解。
希望这些信息对你有所帮助。
matlab求解微分方程
matlab求解微分方程Matlab作为一款多功能的计算机科学软件，具有强大的功能，能够有效地解决复杂数学问题。
其中，Matlab特别擅长求解微分方程。
微分方程是一类重要的数学方程，能够解释物体随时间变化的现象，广泛应用于科学领域，如力学、热力学等领域。
下面，将介绍如何使用Matlab来求解微分方程。
对于微分方程，Matlab提供了一系列的函数来支持求解。
常用的函数包括ODE45、ODE23、ODE113，用于解决传递微分方程，以及dsolve用于解决非传递微分方程。
首先，对于传递微分方程，我们可以使用ODE45、ODE23或ODE113函数求解。
其中，ODE45函数是Matlab中一种常用的函数，它可以求解线性和非线性微分方程，而ODE23和ODE113则只能求解线性方程。
要使用这些函数，首先需要准备解决微分方程的参数，包括：要求解的函数、初始条件和积分区间等，可以使用如下Matlab代码实现：%求解的函数f = @(t,y) y+t^2-1;%始条件t0 = 0;y0 = 1;%分区间tspan = [t0,2];%用函数求解[t,y] = ode45(f,tspan,y0);接着，我们就可以调用ODE45、ODE23或ODE113函数求解传递微分方程，代码如下：%用函数求解[t,y] = ode45(f,tspan,y0);%于求解非线性方程[t,y] = ode23(f,tspan,y0);%于求解精确的线性方程[t,y] = ode113(f,tspan,y0);另外，Matlab也提供了一个dsolve函数，用于求解非传递微分方程。
它可以解决一元微分方程、偏微分方程以及系统微分方程等问题。
要使用dsolve函数，可以这样写：%求解的函数f =D2y+4*Dy+3*y = 0’;%解方程y = dsolve(f);以上，就是Matlab中求解微分方程的基本步骤。
可以看到，Matlab提供了一系列函数来支持求解微分方程，让我们不再受到繁琐的数学推导的束缚，使用起来方便快捷。
MATLAB解微分方程
但能求解析解的常微分方程是有限的，大多数的常
微分方程是给不出解析解的.
例
2 2 y x y
这个一阶微分方程就不能用初等函数及其积
分来表达它的解。
例
y y y (0) 1
x x y e , e 的值仍需插值方法来计算. 的解
7
8
事实上，从实际问题当中抽象出来的微分方程，通
15
同理，在x= xn 处，用差商代替导数： y( xn1 ) y( xn ) y( xn1 ) y( xn ) y( xn ) xn1 xn h 由 得 若记
y( xn ) f ( xn , yn )
y( xn1 ) y( xn ) hf ( xn , yn ) y( xn ) yn , y( xn1 ) yn1
26
Two Body PБайду номын сангаасoblem
(t ) u (t ) / r (t ) 3 u (t ) v (t ) / r (t ) 3 v
r (t ) u (t ) 2 v(t ) 2
u (t ) v(t ) y (t ) u (t ) v ( t )
则上式可记为
yn1 yn hf ( xn , yn )
此 即 为 求 解 初 值 问 题 的 Euler 方 法 ， 又 称 显 式 16 Euler方法。
例: 用Euler方法求解常微分方程初值问题
y 2 y 2 y x y (0) 0. (0 x3 )
值积分、数值微分、泰勒展开等离散化方法，对初 值问题
y f ( x, y ) y ( x0 ) y 0
y
用 Matlab 求解微分方程
方式二 输入：[x,y,z]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x -5*y+3*z','Dz=4*x-4*y+2*z', 't')； x=simple(x) % 将x化简 y=simple(y) z=simple(z) 输出：x = C2/exp(t)+C3*exp(t)^2 y = C2*exp(-t)+C3*exp(2*t)+exp(-2*t)*C1 z = C3*exp(2*t)+exp(-2*t)*C1
2
例 8.5.3 求解下列微分方程组
dx dt 2 x 3 y 3 z dy 4 x 5 y 3z dt dz 4 x 4 y 2 z dt
求通解
方式一 输入： [x,y,z]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x -5*y+3*z','Dz=4*x-4*y+2*z', 't')； 输出：x = C2*exp(-t)+C3*exp(2*t) y = C2*exp(-t)+C3*exp(2*t)+exp(-2*t)*C1 z = C3*exp(2*t)+exp(-2*t)*C1
运行程序，得到如图的结果。图中，y1 的图 形为实线， y2 的图形为“ * ”线， y3 的图形为 1 “+”线。
0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0 2 4 6 8 10 12
例 8.5.6 导弹追踪问题
matlab求微分方程精确解与近似解课件
Runge-Kutta 方法
Euler 法与 R-K法误差比较
Matlab 解初值问题
用 Maltab自带函数 解初值问题
求解析解：dsolve
求数值解：
ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb
dsolve 求解析解
dsolve 的使用
y=dsolve('eq1','eq2', ... ,'cond1','cond2', ... ,'v') 其中 y 为输出， eq1、eq2、...为微分方程，cond1、 cond2、...为初值条件，v 为自变量。
例 1：求微分方程 dy 2xyxex2 的通解，并验证。 dx
>> y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') >> syms x; diff(y)+2*x*y - x*exp(-x^2)
y=y+h*subs(f,{'x','y'},{x,y}); x=x+h; szj=[szj;x,y]; end szj plot(szj(:,1),szj(:,2),'or-')
Euler折线法举例（续）
解析解： y53e3x 2x231/3
解析解
近似解
Runge-Kutta 方法
为了减小误差，可采用以下方法：
例
4：求初值问题
dy dx
2 y 2x2
2x
的数值解，求解范
围为 [0,0.5]
y(0) 1
matlab求解微分方程数值解与解析解
matlab求解微分方程数值解与解析解微分方程是数学中的重要内容，它描述了物理、工程、经济等领域中的许多现象和问题。
在实际应用中，我们经常需要求解微分方程的解析解或数值解。
本文将以Matlab为工具，探讨如何求解微分方程并对比解析解与数值解的差异。
一、引言微分方程是描述自然界中许多现象和问题的数学语言，它包含了未知函数及其导数与自变量之间的关系。
微分方程的求解可以帮助我们了解问题的性质和变化规律，并为实际应用提供参考。
在许多情况下，微分方程的解析解很难求得，这时我们可以利用计算机进行数值求解。
二、微分方程的数值解法1.欧拉法欧拉法是最简单的数值求解微分方程的方法之一。
它通过将微分方程转化为差分方程，然后利用离散的点逼近连续的解。
具体步骤如下：（1）将微分方程转化为差分方程，即用近似的导数代替真实的导数；（2）选择初始条件，即确定初始点的值；（3）选择步长和求解区间，即确定求解的范围和步长；（4）使用迭代公式计算下一个点的值；（5）重复步骤（4），直到达到指定的求解区间。
2.改进的欧拉法欧拉法存在精度较低的问题，为了提高精度，可以使用改进的欧拉法。
改进的欧拉法是通过使用两次导数的平均值来计算下一个点的值，从而提高了数值解的精度。
3.龙格-库塔法龙格-库塔法是一种常用的数值求解微分方程的方法，它通过使用多个点的导数来逼近连续解。
龙格-库塔法的步骤如下：（1）选择初始条件和步长；（2）使用迭代公式计算下一个点的值；（3）计算下一个点的导数；（4）根据导数的值和步长计算下一个点的值；（5）重复步骤（3）和（4），直到达到指定的求解区间。
龙格-库塔法的精度较高，适用于求解一阶和高阶微分方程。
三、微分方程的解析解解析解是指能够用公式或函数表示的方程的解。
有些微分方程具有解析解，可以通过数学方法求得。
例如，一阶线性常微分方程和某些特殊类型的二阶微分方程等。
解析解的优势在于精确性和直观性，能够帮助我们深入理解问题的本质。
matlab求解运动微分方程
matlab求解运动微分方程摘要：1.Matlab 求解运动微分方程的概述2.运动微分方程的定义和例子3.Matlab 求解运动微分方程的基本步骤4.Matlab 求解运动微分方程的实例演示5.Matlab 求解运动微分方程的优点和局限性正文：一、Matlab 求解运动微分方程的概述Matlab 是一种广泛应用于科学计算和工程设计的软件，它具有强大的数值计算和数据分析功能。
在运动学和动力学领域，微分方程是描述物体运动状态和变化规律的重要工具，而Matlab 则可以有效地求解这些微分方程，为研究者提供可靠的理论依据。
二、运动微分方程的定义和例子运动微分方程是一类描述物体运动状态和变化规律的偏微分方程。
例如，牛顿第二定律可以表示为：F=ma，其中F 表示物体所受合外力，m 表示物体的质量，a 表示物体的加速度。
在Matlab 中，我们可以通过符号运算和数值计算的方法求解这类微分方程。
三、Matlab 求解运动微分方程的基本步骤1.准备数据：根据实际问题，确定物体的运动状态，如初始速度、初始位置等，并将这些数据转换为Matlab 可以处理的数值形式。
2.建立模型：根据问题的实际情况，建立相应的微分方程模型，如牛顿第二定律模型、简谐振动模型等。
3.编写程序：利用Matlab 的符号运算和数值计算功能，编写求解微分方程的程序。
4.运行程序：执行编写的程序，得到微分方程的解。
5.分析结果：根据求解结果，对物体的运动状态和变化规律进行分析。
四、Matlab 求解运动微分方程的实例演示假设有一个质量为m 的物体在水平面上受到一个随时间变化的力F(t) 的作用，我们可以通过Matlab 求解牛顿第二定律微分方程：1.定义符号变量：m、F(t)、a(t)2.建立微分方程模型：a(t) = F(t)/m3.编写求解程序：使用Matlab 的ode45 函数求解微分方程4.运行程序：得到物体的加速度随时间的变化关系5.分析结果：根据加速度变化关系，分析物体的运动状态和变化规律五、Matlab 求解运动微分方程的优点和局限性Matlab 求解运动微分方程具有以下优点：1.强大的数值计算能力：Matlab 具有丰富的数值计算函数和方法，能够有效地求解微分方程。
微分方程的Matlab求解ppt课件
y1' y2 y3
y2 ' y1 y3 y3' 0.51y1 y2
y1(0) 0, y2 (0) 1, y3(0) 1
解 1、建立m-文件如下：
function dy=rigid(t,y) dy=zeros(3,1); dy(1)=y(2)*y(3); dy(2)=-y(1)*y(3); dy(3)=-0.51*y(1)*y(2);
s.z x=simple(s.x) %简化结果
y=simple(s.y)
z=simple(s.z)
结 果 为：
x =-(-C1-C2*exp(-3*t)+C2-C3+C3*exp(-3*t))*exp(2*t)
y =(-C1*exp(-4*t)+C1+C2*exp(-4*t)+C2*exp(-3*t)-C2+C3-C3*exp(-3*t))*exp(2*t)
x(0) 0, y(0) 0
返回
2. 模型求解
(1) w=20时,建立m-文件如下: function dy=eq3(t,y) dy=zeros(2,1); dy(1)=20*(10+20*cos(t)-y(1))/sqrt ((10+20*cos(t)-y(1))^2+(20+15*sin(t)-y(2))^2); dy(2)=20*(20+15*sin(t)-y(2))/sqrt ((10+20*cos(t)-y(1))^2+(20+15*sin(t)-y(2))^2);
注意: 1、在解n个未知函数的方程组时，x0和x均为n维向量，
m-文件中的待解方程组应以x的分量形式写成.
2、使用Matlab软件求数值解时，高阶微分方程必须 变换成等价的一阶微分方程组.
实验五__用matlab求解常微分方程
实验五 用matlab 求解常微分方程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 =++++--若上式中的系数ni 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反过来，在许多情况下，一阶微分方程组也可化为高阶方程。
matlab求解微分方程特解
一、概述微分方程是描述自然现象和工程问题的数学工具，其中特解是微分方程的解的一种。
而MATLAB是一种高级技术计算语言和交互式环境，被广泛应用于工程、科学和其他领域。
在MATLAB中求解微分方程特解是非常常见的问题，本文将介绍如何使用MATLAB求解微分方程特解。
二、微分方程特解的概念微分方程的一般形式可表示为：dy/dx = f(x, y)其中y是未知函数，x是自变量，f是已知函数。
微分方程的特解是指满足特定初值条件的解，通常表示为y(x0) = y0，其中x0和y0是已知的初值。
三、MATLAB求解微分方程特解的基本步骤1. 定义微分方程在MATLAB中，首先需要定义微分方程的函数形式。
假设我们要求解的微分方程为dy/dx = x + y，则在MATLAB中可以定义函数形式为：function dydx = myfun(x, y)dydx = x + y;2. 定义初值条件接下来需要定义初值条件，即给定的初始条件。
假设初值条件为y(0)= 1，则在MATLAB中可以定义为：x0 = 0;y0 = 1;3. 求解微分方程通过调用MATLAB中的内置函数ode45，可以求解微分方程的特解。
具体的求解过程为：[t, y] = ode45(myfun, [x0, xf], y0);其中myfun表示微分方程的函数形式，[x0, xf]表示求解的自变量范围，y0表示初值条件，t和y分别为求解得到的自变量和特解。
四、示例下面通过一个具体的示例来演示如何使用MATLAB求解微分方程特解。
假设我们要求解的微分方程为dy/dx = x^2 + y，初值条件为y(0) = 1，求解范围为x从0到5。
在MATLAB中定义微分方程的函数形式为：function dydx = myfun(x, y)dydx = x^2 + y;然后定义初值条件为：x0 = 0;y0 = 1;最后调用ode45函数求解微分方程特解：[t, y] = ode45(myfun, [x0, 5], y0);求解得到的自变量和特解分别存储在t和y中，可以通过绘图或其他方式对特解进行进一步分析。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第五讲 Matlab求解微分方程教学目的：学会用MATLAB求简单微分方程的解析解、数值解，加深对微分方程概念和应用的理解；针对一些具体的问题，如追击问题，掌握利用软件求解微分方程的过程；了解微分方程模型解决问题思维方法及技巧；体会微分方程建摸的艺术性．教学重点：利用机理分析建模微分方程模型，掌握追击问题的建模方法，掌握利用MATLAB求解数值解．教学难点：利用机理分析建模微分方程模型，通过举例，结合图形以及与恰当的假设突破教学难点．1 微分方程相关函数（命令）及简介函数名函数功能Dy 表示y 关于自变量的一阶导数D2y 表示y 关于自变量的二阶导数dsolve('equ1','equ2',…) 求微分方程的解析解，equ1、equ2、…为方程（或条件）simplify(s) 对表达式s 使用maple 的化简规则进行化简[r,how]=simple(s) simple 命令就是对表达式s 用各种规则进行化简，然后用r 返回最简形式，how 返回形成这种形式所用的规则．[T,Y] = solver(odefun,tspan,y0) 求微分方程的数值解，其中的solver为命令ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb 之一，odefun 是显式常微分方程：，在积分⎪⎩⎪⎨⎧==)(),(yt yytfdtdy区间tspan=上，从到，用初始条件求解，],[0ftttfty要获得问题在其他指定时间点上的解，则令,21,,ttttspan=（要求是单调的）．],,,[,210fttttezplot(x,y,[tmin,tmax])符号函数的作图命令．x,y 为关于参数t 的符号函数，[tmin,tmax] 为t 的取值范围．inline() 建立一个内联函数．格式：inline('expr', 'var1','var2',…) ，注意括号里的表达式要加引号．因为没有一种算法可以有效地解决所有的ODE 问题，为此，Matlab 提供了多种求解器Solver，对于不同的ODE 问题，采用不同的Solver．求解器 Solver ODE 类型 特点说明ode45非刚性单步算法；4、5阶Runge-Kutta 方程；累计截断误差达 3)(x ∆大部分场合的首选算法ode23 非刚性单步算法；2、3阶Runge-Kutta方程；累计截断误差达 3)(x ∆使用于精度较低的情形ode113 非刚性 多步法；Adams 算法；高低精度均可到6310~10--计算时间比 ode45 短ode23t适度刚性采用梯形算法适度刚性情形ode15s 刚性多步法；Gear's 反向数值微分；精度中等若 ode45 失效时，可尝试使用ode23s 刚性单步法；2阶 Rosebrock 算法；低精度当精度较低时，计算时间比 ode15s 短ode23tb 刚性 梯形算法；低精度当精度较低时，计算时间比 ode15s 短要特别的是：ode23、ode45 是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的 Matlab 的常用程序，其中：ode23 采用龙格-库塔2 阶算法，用3 阶公式作误差估计来调节步长，具有低等的精度．ode45 则采用龙格-库塔4 阶算法，用5 阶公式作误差估计来调节步长，具有中等的精度．2 求解微分方程的一些例子2.1 几个可以直接用 Matlab 求微分方程精确解的例子：例1：求解微分方程，并加以验证． 22x xe xy dxdy-=+求解本问题的Matlab 程序为： syms x y %line1 y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') %line2 diff(y,x)+2 *x*y-x*exp(-x^2) %line3 simplify(diff(y,x)+2*x*y-x*exp(-x^2)) %line4说明：(1) 行line1是用命令定义x,y 为符号变量．这里可以不写，但为确保正确性，建议写上；(2) 行line2是用命令求出的微分方程的解：1/2*exp(-x^2)*x^2+exp(-x^2)*C1(3) 行line3使用所求得的解．这里是将解代入原微分方程，结果应该为0，但这里给出：-x^3*exp(-x^2)-2*x*exp(-x^2)*C1+2*x*(1/2*exp(-x^2)*x^2+exp(-x^2)*C1)(4) 行line4 用 simplify() 函数对上式进行化简，结果为 0， 表明)(x y y =的确是微分方程的解．例2：求微分方程在初始条件下的特解，并画出解函0'=-+x e y xy e y 2)1(=数的图形．求解本问题的 Matlab 程序为： syms x yy=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') ezplot(y)微分方程的特解为：y=1/x*exp(x)+1/x* exp (1) (Matlab 格式)，即，xe e y x+=此函数的图形如图 1：图1 y 关于x 的函数图象2.2 用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)．例３：求解微分方程初值问题的数值解，求解范围为⎪⎩⎪⎨⎧=++-=1)0(2222y xx y dx dy 区间[0, 0.5]．fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,0.5],1); x'; y'; plot(x,y,'o-') >> x' ans =0.0000 0.0400 0.0900 0.1400 0.1900 0.2400 0.2900 0.3400 0.3900 0.4400 0.4900 0.5000>> y' ans =1.0000 0.9247 0.8434 0.7754 0.7199 0.6764 0.6440 0.6222 0.6105 0.6084 0.6154 0.6179 图形结果为图 2．图2 y 关于x 的函数图像3 常微分在实际中的应用 3.1 导弹追踪问题1、符号说明，，乙舰的速率恒为v 0；设时刻乙舰的坐标为，w t ((),())X t Y t 导弹的坐标为；当零时刻，，．建((),())x t y t ((0),(0))(1,0)X Y =((0),(0))(0,0)x y =立微分方程模型．202,01(0)0,'(0)0v d yk x dx w y y ⎧⎪⎪= =<<⎨⎪⎪==⎩由微分方程模型解得11(1)(1)11(),12(1)2(1)2(1)2(1)k k x x y x k k k k k +-+--=-+-≠+-+-++代入题设的数据，得到导弹的运行轨迹为1/5k = 4655555(1)(1)81224y x x =--+-+当时，即当乙舰航行到点处时被导弹击中. 被击中时间1=x 245=y )245,1(为:. 若v 0=1, 则在t =0.21处被击中. 00245v v y t ==利用MALAB 作图如图３．clear, x=0:0.01:1; y=-5*(1-x).^(4/5)/8+5*(1-x).^(6/5)/12+5/24; plot(x,y,'*') 图３导弹运行轨迹（解析法） 图４两种方法对比的导弹运行轨迹2、数值方法求解．设导弹速率恒为，则得到参数方程为w ⎪⎪⎩⎪⎪⎨⎧--+-=--+-=)()()()()()(2222y Y y Y x X wdt dy x X y Y x X w dt dx因乙舰以速度沿直线运动，设，，因此导弹运0v 1x =01v =5,1,w X Y t ===动轨迹的参数方程为：⎪⎪⎪⎩⎪⎪⎪⎨⎧==--+-=--+-=0)0(,0)0()()()1(5)1()()1(52222y x y t y t x dt dyx y t x dt dxMATLAB 求解数值解程序如下，结果见图４ t0=0,tf=0.21;[t,y]=ode45('eq2',[t0 tf],[0 0]); X=1;Y=0:0.001:0.21;plot(X,Y,'-') plot(y(:,1),y(:,2),'*'),hold onx=0:0.01:1; y=-5*(1-x).^(4/5)/8+5*(1-x).^(6/5)/12+5/24; plot(x,y,'r')很明显数值计算的方法比较简单而且适用． 3.2 蚂蚁追击问题（1）建立平面直角坐标系．以正多边形的中心为原点, 设正多边形的一个顶点为起始点, 连接此点和原点作轴. 根据轴作出相应的轴; 选取足够小x x y 的t ∆进行采样．（2）在每一时刻，计算每只蚂蚁在下一时刻时的坐标．不妨设甲追t t t +∆逐对象是乙，在时间时，甲的坐标为，乙的坐标为．甲在t A 11(,)x y B 22(,)x y t t +∆时在点(如图1), 其坐标为'A ，11(cos ,sin )x v t y v t θθ+∆+∆其中. 同理，依次计算下2121cos ,sin ,x x y yd d dθθ--===一只蚂蚁在时的坐标．通过间隔t ∆进行采样，得到新一轮各个蚂蚁在一个t t +∆新的正多边形位置坐标．（4）重复2）步，直到充分小为止．d （5）连接每只蚂蚁在各时刻的位置，就得到所求的轨迹．11,(,)A x y22,(,)B x y 图1 在采样时间内，相连蚂蚁追击用MALAB 求解并作图，函数zhuJi(x,y)在附录一定义，如图６ t=[1:8]; s=7*exp(t.*2*pi/length(t)*i); x=real(s); y=imag(s);zhuJi(x,y)图6 当蚂蚁为7只时的图形习题1. 求微分方程的通解． 0sin 2')1(2=-+-x xy y x2. 求微分方程的通解． x e y y y x sin 5'2''=+-3. 求微分方程组⎪⎪⎩⎪⎪⎨⎧=-+=++00y x dtdy y x dtdx在初始条件下的特解，并画出解函数的图形．0|,1|00====t t y x ()y f x =4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解)，求解区间为．利用画图来比较两种求解器之间的差异．[0,2]t ∈4 参考文献:[1] Mastering Matlab 6, D. Hanselman, B. Littlefield, 清华大学出版,2002[2] 赵静等编．数学建模与数学实验（第3版）．北京：高等教育出版社．2008. [3] 姜启源编. 数学模型（第二版）．北京：高等教育出版社.1993.[4] 石勇国. 蚂蚁追击问题与等角螺线. 宜宾学院学报. 2008，(6): 23-25. [5] 张伟年，杜正东，徐冰．常微分方程．北京：高等教育出版社．2006．5 附录附录一:zhuji(x,y)的M 文件 function zhuji(x,y) clf v=1; dt=0.05;x(length(x)+1)=x(1);y(length(y)+1)=y(1); plot(x,y,'*') holdfor it=1:1000for i=1:length(x)-1 d=sqrt((x(i)-x(i+1))^2+(y(i)-y(i+1))^2); x(i)=x(i)+v*dt*(x(i+1)-x(i))/d; y(i)=y(i)+v*dt*(y(i+1)-y(i))/d; endplot(x,y,'.') hold onx(length(x))=x(1); y(length(y))=y(1); end。