ODE的求解教学讲解课件

ODE的求解MCMContents第2页，共122页1常微分方程求解问题2微分方程matlab求解3Contents第3页，共122页常微分方程的基本概念及其在一些领域的应用实例;One数值解法：Euler法

Runge-Kutta方法；Three求常微分方程解析解的方法以及如何借助数学软件Matlab来实现;TwoSimulink仿真在求解常微分方程上的应用.

FourFive向量场绘图1.常微分方程概念第4页，共122页由自变量、自变量的未知函数以及函数的导数组成的等式叫做常微分方程。自变量自变量的函数函数的导数常微分方程概念（续）第5页，共122页微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数,所以一个n阶常微分方程的一般形式如下:

这里是的已知函数,而且一定含有,y是x的未知函数.在以下内容中把常微分方程简称为“微分方程”或“方程”.

常微分方程概念（续）第6页，共122页n阶线性方程:

这里为x的已知函数。不是线性微分方程的方程称为非线性方程。

常微分方程概念（续）第7页，共122页常微分方程的解初值问题（主要介绍本类问题，下面详细展开）边值问题（与这类问题相关的Matlab解法可以参考bvp4c函数）………………

2.常微分方程的应用第8页，共122页应用常微分方程理论是基础数学的重要组成部分之一，在反映客观现实世界运动过程的各种量之间的关系中，大量存在满足微分方程关系的数学模型.比如在自然科学、经济、生态、人口以及交通等各个领域,某一个（或某几个）量的函数变化规律常用常微分方程（组）来描述，很多弹性物体振动、人口增长、种群竞争等模型都是使用相应的微分方程来描述的.第9页，共122页Logistic人口模型

SIR传染病模型Volterra种群模型ODE模型:常微分方程的应用（续）第10页，共122页Lorenz族系统三分子化学动力学模型一种三神经元神经网络模型ODE模型:常微分方程的应用(续)例.几个常见的微分方程模型第11页，共122页第12页，共122页第13页，共122页第14页，共122页高等数学中，曾经介绍过求一些常微分方程的解，比如：Bernoulli方程和恰当常微分方程(全微分方程)等等。但是对于复杂的方程，求解过程会比较繁琐，如何应用我们已经熟悉的工具——Matlab来求常微分方程的解呢？4.常微分方程的求解？第15页，共122页第一部分：解析解

对一般的常微分方程（下简称方程），一般来讲是不能够求出其解析解，但是对一类常见的特殊的方程——线性常微分方程而言，很容易求出其解.首先来看一般求解线性方程的命令：其中，可以描述微分方程，也可以描述初始条件.注：matlab用Dky表示y对x的k阶导数。系统默认的自变量为t第16页，共122页求解方程

解例解析解应用dsolve命令y=dsolve(‘D2y=sin(t)’)，得到方程的通解y=-sin(t)+C1*t+C2，其中C1和C2为任意常数.求解方程

解例方程的自变量为x，应指明自变量xy=dsolve(‘D2y=cos(x)’,’x’),得到方程的通解y=-cos(x)+C1*x+C2第17页，共122页解例求一右半平面上曲线使其任一点切线在纵轴上截距等于切点的横坐标.设曲线方程为依题意曲线上任意一点都应满足:

应用Matlab求解y=dsolve(‘Dy=y/x-1’,‘x’)，得到y=-x*log(x)+C1*xxx,y第18页，共122页很显然这是一条定义域为正半轴的曲线，当C1取1和-1时图像分别如图所示，可以看出，Matlab能够很方便的求出一些并不直观的常微分方程的解析解.

第19页，共122页

