ODE的求解 MCM 2012.10.29 Edit by niuben

Contents 1 2 常微分方程求解问题 3 微分方程 matlab 求解

Contents One Two Five 向量场绘图 Three Four 求常微分方程解析解的方法以及如何 常微分方程的基本概念及其在一些领 域的应用实例; Two 求常微分方程解析解的方法以及如何 借助数学软件Matlab来实现; Five 向量场绘图 数值解法：Euler法 Runge-Kutta方法； Simulink仿真在求解常微分方程上 的应用. Three Four

1.常微分方程概念 由自变量、自变量的未知函数以及函数的导数组成的等式叫做常微分方程。 函 数 的 导 自 变 量 自 变 量 的 函 数

常微分方程概念（续） 微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数,所以一个n阶常微分方程的一般形式如下: 这里 是 的已知函数, 而且一定含有 , y是x的未知函数. 在以下内容中把常微分方程简称为“微分方程”或“方程”.

常微分方程概念（续） n阶线性方程 : 这里 为x的已知函数 。不是线性微分方程的方程称为非线性方程。

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

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

常微分方程的应用（续） Logistic人口模型 ODE 模型: SIR传染病模型 Volterra种群模型

常微分方程的应用(续) ODE 模型: Lorenz族系统 三分子化学动力学模型 一种三神经元神经网络模型

例.几个常见的微分方程模型

4.常微分方程的求解？ 高等数学中，曾经介绍过求一些常微分方程的解，比如：Bernoulli方程和恰当常微分方程(全微分方程)等等。但是对于复杂的方程，求解过程会比较繁琐，如何应用我们已经熟悉的工具——Matlab来求常微分方程的解呢？

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

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

y=-x*log(x)+C1*x 求一右半平面上曲线使其任一点切线在纵轴上截距等于切点的横坐标. 例 解 设曲线方程为 x,y 解 设曲线方程为 依题意曲线上任意一点都应满足: 应用Matlab求解 y=dsolve(‘Dy=y/x-1’,‘x’)，得到 y=-x*log(x)+C1*x

很显然这是一条定义域为正半轴的曲线，当C1取1和-1时图像分别如图所示，可以看出，Matlab能够很方便的求出一些并不直观的常微分方程的解析解.

例 解 高阶 在一个控制系统中，已知控制信号为： 求解如下线性常微分方程： 若已知y(0)=3,Dy(0)=2，求其特解. 首先来计算方程右侧控制项, 键入以下程序: >> syms t; 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).

例 解 当一个系统含有多个依赖于某自变量的函数的时候，就需要用常微分方程组来描述，对一般的线性微分方程组，仍然可以用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给出的任意常数.

例：高等数学中的例题 例2.P278 一电路（电源串联一电阻R和一电感L），电源电动势为E=Emsin wt，电阻R和电感L为常量，求电流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)/(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))得到

例：高等数学中的例题 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 与用变换和特征方程的方法求解相比……

部分非线性方程仍可用dsolve求解 求下述方程组的解析解 >> x=dsolve('Dx=x*(1-x^2)') 例 非线性方程 解 >> 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求解的例子 求下述方程组的解析解（Van der pol 方程） >> syms y mu >> y=dsolve('D2y+mu*(y^2-1)*Dy+y=0') 解 运行后得到： Warning: Compact, analytic solution could not be found. It is recommended that you apply PRETTY to the output. Try mhelp dsolve, mhelp RootOf, mhelp DESol, or mhelp allvalues for more information. ………… 在这里，dsolve命令不能够求出方程的解， Matlab给出了提示.

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

首先把一个高阶微分方程化简为一阶微分方程. 引入新变量替换 m+n维 的 一 阶 微 分 方 程

将Van der pol 方程化为一阶微分方程组 例 将Van der pol 方程化为一阶微分方程组 解 数值解 常微分方程数值解的提法一般是针对一个一阶方程(组)和初始条件来说的， 比如有方程：

求方程数值解的常规方法 有Euler法 变步长Euler法 Runge-Kutta法 Runge-Kutta-Felhberg法等 Adams-Bashforth-Moulton法 主要介绍Euler法和Runge-Kutta法.

