1、ODE的求解Contents第2页，共122页1常微分方程求解问题常微分方程求解问题2微分方程微分方程 matlabmatlab 求解求解3Contents第3页，共122页常微分方程的基本概念及其在一些领常微分方程的基本概念及其在一些领域的应用实例域的应用实例;One数值解法：数值解法：EulerEuler法法 Runge-KuttaRunge-Kutta方法；方法；Three求常微分方程解析解的方法以及如何求常微分方程解析解的方法以及如何借助数学软件借助数学软件MatlabMatlab来实现来实现;TwoSimulinkSimulink仿真在求解常微分方程上仿真在求解常微分方程上的应用的应
2、用.FourFiveFive向量场绘图向量场绘图1.常微分方程概念第4页，共122页 由自变量、自变量的未知函数以及函数的导数组成的由自变量、自变量的未知函数以及函数的导数组成的等式叫做常微分方程。等式叫做常微分方程。,x tf t x t自自变变量量自自变变量量的的函函数数函函数数的的导导数数常微分方程概念（续）第5页，共122页 微分方程中出现的未知函数最高阶导数的阶微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数数称为微分方程的阶数,所以一个所以一个n n阶常微分阶常微分方程的一般形式如下方程的一般形式如下:,0nndyd yF x ydxdx 这里这里 是是 的已知函数的已知
3、函数,而且一定含而且一定含有有 ,y y是是x x的未知函数的未知函数.在以下内容中把常微在以下内容中把常微分方程简称为分方程简称为“微分方程微分方程”或或“方程方程”.,nndyd yF x ydxdx,nndyd yx ydxdxnndydx常微分方程概念（续）第6页，共122页n n阶线性方程阶线性方程 :这里这里 为为x x的已知函数的已知函数 。不是线性微。不是线性微分方程的方程称为非线性方程。分方程的方程称为非线性方程。1111nnnnnnd ydydyaxaxax yf xdxdxdx 1,naxaxf x常微分方程概念（续）第7页，共122页常微分方程的解常微分方程的解初值问题
4、（主要介绍本类问题，下面详细展开）初值问题（主要介绍本类问题，下面详细展开）边值问题（与这类问题相关的边值问题（与这类问题相关的MatlabMatlab解法可以解法可以参考参考bvp4cbvp4c函数）函数）2.常微分方程的应用第8页，共122页应用应用常微分方程理论是基础数学的重要组成部分之一，在反映客常微分方程理论是基础数学的重要组成部分之一，在反映客观现实世界运动过程的各种量之间的关系中，大量存在满足观现实世界运动过程的各种量之间的关系中，大量存在满足微分方程关系的数学模型微分方程关系的数学模型.比如在比如在自然科学、经济、生态、自然科学、经济、生态、人人口以及交通口以及交通等各个领域等
5、各个领域,某一个（或某几个）量的函数变化某一个（或某几个）量的函数变化规规律常用常微分方程（组）来描述，很多弹性物体振动、人口律常用常微分方程（组）来描述，很多弹性物体振动、人口增长、种群竞争等模型都是使用相应的微分方程来描述的增长、种群竞争等模型都是使用相应的微分方程来描述的.第9页，共122页1mNNrNN 000;0 xkxylxykxyxxynx xx abxcyyy dexfy常微分方程的应用（续）常微分方程的应用（续）第10页，共122页xa yxyxzcxyzxybz 221xABxx yyBxx y tanhtanhtanhxxazyyaxzzay 常微分方程的应用常微分方程的
6、应用(续续)例.几个常见的微分方程模型第11页，共122页第12页，共122页第13页，共122页第14页，共122页 高等数学中，曾经介绍过求一些常微分方程的解，比如：和等等。但是对于复杂的方程，求解过程会比较繁琐，如何应用我们已经熟悉的工具Matlab来求常微分方程的解呢？第15页，共122页第一部分：解析解 对一般的常微分方程（下简称方程），一般来讲是不能够对一般的常微分方程（下简称方程），一般来讲是不能够求出其解析解，但是对一类求出其解析解，但是对一类的特殊的方程的特殊的方程线性常微分方线性常微分方程而言，很容易求出其解程而言，很容易求出其解.首先来看一般求解线性方程的命令：首先来看一
7、般求解线性方程的命令：if其中，其中，可以描述微分方程，也可以描述初始条件可以描述微分方程，也可以描述初始条件.12dsolve(,)myf ff12dsolve(,)myf ffx第16页，共122页求解方程求解方程 解解例例解析解解析解 22sind ytdt应用应用dsolvedsolve命令命令y=y=dsolvedsolve(D2y=sin(t)(D2y=sin(t)，得到方程的通解得到方程的通解y=-sin(t)+C1y=-sin(t)+C1*t+C2t+C2，其中其中C1C1和和C2C2为任意常数为任意常数.求解方程求解方程 解解例例方程的自变量为方程的自变量为 x x，应指明自
8、变量，应指明自变量x xy=y=dsolvedsolve(D2y=(D2y=coscos(x),x),(x),x),得到方程的通解得到方程的通解y=-y=-coscos(x)+C1(x)+C1*x+C2 x+C2 22cosd yxdx第17页，共122页解解例例求一右半平面上曲线使其任一点切线在纵轴上截距等求一右半平面上曲线使其任一点切线在纵轴上截距等于切点的横坐标于切点的横坐标.设曲线方程为 ,yf x依题意曲线上任意一点都应满足:dyyxxdx应用Matlab求解 y=dsolve(Dy=y/x-1,x)，得到y=-x*log(x)+C1*xxx,y第18页，共122页很显然这是一条定义
9、域为正半轴的曲线，当很显然这是一条定义域为正半轴的曲线，当C1C1取取1 1和和-1-1时图像时图像分别如图所示，可以看出，分别如图所示，可以看出，MatlabMatlab能够很方便的求出一些并不能够很方便的求出一些并不直观的常微分方程的解析解直观的常微分方程的解析解.1lnyxxC x 11C 11C 第19页，共122页 在一个控制系统中，已知控制信号为：在一个控制系统中，已知控制信号为：求解如下线性常微分方程：求解如下线性常微分方程：例例 cos()tu tet 232y ty tyu tu t若已知若已知y(0)=3,Dy(0)=2y(0)=3,Dy(0)=2，求其特解，求其特解.首先
10、来计算方程右侧控制项首先来计算方程右侧控制项,键入以下程序键入以下程序:syms t;syms t;u=exp(-t)u=exp(-t)*cos(t);cos(t);diff(u,t)+2diff(u,t)+2*u%u%计算控制项的值计算控制项的值得到控制项为得到控制项为exp(-t)exp(-t)*cos(t)-exp(-t)cos(t)-exp(-t)*sin(t)sin(t)为求方程的解，键入以下命令为求方程的解，键入以下命令:y=dsolve(D2y+2y=dsolve(D2y+2*Dy-3Dy-3*y y=exp(-t)=exp(-t)*cos(t)-exp(-t)cos(t)-ex
11、p(-t)*sin(t);%sin(t);%求解求解得到最后结果如下得到最后结果如下:y=-1/5y=-1/5*exp(-t)exp(-t)*cos(t)+1/5cos(t)+1/5*exp(-t)exp(-t)*sin(t)+C1sin(t)+C1*exp(t)+C2exp(t)+C2*exp(-3exp(-3*t).t).下面来求在已有初始条件下方程的特解下面来求在已有初始条件下方程的特解,键入键入y=dsolve(y=dsolve(D2y+2D2y+2*Dy-3Dy-3*y y=exp(-t)=exp(-t)*cos(t)-exp(-t)cos(t)-exp(-t)*sin(t),y(0
12、)=3,Dy(0)=2);sin(t),y(0)=3,Dy(0)=2);%注意初始条件的写法注意初始条件的写法.运行得到运行得到:y=-1/5y=-1/5*exp(-t)exp(-t)*cos(t)+1/5cos(t)+1/5*exp(-t)exp(-t)*sin(t)+14/5sin(t)+14/5*exp(t)+2/5exp(t)+2/5*exp(-3exp(-3*t).t).解解高阶第20页，共122页解解例例求下述方程组的解析解求下述方程组的解析解 当一个系统含有多个依赖于某自变量的函数的时候当一个系统含有多个依赖于某自变量的函数的时候，就需要用常微分方程组就需要用常微分方程组来描述来
13、描述，对一般的线性微分方程组对一般的线性微分方程组，仍然可以用仍然可以用dsolvedsolve来求解来求解.22454ttx tx tx ty tey tx ty te求解过程如下求解过程如下：x,y=dsolve(D2y+2 x,y=dsolve(D2y+2*Dx=x+2Dx=x+2*y-exp(-t),Dy=4y-exp(-t),Dy=4*x+5x+5*y+4y+4*exp(-t)exp(-t)%注意这里得到的注意这里得到的结果应为向量结果应为向量x,y.x,y.得到结果：得到结果：x=-12x=-12*(8(8*exp(-t)+7exp(-t)+7*C1C1*(11/12+1/12(1
14、1/12+1/12*193(1/2)193(1/2)*exp(1/12exp(1/12*(11+193(1/2)(11+193(1/2)*t)t)+7+7*C2C2*(11/12-1/12(11/12-1/12*193(1/2)193(1/2)*exp(-1/12exp(-1/12*(-11+193(1/2)(-11+193(1/2)*t)/(23+193(1/2)/(-t)/(23+193(1/2)/(-23+193(1/2)+6023+193(1/2)+60*(-8(-8*exp(-t)+7exp(-t)+7*C1C1*exp(1/12exp(1/12*(11+193(1/2)(11+19
15、3(1/2)*t)+7t)+7*C2C2*exp(-1/12exp(-1/12*(-(-11+193(1/2)11+193(1/2)*t)/(23+193(1/2)/(-23+193(1/2)-exp(-t)t)/(23+193(1/2)/(-23+193(1/2)-exp(-t)y=-48y=-48*(-8(-8*exp(-t)+7exp(-t)+7*C1C1*exp(1/12exp(1/12*(11+193(1/2)(11+193(1/2)*t)+7t)+7*C2C2*exp(-1/12exp(-1/12*(-11+193(1/2)(-11+193(1/2)*t)/(23+193(1/2)
16、/(-23+193(1/2)t)/(23+193(1/2)/(-23+193(1/2)续上例结果很复杂，对其进行化简：x=simple(x)y=simple(y)得到结果：x=5/7*exp(-t)-9/48*C1*exp(11/12*t+1/12*t*193(1/2)+1/48*C1*exp(11/12*t+1/12*t*193(1/2)*193(1/2)-49/48*C2*exp(11/12*t-1/12*t*193(1/2)-1/48*C2*exp(11/12*t-1/12*t*193(1/2)*193(1/2)y=-8/7*exp(-t)+C1*exp(11/12*t+1/12*t*1
17、93(1/2)+C2*exp(11/12*t-1/12*t*193(1/2)还可以得到近似的结果：x=vpa(x，4)%保留4位有效数字y=vpa(y，4)x=.7143*exp(-1.*t)-.7314*C1*exp(2.074*t)-1.310*C2*exp(-.2410*t)y=-1.143*exp(-1.*t)+.9998*C1*exp(2.074*t)+.9998*C2*exp(-.2410*t)其中，C1、C2为Matlab给出的任意常数.第21页，共122页例：高等数学中的例题 例2.P278一电路（电源串联一电阻R和一电感L），电源电动势为E=Emsin wt，电阻R和电感L为
18、常量，求电流I（t）.解：显然电流满足如下常微分方程：omega用w代替DI+R/L*I=Em/L*sin wt设电路接通时刻为t0=0，于是有初始条件I(0)=0在Matlab下求解syms R L Em w;dsolve(DI+R/L*I=Em/L*sin(w*t),I(0)=0)；得到ans=exp(-R/L*t)/(R2+w2*L2)*Em*w*L-Em*(w*L*cos(w*t)-R*sin(w*t)/(R2+w2*L2)键入pretty(simple(ans)得到sinmERIItLL第22页，共122页例：高等数学中的例题 P318求Euler方程的通解求解命令dsolve(x3
19、*D3y+x2*D2y-4*x*Dy=3*x2,x)得到通解1/3*C1*x3-1/2*x2-C2/x+C3与用变换和特征方程的方法求解相比32243x yx yxyx第23页，共122页第24页，共122页解解例例求下述方程组的解析解求下述方程组的解析解 x=x=dsolvedsolve(DxDx=x=x*(1-x2)(1-x2)解得两组通解如下：解得两组通解如下：x=x=1/(1+exp(-2 1/(1+exp(-2*t)t)*C1)(1/2)C1)(1/2)-1/(1+exp(-2 -1/(1+exp(-2*t)t)*C1)(1/2)C1)(1/2)2(1)x tx tx t非线性方程非
20、线性方程部分非线性方程仍可用部分非线性方程仍可用dsolvedsolve求解求解第25页，共122页一个不能用一个不能用dsolvedsolve求解的例子求解的例子解解例例求下述方程组的解析解（求下述方程组的解析解（Van der pol Van der pol 方程）方程）symssyms y mu y mu y=y=dsolvedsolve(D2y+mu(D2y+mu*(y2-1)(y2-1)*Dy+yDy+y=0)=0)运行后得到：运行后得到：Warning:Compact,analytic solution could not be found.It is Warning:Compac
21、t,analytic solution could not be found.It is recommended that you apply PRETTY to the output.Try recommended that you apply PRETTY to the output.Try mhelpmhelp dsolvedsolve,mhelpmhelp RootOfRootOf,mhelpmhelp DESolDESol,or,or mhelpmhelp allvaluesallvalues for more for more information.information.210
22、y tyy ty在这里，在这里，dsolvedsolve命令不能够求出方程的解，命令不能够求出方程的解，MatlabMatlab给出了提示给出了提示.第26页，共122页方程无解析解方程无解析解数值解数值解 解的性质解的性质 考察方程随（几）个参考察方程随（几）个参数的变化，以及解是如数的变化，以及解是如何变化的何变化的第二部分：数值解 解的图像解的图像第27页，共122页 1111(,)(,)mmnnmnxf t x xxy yyyg t x xxy yy引入新变量替换引入新变量替换 11121,mnmmm nxx xxxxxyxy 首先把一个高阶微分方程化简为一阶微分方程.12231212
23、12(,)(,)mm nmmm nm nxxxxxf t x xxxxxg t x xxm m+n n维维的的一一阶阶微微分分方方程程第28页，共122页 210y tyy ty解解例例将将Van der pol Van der pol 方程化为一阶微分方程组方程化为一阶微分方程组 1222121(1)xxxxxx 常微分方程数值解的提法一般是针对一个一阶方常微分方程数值解的提法一般是针对一个一阶方程程(组组)和初始条件来说的，和初始条件来说的，比如有方程：比如有方程：00,().x tf t x tx tx12,xy xy 数值解数值解第29页，共122页求方程数值解的常规方法求方程数值解的
24、常规方法有有EulerEuler法法变步长变步长EulerEuler法法Runge-KuttaRunge-Kutta法法Runge-Kutta-FelhbergRunge-Kutta-Felhberg法等法等Adams-Adams-BashforthBashforth-Moulton-Moulton法法 主要介绍主要介绍EulerEuler法和法和Runge-KuttaRunge-Kutta法法.第30页，共122页,f t xEulerEuler方法方法在小区间在小区间 上用差商来代替方程中导数，上用差商来代替方程中导数，EulerEuler方法可以分为向前和向后方法可以分为向前和向后Eul
25、erEuler方法方法.1,nnt t如果如果 中的中的 t t 在在 上选取区间的左端上选取区间的左端点，则有向前点，则有向前EulerEuler公式：公式：1,nnt t1,0,1,nnnnxxhf txnx 1nnx tx th第31页，共122页0(0)xx,nnhf t xt0t1t2近似解近似解xth1x()x h解析解解析解影响计算的速影响计算的速 度与精度度与精度;迭代次数的增加会带来迭代次数的增加会带来 较多的累积误差较多的累积误差;向前向前EulerEuler公式的公式的 精度并不很高精度并不很高提高精度的办法提高精度的办法:使用变步长使用变步长 改进改进EulerEule
26、r公式等公式等1 1阶精度阶精度第32页，共122页 Runge-KuttaRunge-Kutta方法方法EulerEuler方法的基本思想，即差商代替导数，很自然的想方法的基本思想，即差商代替导数，很自然的想法是在区间内多取几个点，将它们的斜率加权平均作法是在区间内多取几个点，将它们的斜率加权平均作为导数近似值，这就是为导数近似值，这就是Runge-KuttaRunge-Kutta方法的思想方法的思想.11 122121,nnnnnnxxhkkkf txkf th xhk1221,0.5,其中：第33页，共122页1123412132431226,0.5,0.50.5,0.5,nnnnnnn
27、nnnxxh kkkkkf txkf th xhkkf th xhkkf th xhk第34页，共122页在在MatlabMatlab 中，数值解求解函数有：中，数值解求解函数有：t,xt,x=ode45=ode45（或（或ode23ode23）(Fun,t_0,t_f,x_0)(Fun,t_0,t_f,x_0)t,xt,x=ode45 (Fun,t_0,t_f,x_0,options)=ode45 (Fun,t_0,t_f,x_0,options)t,xt,x=ode45 (Fun,t_0,t_f,x_0,options,p_1,=ode45 (Fun,t_0,t_f,x_0,options
28、,p_1,）MatlabMatlab函数可以由指定的函数可以由指定的M-M-文件给出文件给出;t_0,t_ft_0,t_f为自变量取值区间为自变量取值区间;optionsoptions可给出误差限可给出误差限(缺省时相对误差缺省时相对误差1e-3,1e-3,绝对误差绝对误差1e-6);1e-6);具体形式为：具体形式为：options=options=odesetodeset(reltol,rt,abstol,atreltol,rt,abstol,at).).了解了解例.数学摆的求解在数学摆模型中考虑m=1kg,l=1m,g=9.8m/s2,初始角度theta_0=pi/12，初速度为0，阻尼
29、系数lambda=0.1根据前面描述无阻尼数学摆的方程可以化为如下方程组：第35页，共122页要求这个方程的数值解首先编写要求这个方程的数值解首先编写M-M-文件文件%-fun.m-function xdot=fun(t,theta)g=9.8;l=1;xdot=theta(2);-g*sin(theta(1)/l;%theta(1)为 theta(2)为%-fun.m-键入命令：theta0=pi/12,0;t,theta=ode23(fun,0,10,theta0);%求解，注意初值的选取，theta0一定是2维向量plot(t,theta(:,1);%绘出积分曲线图t,theta(:,1
30、)%以矩阵形式输出数值解，第一列为时间，第二列为theta值 sin/gl sin/gl 可以得到数值解图像第36页，共122页数值解可以列表如下：t theta 0.0000 0.2618 0.0002 0.2618 省略 9.7542 0.1381 9.8337 0.1872 9.9278 0.230410.0000 0.2503键入：theta0=pi/12,0;t,theta=ode23(fun2,0,100,theta0);plot(t,theta(:,1);第37页，共122页考虑阻尼后数学摆的方程可以化为如下方程组：考虑阻尼后数学摆的方程可以化为如下方程组：sin/glm 为求其
31、数值解，首先编写M-文件如下：%-fun2.m-function xdot=fun2(t,theta)g=9.8;l=1;m=1;lambda=0.1;xdot=theta(2);-g*sin(theta(1)/l-lambda*theta(2)/m;%-fun2.m-第38页，共122页可以画出数值解的图像如下可以画出数值解的图像如下：可以看出，在阻力作用下，数学摆的运动最后将趋于静止.例：Voterra系统%vot.mfunction xdot=vot(t,x)a=1;c=-.1;d=-.5;e=.02;f=0;b=0;xdot=x(1)*(a+b*x(1)+c*x(2);x(2)*(d+
32、e*x(1)+f*x(2);%vot.mt,x=ode45(vot,0,300,x0)plot(x(:,1),x(:,2)第39页，共122页xx abxcyyy dexfy第40页，共122页%vot.mvot.mfunction xdot=vot(t,x)a=1;c=-.1;d=-.5;e=.02;f=-.001;b=-.001;xdot=x(1)*(a+b*x(1)+c*x(2);x(2)*(d+e*x(1)+f*x(2);%vot.mvot.m t,xt,x=ode45(=ode45(votvot,0,300,x0),0,300,x0)plot(x(:,1),x(:,2)plot(x(
33、:,1),x(:,2)与上图相比，系统的性质发生了变化第41页，共122页首先编写首先编写M-M-文件文件lorenzeq.mlorenzeq.m%-%-lorenzeqlorenzeq.m-.m-function xdot=lorenzeq(t,x)function xdot=lorenzeq(t,x)xdotxdot=-8/3=-8/3*x(1)+x(2)x(1)+x(2)*x(3);-10 x(3);-10*x(2)+10 x(2)+10*x(3);-x(1)x(3);-x(1)*x(2)+28x(2)+28*x(2)-x(3);x(2)-x(3);%-%-lorenzeqlorenze
34、q.m-.m-解解例例112322331223xxx xxxxxx xxx 画出画出LorenzLorenz系统系统的数值解曲线的数值解曲线 101238/3,10,28,000,010 xxx在命令窗口中运行在命令窗口中运行 t_finalt_final=100;x0=0;0;1e-10;=100;x0=0;0;1e-10;t,xt,x=ode45(=ode45(lorenzeqlorenzeq,0,t_final,x0);,0,t_final,x0);可以在左侧可以在左侧workspaceworkspace中看到输出的数值解如下：中看到输出的数值解如下：第42页，共122页续上例续上例第4
35、3页，共122页应用应用plotplot（t t，x x）画图得到三条曲线如下图）画图得到三条曲线如下图.(.(分别是分别是t-t-x x1 1，t-xt-x2 2，t-xt-x3 3曲线曲线)续上例续上例第44页，共122页也可绘制三维曲线（也可绘制三维曲线（x x1 1-x-x2 2-x-x3 3曲线）如下：曲线）如下：figure;plot3(x(:,1),x(:,2),x(:,3);figure;plot3(x(:,1),x(:,2),x(:,3);得到了得到了LorenLoren系统具有混沌性质的解曲线系统具有混沌性质的解曲线.应用应用comet3(x(:,1),x(:,2),x(:
36、,3);comet3(x(:,1),x(:,2),x(:,3);可看动态图像，可以看到可看动态图像，可以看到LorenzLorenz系统具有混沌性系统具有混沌性的解的解.续上例续上例第45页，共122页解解例例*针对上述方程，也可在编程时将参数设计在程序中，通过针对上述方程，也可在编程时将参数设计在程序中，通过主窗口下赋不同值，求不同参数下微分方程组的解，此时主窗口下赋不同值，求不同参数下微分方程组的解，此时 用于站位，不可省略用于站位，不可省略.首先编写首先编写m-m-文件文件lorenzleq.mlorenzleq.m%-%-lorenzleq.mlorenzleq.m-function
37、function xdotxdot=lorenzleqlorenzleq(t,x,flag,beta,rho,sigmat,x,flag,beta,rho,sigma)xdotxdot=-beta=-beta*x(1)+x(2)x(1)+x(2)*x(3);-rhox(3);-rho*x(2)+rhox(2)+rho*x(3);-x(1)x(3);-x(1)*x(2)+sigmax(2)+sigma*x(2)-x(3);x(2)-x(3);%-%-lorenzleq.mlorenzleq.m-在窗口中可以对不同变量赋值，t_final=100;x0=0;0;1e-10;b1=8/3;r1=10
38、;s1=28;%注意：这里的参数名不要求与M-文件中一致 t,x=ode45(lorenzleq,0,t_final,x0,b1,r1,s1);%用来占位subplot(1,3,1);plot(t,x(:,1);subplot(1,3,2);plot(t,x(:,2);subplot(1,3,3);plot(t,x(:,3);得到的图像与前面相同.变化参数可以不用修改程序主体得到不同的数值解.例如改变参数值为b1=8;r1=5;s1=5.重新运行:t,x=ode45(lorenzleq,0,t_final,x0,b1,r1,s1);subplot(1,3,1);plot(t,x(:,1);su
39、bplot(1,3,2);plot(t,x(:,2);subplot(1,3,3);plot(t,x(:,3);第46页，共122页第47页，共122页可以看出，变化Lorenz系统的系数，方程的混沌现象消失了，所以该方法可以便于比较一个系统的解是如何随参数变化的.注意，该方法可以由定义全局变量的方法进行，首先编写M-文件，在M-文件中定义全局变量如下：%-lorenzlleq.m-function xdot=lorenzlleq(t,x)global beta rho sigmaxdot=-beta*x(1)+x(2)*x(3);-rho*x(2)+rho*x(3);-x(1)*x(2)+s
40、igma*x(2)-x(3);%-lorenzlleq.m-键入程序global beta rho sigmabeta=8/3;rho=10;sigma=28;%为参数赋值x0=0;0;1e-10;t,x=ode45(lorenzleq,0,100,x0)plot(t,x)也可得到相同结果.第48页，共122页第49页，共122页向量场的绘制*前面曾提到过，Matlab中提供了命令quiver来绘制向量场，具体格式为quiver（x,y,u,v），x和y表示向量起点，u和v确定向量的方向。例例绘出 所确定的向量场，以及过（-4，-1.000001），（-4，-1），（-4，-0.8），（-4，
41、1），（-4，4）的5条积分曲线.21x txt 首先编写M-文件ode.m%-ode.m-function xdot=ode1(t,x)xdot=1-x2;%-ode.m-第50页，共122页键入：clearc=0.01;x0=-4:.4:4y0=-4:.4:4;x,y=meshgrid(x0,y0);%设定绘制向量场的点d=sqrt(12+(1-y.2).2);u=c./d;%与v一同决定向量场方向v=c*(1-y.2)./d;hold onquiver(x,y,u,v);%绘制向量场t,x=ode23(ode1,-4,3,-1.000001);plot(t,x,-r);t,x=ode23
42、(ode1,-4,3,-1);plot(t,x,-r);t,x=ode23(ode1,-4,3,-0.8);plot(t,x,-r);t,x=ode23(ode1,-4,3,1);plot(t,x,-r);t,x=ode23(ode1,-4,3,4);plot(t,x,-r);hold off第51页，共122页最后可以得到如下图结果，其中有五条积分曲线，分别为过初值（-4，-1.000001），（-4，-1），（-4，-0.8），（-4，1），（-4，4）的微分方程的解，每一条解都是初始值点沿向量场方向运动得到的，具体地说，以y（0）大于-1为初值的解沿向量场走向最终将趋近于1，其他的所有解
43、显然会趋近于负无穷大.例：Logistic 模型决定的向量场第52页，共122页1,50,.1mmNNrN NrN（续上例）参考程序clearc=0.01;x0=0:4:100;y0=0:4:100;x,y=meshgrid(x0,y0);%设定绘制向量场的点d=sqrt(1+(0.1*(1-y/50).*y).2);u=c./d;%与v一同决定向量场方向v=c*(0.1*(1-y/50).*y)./d;hold onquiver(x,y,u,v);第53页，共122页第54页，共122页Simulink6 交互式仿真 应用应用MatlabMatlab 功能强大的仿真工具箱，我们可以从功能强大
44、的仿真工具箱，我们可以从 系统的角度对一个常微系统的角度对一个常微分方程（组）进行数值模拟，具体方法介绍如下分方程（组）进行数值模拟，具体方法介绍如下：第一步，启动第一步，启动SimulinkSimulink模块可以在主窗口下输入模块可以在主窗口下输入 simulink simulink 或或Simulink3Simulink3启动启动或者在或者在MatlabMatlab主窗口下从按钮主窗口下从按钮 启动启动.启动后可以看到如下界面：启动后可以看到如下界面：第55页，共122页Simulink6 交互式仿真编写编写simulinksimulink程序与程序与M-fileM-file不同，是在仿
45、真编辑窗口中完成的，打开编辑器的方法有：不同，是在仿真编辑窗口中完成的，打开编辑器的方法有：1.1.主窗口下主窗口下 File-New-New modelFile-New-New model2.Simulink2.Simulink界面下点击界面下点击 Create a new modelCreate a new model第56页，共122页1.1.对一个常微分方程，通常的办法是将其转化为相应对一个常微分方程，通常的办法是将其转化为相应的积分方程来进行模拟，也就是说，对系统的积分方程来进行模拟，也就是说，对系统 它等价于它等价于 积分方程：积分方程：0,.ttx tf s x s 00,().
46、x tf t x tx txSimulink6 交互式仿真下面来看对一个微分方程（组）的仿真计算第57页，共122页2.2.SimulinkSimulink中常用的模块简介如下：中常用的模块简介如下：积分模块积分模块求导模块求导模块绝对值绝对值积积系数系数(增益增益)求和求和导入工作间导入工作间自定义函数自定义函数第58页，共122页3.3.各模块都可以自由调整相关参数。各模块都可以自由调整相关参数。初始值初始值比如积分模比如积分模块中可以设块中可以设定初始值定初始值 第59页，共122页 系数模块中可以设定值系数模块中可以设定值 系数值系数值第60页，共122页4.4.设定仿真环境设定仿真环
47、境 仿真时长仿真时长误差限误差限第61页，共122页Start 第四步，运行第四步，运行(Ctr+TCtr+T或运行图标或运行图标)第五步，输出结果（数组（第五步，输出结果（数组（to workspaceto workspace模块）、模块）、绘图等），主要命令（绘图等），主要命令（plot,plot3,subplotplot,plot3,subplot）.例.简单的差分方程模拟斐波那契数列(第一个月1对小兔子开始，每个月的兔子数量)%rabbit.mdl第62页，共122页例 cos(),00yxy 常微分方程第63页，共122页例 求解Dy=-alpha*y,y(0)=100,0100yy
48、 y 常微分方程第64页，共122页yDy输出仿真时间第65页，共122页第66页，共122页例：求解例：求解Van Van derder polpol 方程方程 1222121(1)xxxxxx 第67页，共122页=-0.5续上例续上例第68页，共122页=-0.5 ，方程形式简化为，方程形式简化为 ，解为，解为 1221xxxx 22121xx续上例续上例第69页，共122页=0.5，周期解，周期解续上例续上例例.求数学摆的数值解（无阻尼）第70页，共122页其中eta初始值为0，theta初始值为pi/12.Fcn模块对应第二个方程右端形式为-9.8*sin(u)/1.仿真参数按如下设
49、置：第71页，共122页例.求数学摆的数值解（有阻尼）第72页，共122页与上例相比只需相应加入 一项./m第73页，共122页求解*：Dx=-0.1*t*x0.1xtx 第74页，共122页作业II1.给出适当参数绘制出SIR模型的解并给出适当的解释;2.#熟悉ODE的解析解与数值解的Matlab求解方法；3.了解一维方程向量场绘图；自己给出适当方程，绘图、观察积分曲线的走向；4.（1）求解Van de pol方程，当$mu=1$时；（2）*尝试在matlab下输入以下命令：t,y=ode45(vdp1,0 20,2 0);%句柄绘图，vdp1是一个函数 Subplot(1,2,1);plot(t,y(:,1);Subplot(1,2,2);plot(y(:,1),y(:,2);观察结果，解释.Hint:help vdp15.#试着针对Lorenz方程组搭建SImulink仿真模块，完成本文中例题相同的工作；6.#练习微分方程的求解（可以在高等数学中找方程或自己随意给出）。7.在练习数学建模基础的同时，查阅文献，找出至少5中人口（或种群）微分方程模型，并求解.第75页，共122页Matlab 学习1.找一本资料2.help3.Internet4.type第76页，共122页77Thank you!Thank you!