在一个控制系统中，已知控制信号为：求解如下线性常微分方程：例若已知y(0)=3,Dy(0)=2，求其特解.首先来计算方程右侧控制项,键入以下程序:>>symst;u=exp(-t)*cos(t);diff(u,t)+2*u%计算控制项的值得到控制项为exp(-t)*cos(t)-exp(-t)*sin(t)为求方程的解，键入以下命令:y=dsolve(‘D2y+2*Dy-3*y=exp(-t)*cos(t)-exp(-t)*sin(t)’);%求解得到最后结果如下:y=-1/5*exp(-t)*cos(t)+1/5*exp(-t)*sin(t)+C1*exp(t)+C2*exp(-3*t).下面来求在已有初始条件下方程的特解,键入>>y=dsolve('D2y+2*Dy-3*y=exp(-t)*cos(t)-exp(-t)*sin(t)','y(0)=3','Dy(0)=2');%注意初始条件的写法.运行得到:y=-1/5*exp(-t)*cos(t)+1/5*exp(-t)*sin(t)+14/5*exp(t)+2/5*exp(-3*t).

解高阶第20页，共122页解例求下述方程组的解析解

当一个系统含有多个依赖于某自变量的函数的时候，就需要用常微分方程组来描述，对一般的线性微分方程组，仍然可以用dsolve来求解.