Euler方法 1. 基本思想: 在小区间 上用差商来代替方程中导数，Euler方法可以分为向前和向后Euler方法. 如果 中的 t 在 上选取区间的左端点，则有向前Euler公式：

2. 几何意义: 1阶精度 h h的大小影响计算的速 度与精度; 迭代次数的增加会带来 较多的累积误差; 向前Euler公式的 精度并不很高 提高精度的办法: 使用变步长 改进Euler公式等 1阶精度 x 近似解 解析解 h t t0 t1 t2

Runge-Kutta方法 基本思想: Euler方法的基本思想，即差商代替导数，很自然的想法是在区间内多取几个点，将它们的斜率加权平均作为导数近似值，这就是Runge-Kutta方法的思想. 2阶R-K公式:

4阶R-K公式: 公式复杂， 如何编程实现？？？

[t,x]=ode45（或ode23） (Fun,[t_0,t_f],x_0) 在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-文件 在数学摆模型中考虑m=1kg,l=1m,g=9.8m/s^2,初始角度\theta_0=pi/12，初速度为0，阻尼系数\lambda=0.1 根据前面描述无阻尼数学摆的方程可以化为如下方程组： 要求这个方程的数值解首先编写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)为 键入命令： >>theta0=[pi/12,0]; [t,theta]=ode23('fun',[0,10],theta0); %求解，注意初值的选取，theta0一定是2维向量 plot(t,theta(:,1)); %绘出积分曲线图 [t,theta(:,1)] %以矩阵形式输出数值解，第一列为时间，第二列为theta值

可以得到数值解图像 数值解可以列表如下： t theta 0.0000 0.2618 0.0002 0.2618 …… 省略 …… 9.7542 0.1381 9.8337 0.1872 9.9278 0.2304 10.0000 0.2503

考虑阻尼后数学摆的方程可以化为如下方程组： 为求其数值解，首先编写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]; 键入：theta0=[pi/12,0];[t,theta]=ode23('fun2',[0,100],theta0); plot(t,theta(:,1));

可以画出数值解的图像如下： 可以看出，在阻力作用下，数学摆的运动最后将趋于静止.

例：Voterra系统 %……………………vot.m…………………… function 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+e*x(1)+f*x(2))]; [t,x]=ode45('vot',[0,300],x0) plot(x(:,1),x(:,2))

%……………………vot.m…………………… function 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))]; [t,x]=ode45('vot',[0,300],x0) plot(x(:,1),x(:,2)) 与上图相比，系统的性质发生了变化

例 解 画出Lorenz系统 的数值解曲线 首先编写M-文件lorenzeq.m function xdot=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)]; 在命令窗口中运行 >> t_final=100;x0=[0;0;1e-10]; [t,x]=ode45('lorenzeq',[0,t_final],x0); 可以在左侧workspace中看到输出的数值解如下：

续上例

应用plot（t，x）画图得到三条曲线如下图.(分别是t-x1，t-x2，t-x3曲线) 续上例 应用plot（t，x）画图得到三条曲线如下图.(分别是t-x1，t-x2，t-x3曲线)

续上例 也可绘制三维曲线（x1-x2-x3曲线）如下： >> figure;plot3(x(:,1),x(:,2),x(:,3)); 得到了Loren系统具有混沌性质的解曲线. 续上例 应用comet3(x(:,1),x(:,2),x(:,3));可看动态图像，可以看到Lorenz系统具有混沌性的解.

例* 解 针对上述方程，也可在编程时将参数设计在程序中，通过主窗口下赋不同值，求不同参数下微分方程组的解，此时 [ ] 用于站位，不可省略. 首先编写m-文件lorenzleq.m %----------------------------------------- lorenzleq.m--------------------------------------- function xdot=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)]; 在窗口中可以对不同变量赋值， >> 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);

