1 / 76

ODE 的求解

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

britain
Download Presentation

ODE 的求解

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. ODE的求解 MCM 2012.10.29 Edit by niuben

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

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

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

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

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

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

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

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

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

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

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

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

  14. 解析解 求解方程 例 应用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 解

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

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

  17. 高阶 在一个控制系统中，已知控制信号为： 例 求解如下线性常微分方程： 若已知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).

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

  19. 续上例 • 结果很复杂，对其进行化简： • >>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给出的任意常数.

  20. 例：高等数学中的例题 例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))得到

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

  22. 部分非线性方程仍可用dsolve求解 求下述方程组的解析解 例 非线性方程 解 >> x=dsolve('Dx=x*(1-x^2)') 解得两组通解如下： x = [ 1/(1+exp(-2*t)*C1)^(1/2)] [ -1/(1+exp(-2*t)*C1)^(1/2)]

  23. 一个不能用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给出了提示.

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

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

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

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

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

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

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

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

  32. 在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). 了解

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

  34. 可以得到数值解图像 数值解可以列表如下： 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

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

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

  37. 例：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))]; • %……………………vot.m…………………… • [t,x]=ode45('vot',[0,300],x0) • plot(x(:,1),x(:,2))

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

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

  40. 续上例

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

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

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

  44. 可以看出，变化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)]; %----------------------------------------- lorenzlleq.m--------------------------------------- 键入程序 >>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) 也可得到相同结果.

  45. 向量场的绘制* 前面曾提到过，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; %-----------------------------------------ode.m---------------------------------------

  46. 键入： 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

More Related