求解过程如下：>>[x,y]=dsolve('D2y+2*Dx=x+2*y-exp(-t)','Dy=4*x+5*y+4*exp(-t)')%注意这里得到的结果应为向量[x,y].得到结果：x=-12*(8*exp(-t)+7*C1*(11/12+1/12*193^(1/2))*exp(1/12*(11+193^(1/2))*t)+7*C2*(11/12-1/12*193^(1/2))*exp(-1/12*(-11+193^(1/2))*t))/(23+193^(1/2))/(-23+193^(1/2))+60*(-8*exp(-t)+7*C1*exp(1/12*(11+193^(1/2))*t)+7*C2*exp(-1/12*(-11+193^(1/2))*t))/(23+193^(1/2))/(-23+193^(1/2))-exp(-t)y=-48*(-8*exp(-t)+7*C1*exp(1/12*(11+193^(1/2))*t)+7*C2*exp(-1/12*(-11+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*193^(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=Emsinwt，电阻R和电感L为常量，求电流I（t）.解：显然电流满足如下常微分方程：\omega用w代替DI+R/L*I=Em/L*sinwt设电路接通时刻为t0=0，于是有初始条件I(0)=0在Matlab下求解symsRLEmw;dsolve(‘DI+R/L*I=Em/L*sin(w*t)’,‘I(0)=0’)；得到ans=exp(-R/L*t)/(R^2+w^2*L^2)*Em*w*L-Em*(w*L*cos(w*t)-R*sin(w*t))/(R^2+w^2*L^2)键入pretty(simple(ans))得到第22页，共122页例：高等数学中的例题P318求Euler方程的通解求解命令dsolve('x^3*D3y+x^2*D2y-4*x*Dy=3*x^2','x')得到通解1/3*C1*x^3-1/2*x^2-C2/x+C3与用变换和特征方程的方法求解相比……第23页，共122页第24页，共122页解例求下述方程组的解析解

>>x=dsolve('Dx=x*(1-x^2)')

解得两组通解如下：

x=[1/(1+exp(-2*t)*C1)^(1/2)][-1/(1+exp(-2*t)*C1)^(1/2)]非线性方程部分非线性方程仍可用dsolve求解第25页，共122页一个不能用dsolve求解的例子解例求下述方程组的解析解（Vanderpol方程）>>symsymu>>y=dsolve('D2y+mu*(y^2-1)*Dy+y=0')

运行后得到：Warning:Compact,analyticsolutioncouldnotbefound.ItisrecommendedthatyouapplyPRETTYtotheoutput.Trymhelpdsolve,mhelpRootOf,mhelpDESol,ormhelpallvaluesformoreinformation.…………在这里，dsolve命令不能够求出方程的解，Matlab给出了提示.第26页，共122页方程无解析解数值解

解的性质考察方程随（几）个参数的变化，以及解是如何变化的第二部分：数值解

解的图像第27页，共122页引入新变量替换

首先把一个高阶微分方程化简为一阶微分方程.m+n维的一阶微分方程第28页，共122页解例将Vanderpol方程化为一阶微分方程组常微分方程数值解的提法一般是针对一个一阶方程(组)和初始条件来说的，比如有方程：数值解第29页，共122页求方程数值解的常规方法有Euler法变步长Euler法Runge-Kutta法Runge-Kutta-Felhberg法等Adams-Bashforth-Moulton法

主要介绍Euler法和Runge-Kutta法.第30页，共122页Euler方法在小区间上用差商来代替方程中导数，Euler方法可以分为向前和向后Euler方法.如果

中的t在上选取区间的左端点，则有向前Euler公式：1.基本思想:第31页，共122页t0t1t2近似解xth2.几何意义:解析解h的大小影响计算的速

度与精度;迭代次数的增加会带来

较多的累积误差;向前Euler公式的

精度并不很高提高精度的办法:使用变步长改进Euler公式等1阶精度第32页，共122页

Runge-Kutta方法Euler方法的基本思想，即差商代替导数，很自然的想法是在区间内多取几个点，将它们的斜率加权平均作为导数近似值，这就是Runge-Kutta方法的思想.基本思想:2阶R-K公式:第33页，共122页4阶R-K公式:公式复杂，如何编程实现？？？第34页，共122页在Matlab中，数值解求解函数有：[t,x]=ode45（或ode23）(Fun,[t_0,t_f],x_0)[t,x]=ode45(Fun,[t_0,t_f],x_0,options)[t,x]=ode45(Fun,[t_0,t_f],x_0,options,p_1,…）Matlab函数可以由指定的M-文件给出;[t_0,t_f]为自变量取值区间;options可给出误差限(缺省时相对误差1e-3,绝对误差1e-6);具体形式为：options=odeset(‘reltol’,rt,’abstol’,at).了解例.数学摆的求解在数学摆模型中考虑m=1kg,l=1m,g=9.8m/s^2,初始角度\theta_0=pi/12，初速度为0，阻尼系数\lambda=0.1根据前面描述无阻尼数学摆的方程可以化为如下方程组：第35页，共122页要求这个方程的数值解首先编写M-文件%-----------------------------------------fun.m---------------------------------------functionxdot=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)] %以矩阵形式输出数值解，第一列为时间，第二列为theta值

可以得到数值解图像第36页，共122页数值解可以列表如下：

ttheta

0.00000.26180.00020.2618……省略……9.75420.13819.83370.18729.92780.230410.00000.2503键入：theta0=[pi/12,0];[t,theta]=ode23('fun2',[0,100],theta0);plot(t,theta(:,1));第37页，共122页考虑阻尼后数学摆的方程可以化为如下方程组：为求其数值解，首先编写M-文件如下：%-----------------------------------------fun2.m---------------------------------------functionxdot=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.m……functionxdot=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+e*x(1)+f*x(2))];%……vot.m……[t,x]=ode45('vot',[0,300],x0)plot(x(:,1),x(:,2))第39页，共122页第40页，共122页%……vot.m……functionxdot=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.m……[t,x]=ode45('vot',[0,300],x0)plot(x(:,1),x(:,2))与上图相比，系统的性质发生了变化第41页，共122页首先编写M-文件lorenzeq.m

%-----------------------------------------lorenzeq.m------------------------------------functionxdot=lorenzeq(t,x)xdot=[-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3)];%-----------------------------------------lorenzeq.m------------------------------------解例

画出Lorenz系统的数值解曲线在命令窗口中运行

>>t_final=100;x0=[0;0;1e-10];[t,x]=ode45('lorenzeq',[0,t_final],x0);可以在左侧workspace中看到输出的数值解如下：第42页，共122页续上例第43页，共122页应用plot（t，x）画图得到三条曲线如下图.(分别是t-x1，t-x2，t-x3曲线)续上例第44页，共122页也可绘制三维曲线（x1-x2-x3曲线）如下：>>figure;plot3(x(:,1),x(:,2),x(:,3));得到了Loren系统具有混沌性质的解曲线.应用comet3(x(:,1),x(:,2),x(:,3));可看动态图像，可以看到Lorenz系统具有混沌性的解.续上例第45页，共122页解例*针对上述方程，也可在编程时将参数设计在程序中，通过主窗口下赋不同值，求不同参数下微分方程组的解，此时[]用于站位，不可省略.