可以看出，变化Lorenz系统的系数，方程的混沌现象消失了，所以该方法可以便于比较一个系统的解是如何随参数变化的. 注意，该方法可以由定义全局变量的方法进行，首先编写M-文件，在M-文件中定义全局变量如下： %----------------------------------------- lorenzlleq.m--------------------------------------- function xdot=lorenzlleq(t,x) global 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)]; 键入程序 >>global beta rho sigma beta=8/3;rho=10;sigma=28;%为参数赋值 x0=[0;0;1e-10]; [t,x]=ode45('lorenzleq',[0,100],x0) plot(t,x) 也可得到相同结果.

向量场的绘制* 例 前面曾提到过，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--------------------------------------- function xdot=ode1(t,x) xdot=1-x^2;

键入： clear c=0.01; x0=-4:.4:4 y0=-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; hold on quiver(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'); hold off

最后可以得到如下图结果，其中有五条积分曲线，分别为过初值（-4，-1. 000001），（-4，-1），（-4，-0 最后可以得到如下图结果，其中有五条积分曲线，分别为过初值（-4，-1.000001），（-4，-1），（-4，-0.8），（-4，1），（-4，4）的微分方程的解，每一条解都是初始值点沿向量场方向运动得到的，具体地说，以y（0）大于-1为初值的解沿向量场走向最终将趋近于1，其他的所有解显然会趋近于负无穷大.

例：Logistic 模型决定的向量场 Logistic人口模型

（续上例）参考程序 clear c=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 on quiver(x,y,u,v);

Simulink6 交互式仿真 应用Matlab 功能强大的仿真工具箱，我们可以从 系统的角度对一个常微分方程（组）进行数值模拟，具体方法介绍如下： 第一步，启动Simulink模块可以在主窗口下输入 simulink 或Simulink3启动或者在Matlab主窗口下从按钮 启动. 启动后可以看到如下界面：

Simulink6 交互式仿真 编写simulink程序与M-file不同，是在仿真编辑窗口中完成的，打开编辑器的方法有： 1.主窗口下 File-New-New model 2.Simulink界面下点击 Create a new model

Simulink6 交互式仿真 对一个常微分方程，通常的办法是将其转化为相应的积分方程来进行模拟，也就是说，对系统 它等价于 积分方程： 下面来看对一个微分方程（组）的仿真计算 对一个常微分方程，通常的办法是将其转化为相应的积分方程来进行模拟，也就是说，对系统 它等价于 积分方程：

2. Simulink中常用的模块简介如下： 积分模块 求导模块 绝对值 积 系数(增益) 求和 导入工作间 自定义函数

3. 各模块都可以自由调整相关参数。 比如积分模块中可以设定初始值 初始值

系数模块中可以设定值 系数值

4.设定仿真环境 仿真时长 误差限

第五步，输出结果（数组（to workspace模块）、绘图等），主要命令（plot,plot3,subplot…）. 第四步，运行(Ctr+T或运行图标) Start 第五步，输出结果（数组（to workspace模块）、绘图等），主要命令（plot,plot3,subplot…）.

例.简单的差分方程 模拟斐波那契数列(第一个月1对小兔子开始，每个月的兔子数量) %rabbit.mdl

例

例 求解Dy=-\alpha*y,y(0)=100 y 仿真时间 Dy 输出

例：求解Van der pol 方程

续上例

续上例 ，方程形式简化为 ，解为

续上例 ，周期解

例.求数学摆的数值解（无阻尼） 其中eta初始值为0，theta初始值为pi/12. Fcn模块对应第二个方程右端形式为-9.8*sin(u)/1. 仿真参数按如下设置：

例.求数学摆的数值解（有阻尼） 与上例相比只需相应加入 一项.

求解*：Dx=-0.1*t*x

作业II 1.给出适当参数绘制出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 vdp1 5.#试着针对Lorenz方程组搭建SImulink仿真模块，完成本文中例题相同的工作； 6.#练习微分方程的求解（可以在高等数学中找方程或自己随意给出）。 7.在练习数学建模基础的同时，查阅文献，找出至少5中人口（或种群）微分方程模型，并求解. ……

Matlab 学习 1.找一本资料 2. help 3. Internet 4. type 谢谢！