首先编写m-文件lorenzleq.m%-----------------------------------------lorenzleq.m---------------------------------------functionxdot=lorenzleq(t,x,flag,beta,rho,sigma)xdot=[-beta*x(1)+x(2)*x(3);-rho*x(2)+rho*x(3);-x(1)*x(2)+sigma*x(2)-x(3)];%-----------------------------------------lorenzleq.m---------------------------------------

在窗口中可以对不同变量赋值，>>t_final=100;x0=[0;0;1e-10];b1=8/3;r1=10;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));subplot(1,3,2);plot(t,x(:,2));subplot(1,3,3);plot(t,x(:,3));第46页，共122页第47页，共122页可以看出，变化Lorenz系统的系数，方程的混沌现象消失了，所以该方法可以便于比较一个系统的解是如何随参数变化的.注意，该方法可以由定义全局变量的方法进行，首先编写M-文件，在M-文件中定义全局变量如下：%-----------------------------------------lorenzlleq.m---------------------------------------functionxdot=lorenzlleq(t,x)globalbetarhosigmaxdot=[-beta*x(1)+x(2)*x(3);-rho*x(2)+rho*x(3);-x(1)*x(2)+sigma*x(2)-x(3)];%-----------------------------------------lorenzlleq.m---------------------------------------键入程序>>globalbetarhosigmabeta=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，1），（-4，4）的5条积分曲线.首先编写M-文件ode.m%-----------------------------------------ode.m---------------------------------------functionxdot=ode1(t,x)xdot=1-x^2;%-----------------------------------------ode.m---------------------------------------第50页，共122页键入：clearc=0.01;x0=-4:.4:4y0=-4:.4:4;[x,y]=meshgrid(x0,y0);%设定绘制向量场的点d=sqrt(1^2+(1-y.^2).^2);u=c./d; %与v一同决定向量场方向v=c*(1-y.^2)./d;holdonquiver(x,y,u,v); %绘制向量场[t,x]=ode23('ode1',[-4,3],-1.000001);plot(t,x,'-r');[t,x]=ode23('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');holdoff第51页，共122页最后可以得到如下图结果，其中有五条积分曲线，分别为过初值（-4，-1.000001），（-4，-1），（-4，-0.8），（-4，1），（-4，4）的微分方程的解，每一条解都是初始值点沿向量场方向运动得到的，具体地说，以y（0）大于-1为初值的解沿向量场走向最终将趋近于1，其他的所有解显然会趋近于负无穷大.例：Logistic模型决定的向量场第52页，共122页Logistic人口模型

（续上例）参考程序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;holdonquiver(x,y,u,v);第53页，共122页第54页，共122页Simulink6交互式仿真

应用Matlab功能强大的仿真工具箱，我们可以从系统的角度对一个常微分方程（组）进行数值模拟，具体方法介绍如下：

第一步，启动Simulink模块可以在主窗口下输入simulink或Simulink3启动或者在Matlab主窗口下从按钮启动.启动后可以看到如下界面：

第55页，共122页Simulink6交互式仿真编写simulink程序与M-file不同，是在仿真编辑窗口中完成的，打开编辑器的方法有：1.主窗口下File-New-Newmodel2.Simulink界面下点击Createanewmodel第56页，共122页对一个常微分方程，通常的办法是将其转化为相应的积分方程来进行模拟，也就是说，对系统

它等价于积分方程：Simulink6交互式仿真下面来看对一个微分方程（组）的仿真计算第57页，共122页2.Simulink中常用的模块简介如下：积分模块求导模块绝对值积系数(增益)求和导入工作间自定义函数第58页，共122页3.各模块都可以自由调整相关参数。初始值比如积分模块中