理学数学建模matlab.pptx

《理学数学建模matlab.pptx》由会员分享，可在线阅读，更多相关《理学数学建模matlab.pptx（196页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、第1页，共122页常见的问题常见的问题1 分析题目以及选题2 方法的选择3 模型的体现4 对问题求解和软件使用5 论文写作和格式、排版6 其他第1页/共196页第2页，共122页1 分析题目以及选题分析题目以及选题越熟悉的领域越好？X（学经济的就一定要选择*题？）觉得越简单越好？X（*题感觉太难。）感兴趣很重要挖掘内部的数学问题（A题中的数学）抓住主要问题（不要跑题,减速带的设计、不是分析特定类型）要有独到的见解和创新的思路（只要讨论量与量的关系就回归、拟合）要能够根据问题合理安排时间（无法完成题目）第2页/共196页第3页，共122页2 方法的选择方法的选择对问题进行数学描述。（比如A题）要

2、有大致方向，不要直接去查题目的关键词，比如直接搜索“股指预测”等等。平时要积累，比赛时要多查资料多思考。组内讨论。完整学习已有方法，关键步骤要知道为什么（一些特殊的回归模型）。要结合自身条件选择方法，要可求解（偏微）。对以有方法，要能够结合问题特点进行修正（规划问题、最短路等等，东三省D题）。第3页/共196页第4页，共122页3 模型的体现模型的体现要有完整的建模过程，达到让人看懂（看不懂,则白做）细致的问题分析（由此给出建立模型的依据，数学建模绝不简单的是用计算机进行数据分析，比如C题）精确的提炼出所需变量(A题中需要分析的量)问题内在机理（变量之间的关系）选择特定方法的原因（比如一般的统

3、计方法的使用都依赖于强烈的问题的背景和已知的统计学信息）模型的数学表达（便于下文中应用数学技巧处理和求解）必要的解释(补充一些说明使问题叙述更清晰)第4页/共196页第5页，共122页4 模型求解和软件使用模型求解和软件使用模型的求解要能够完整的回答题目中问题。要检验结果的合理性（不一定要体现在论文中）。要有结果分析。常见软件的使用（C，Matlab，Lingo，Mathematica等等）。学会使用可以查到的程序（读懂很重要）。会修改和改正他人的程序（包括改正错误）。结果的表述（不要罗列大量的数据表格）。第5页/共196页第6页，共122页5 论文写作和排版论文写作和排版写些什么内容（说明文

4、？）正文层次段落设置公式/Mathtype（与文字混排要注意行距、标号、对齐）提要式语句和结论性语句图、表（各种软件图的保存、图名、表头、边框）引用的工作要明示、参考文献第6页/共196页第7页，共122页6 其他其他如何查资料队友间合作讨论问题（达成一致意见）良好的写作和编程习惯任务分配（不易完全分块、建模要共同完成）第7页/共196页第8页，共122页Contents1常微分方程求解问题常微分方程求解问题2微分方程微分方程 matlab 求解求解3第8页/共196页第9页，共122页Contents常微分方程的基本概念及其在一些领常微分方程的基本概念及其在一些领域的应用实例域的应用实例;O

5、ne数值解法：数值解法：Euler法法 Runge-Kutta方法；方法；Three求常微分方程解析解的方法以及如何求常微分方程解析解的方法以及如何借助数学软件借助数学软件Matlab来实现来实现;TwoSimulink仿真在求解常微分方程上仿真在求解常微分方程上的应用的应用.FourFive向量向量场绘图场绘图第9页/共196页第10页，共122页1.常微分方程概念常微分方程概念 由自变量、自变量的未知函数以及函数的导数组成由自变量、自变量的未知函数以及函数的导数组成的等式叫做常微分方程。的等式叫做常微分方程。自自变变量量自自变变量量的的函函数数函函数数的的导导数数第10页/共196页第11

6、页，共122页常微分方程概念（续）常微分方程概念（续）微分方程中出现的未知函数最高阶导数的阶微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数数称为微分方程的阶数,所以一个所以一个n阶常微分阶常微分方程的一般形式如下方程的一般形式如下:这里这里 是是 的已知函数的已知函数,而且一而且一定含有定含有 ,y是是x的未知函数的未知函数.在以下内容中在以下内容中把常微分方程简称为把常微分方程简称为“微分方程微分方程”或或“方程方程”.第11页/共196页第12页，共122页常微分方程概念（续）常微分方程概念（续）n阶线性方程阶线性方程:这里这里 为为x的已知函数的已知函数。不是线。不是线性微分

7、方程的方程称为非线性方程。性微分方程的方程称为非线性方程。第12页/共196页第13页，共122页常微分方程概念（续）常微分方程概念（续）常微分方程的解常微分方程的解初值问题（主要介绍本类问题，下面详细展开）初值问题（主要介绍本类问题，下面详细展开）边值问题（与这类问题相关的边值问题（与这类问题相关的Matlab解法可以解法可以参考参考bvp4c函数）函数）第13页/共196页第14页，共122页2.常微分方程的应用常微分方程的应用应用用常微分方程理论是基础数学的重要组成部分之一，在反映客常微分方程理论是基础数学的重要组成部分之一，在反映客观现实世界运动过程的各种量之间的关系中，大量存在满足观

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

9、染病模型传染病模型传染病模型传染病模型VolterraVolterraVolterraVolterra种群模型种群模型种群模型种群模型常微分方程的应用（续）常微分方程的应用（续）常微分方程的应用（续）常微分方程的应用（续）第15页/共196页第16页，共122页LorenzLorenzLorenzLorenz族系统族系统族系统族系统三分子化学动力学模型三分子化学动力学模型三分子化学动力学模型三分子化学动力学模型 一种三神经元神经网络模型一种三神经元神经网络模型一种三神经元神经网络模型一种三神经元神经网络模型常微分方程的应用常微分方程的应用常微分方程的应用常微分方程的应用(续续续续)第16页/共

10、196页第17页，共122页例例.几个常见的微分方程模型几个常见的微分方程模型第17页/共196页第18页，共122页第18页/共196页第19页，共122页第19页/共196页第20页，共122页 高等数学中，曾经介绍过求一些常微分方程的解，比如：BernoulliBernoulli方程方程和恰当常微分方程恰当常微分方程(全微分方程全微分方程)等等。但是对于复杂的方程，求解过程会比较繁琐，如何应用我们已经熟悉的工具Matlab来求常微分方程的解呢？4.常微分方程的求解？常微分方程的求解？第20页/共196页第21页，共122页 对一般的常微分方程（下简称方程），一般来讲是不能够对一般的常微分

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

12、页第22页，共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 第22页/共196页第23页，共122页解解例例求一右半平面上曲线使其任一点切线在纵轴上截距等求一右半平面上曲线使其任一点切线在纵轴上截距等于切点的横坐标于切点的横坐标

13、.设曲线方程为 依题意曲线上任意一点都应满足:应用Matlab求解 y=dsolve(Dy=y/x-1,x)，得到y=-x*log(x)+C1*xxx,y第23页/共196页第24页，共122页很显然这是一条定义域为正半轴的曲线，当很显然这是一条定义域为正半轴的曲线，当C1取取1和和-1时图像时图像分别如图所示，可以看出，分别如图所示，可以看出，Matlab能够很方便的求出一些并能够很方便的求出一些并不直观的常微分方程的解析解不直观的常微分方程的解析解.第24页/共196页第25页，共122页 在一个控制系统中，已知控制信号为：在一个控制系统中，已知控制信号为：求解如下线性常微分方程：求解如下

14、线性常微分方程：例例若已知若已知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*e

15、xp(-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).解解高阶第25页/共196页第26页，共122页解解例例求下述方程组的解析解求下述方程组的解析解 当一个系统含有多个依赖

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

17、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)第26页/共196页第27页，共122页

18、续上例续上例结果很复杂，对其进行化简：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)还可

19、以得到近似的结果：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给出的任意常数.第27页/共196页第28页，共122页例：高等数学中的例题例：高等数学中的例题 例例2.P278一电路（电源串联一电阻R和一电感L），电源电动势为E=Emsin wt，电阻R和电感L为常量，求电流I（t）.解：显然电流满足

20、如下常微分方程：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)得到第28页/共196页第29页，共122页例：高等数学中的例题例：高等数学中的例题 P318求Euler方程的通解求解命令dsolve(x3*D3y+x2*D2

21、y-4*x*Dy=3*x2,x)得到通解1/3*C1*x3-1/2*x2-C2/x+C3与用变换和特征方程的方法求解相比第29页/共196页第30页，共122页解解例例求下述方程组的解析解求下述方程组的解析解 x=dsolve(Dx=x*(1-x2)解得两组通解如下：解得两组通解如下：x=1/(1+exp(-2*t)*C1)(1/2)-1/(1+exp(-2*t)*C1)(1/2)非线性方程非线性方程部分非线性方程仍可用部分非线性方程仍可用部分非线性方程仍可用部分非线性方程仍可用dsolvedsolve求解求解求解求解第30页/共196页第31页，共122页一个不能用一个不能用dsolve求解

22、的例子求解的例子解解例例求下述方程组的解析解（求下述方程组的解析解（Van der pol Van der pol 方程）方程）syms y mu y=dsolve(D2y+mu*(y2-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 informa

23、tion.在这里，在这里，dsolve命令不能够求出方程的解，命令不能够求出方程的解，Matlab给出了提示给出了提示.第31页/共196页第32页，共122页方程无解析解方程无解析解数值解数值解 解的性质解的性质 考察方程随（几）个参考察方程随（几）个参数的变化，以及解是如数的变化，以及解是如何变化的何变化的第二部分：数值解 解的图像解的图像第32页/共196页第33页，共122页引入新变量替换引入新变量替换 首先把一个高阶微分方程化简为一阶微分方程.m+n维维的的一一阶阶微微分分方方程程第33页/共196页第34页，共122页解解例例将将Van der pol 方程化为一阶微分方程组方程化

24、为一阶微分方程组 常微分方程数值解的提法一般是针对一个一阶方常微分方程数值解的提法一般是针对一个一阶方程程(组组)和初始条件来说的，和初始条件来说的，比如有方程：比如有方程：数值解数值解第34页/共196页第35页，共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为自变量取

25、值区间为自变量取值区间;options可给出误差限可给出误差限(缺省时相对误差缺省时相对误差1e-3,绝对误差绝对误差1e-6);具体形式为：具体形式为：options=odeset(reltol,rt,abstol,at).了解了解第35页/共196页第36页，共122页例例.数学摆的求解数学摆的求解在数学摆模型中考虑m=1kg,l=1m,g=9.8m/s2,初始角度theta_0=pi/12，初速度为0，阻尼系数lambda=0.1根据前面描述无阻尼数学摆的方程可以化为如下方程组：要求这个方程的数值解首先编写要求这个方程的数值解首先编写M-文件文件%-fun.m-function xdot

26、=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页/共196页第37页，共122页可以得到数值解图像数值解可以列表如下：t theta 0.0000 0.2618 0.0002 0.261

27、8 省略 9.7542 0.1381 9.8337 0.1872 9.9278 0.230410.0000 0.2503第37页/共196页第38页，共122页键入：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)/

28、l-lambda*theta(2)/m;%-fun2.m-第38页/共196页第39页，共122页可以画出数值解的图像如下可以画出数值解的图像如下：可以看出，在阻力作用下，数学摆的运动最后将趋于静止.第39页/共196页第40页，共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+e*x(1)+f*x(2);%vot.mt,x=ode45(vot,0,300,x0)plot(x(:,1),x(:,2)第40页/共196页第

29、41页，共122页%vot.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.mt,x=ode45(vot,0,300,x0)plot(x(:,1),x(:,2)与上图相比，系统的性质发生了变化第41页/共196页第42页，共122页首先编写首先编写M-文件文件lorenzeq.m%-lorenzeq.m-function xdot=lorenzeq(t,x)xdot=-8/3*x(1)+x(2)*x(3);-10*

30、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页/共196页第43页，共122页续上例续上例第43页/共196页第44页，共122页应用应用plot（t，x）画图得到三条曲线如下图）画图得到三条曲线如下图.(分别是分别是t-x1，t-x2，

31、t-x3曲线曲线)续上例续上例第44页/共196页第45页，共122页也可绘制三维曲线（也可绘制三维曲线（x1-x2-x3曲线）如下：曲线）如下：figure;plot3(x(:,1),x(:,2),x(:,3);得到了得到了Loren系统具有混沌性质的解曲线系统具有混沌性质的解曲线.应用应用comet3(x(:,1),x(:,2),x(:,3);可看动态图像，可以看到可看动态图像，可以看到Lorenz系统具有系统具有混沌性的解混沌性的解.续上例续上例第45页/共196页第46页，共122页解解例例*针对上述方程，也可在编程时将参数设计在程序中，通过针对上述方程，也可在编程时将参数设计在程序中

32、，通过主窗口下赋不同值，求不同参数下微分方程组的解，此时主窗口下赋不同值，求不同参数下微分方程组的解，此时 用于站位，不可省略用于站位，不可省略.首先编写首先编写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

33、=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);subpl

34、ot(1,3,2);plot(t,x(:,2);subplot(1,3,3);plot(t,x(:,3);第46页/共196页第47页，共122页第47页/共196页第48页，共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*

35、x(3);-x(1)*x(2)+sigma*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页/共196页第49页，共122页隐式微分方程求解隐式微分方程求解 ode15iF【t，x（t），Dx（t）】=0 x(0)=x0;Dx(0)=dx0例：x2+Dx2=1;x0=0解:x0=0dx0=1x0F=1%保留x0dx0F=0%重新计算dx0f=(t,x,dx

36、)x2+dx2-1;x0,dx0=decic(f,0,x0,x0F,dx0,dx0F)%0为t0t,x=ode15i(f,0,1,x0,dx0)plot(t,x)第49页/共196页第50页，共122页微分代数方程的求解微分代数方程的求解方程组中包含微分方程和一般方程例：dx1=-0.2*x1+x2*x3+0.3*x1*x2;dx2=2*x1*x2-5*x2*x3-2*x22 0=x1+x2+x3-1%(x1+x2+x3=1)解：f=(t,x)-0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2);2*x(1)*x(2)-5*x(2)*x(3)-2*x(2)2;x(1)+x(2)+

37、x(3)-1M=1 0 0;0 1 0;0 0 0options=odeset;options.Mass=Mx0=0.8 0.1 0.1t,x=ode15s(f,0,20,x0,options);plot(t,x)第50页/共196页第51页，共122页边值问题：help bvp4c偏微分方程：help pdepe OR pdetool自学第51页/共196页第52页，共122页向量场向量场第52页/共196页第53页，共122页向量场的绘制*前面曾提到过，Matlab中提供了命令quiver来绘制向量场，具体格式为quiver（x,y,u,v），x和y表示向量起点，u和v确定向量的方向。例例

38、绘出 所确定的向量场，以及过（-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-x2;%-ode.m-第53页/共196页第54页，共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);%绘制向量场

39、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第54页/共196页第55页，共122页最后可以得到如下图结果，其中有五条积分曲线，分别为过初值（-4，-1.000001），（-4，-1），（-4，-0.8），（-4，1），（-4，4）的微分方程的

40、解，每一条解都是初始值点沿向量场方向运动得到的，具体地说，以y（0）大于-1为初值的解沿向量场走向最终将趋近于1，其他的所有解显然会趋近于负无穷大.第55页/共196页第56页，共122页例：例：Logistic 模型决定的向量场模型决定的向量场LogisticLogistic人口模型人口模型人口模型人口模型 第56页/共196页第57页，共122页（续上例）参考程序（续上例）参考程序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一同

41、决定向量场方向v=c*(0.1*(1-y/50).*y)./d;hold onquiver(x,y,u,v);第57页/共196页第58页，共122页Simulink6 交互式仿真 应用应用Matlab 功能强大的仿真工具箱，我们可以从功能强大的仿真工具箱，我们可以从 系统的角度对一个常系统的角度对一个常微分方程（组）进行数值模拟，具体方法介绍如下微分方程（组）进行数值模拟，具体方法介绍如下：第一步，启动第一步，启动Simulink模块模块可以在主窗口下输入可以在主窗口下输入 simulink 或或Simulink3启动或者在启动或者在Matlab主窗口下从按钮主窗口下从按钮 启动启动.启动后

42、可以看到如下界面：启动后可以看到如下界面：第58页/共196页第59页，共122页Simulink6 交互式仿真编写编写simulink程序与程序与M-file不同，是在仿真编辑窗口中完成的，打开编辑器的方法有：不同，是在仿真编辑窗口中完成的，打开编辑器的方法有：1.主窗口下主窗口下 File-New-New model2.Simulink界面下点击界面下点击 Create a new model第59页/共196页第60页，共122页1.对一个常微分方程，通常的办法是将其转化为相应对一个常微分方程，通常的办法是将其转化为相应的积分方程来进行模拟，也就是说，对系统的积分方程来进行模拟，也就是说

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

44、64页，共122页4.4.设定仿真环境设定仿真环境 仿真时长仿真时长误差限误差限第64页/共196页第65页，共122页Start 第四步，运行第四步，运行(Ctr+T或运行图标或运行图标)第五步，输出结果（数组（第五步，输出结果（数组（to workspace模块）、模块）、绘图等），主要命令（绘图等），主要命令（plot,plot3,subplot）.第65页/共196页第66页，共122页例例.简单的差分方程简单的差分方程x_n=x_n-1+x_n-2模拟斐波那契数列(第一个月1对小兔子开始，每个月的兔子数量)%rabbit.mdl第66页/共196页第67页，共122页例例第67页/共

45、196页第68页，共122页例例 求解求解Dy=-alpha*y,y(0)=100yDy输出仿真时间第68页/共196页第69页，共122页第69页/共196页第70页，共122页切换切换ODE求解求解例：假设数学摆运动的平面内阻尼分为两部分Theta0时，lambda=0.1Theta0时，lambda=0.5第70页/共196页第71页，共122页脉冲方程脉冲方程 Dx=-.5x+pulse周期方波第71页/共196页第72页，共122页随机参数微分方程随机参数微分方程Dx=-.5x+X X随机数随机数注：其他分布随机数可以用统计方法得到第72页/共196页第73页，共122页例：求解例：

46、求解Van der pol 方程方程 第73页/共196页第74页，共122页续上例续上例第74页/共196页第75页，共122页 ，方程形式简化为，方程形式简化为 ，解为，解为 续上例续上例第75页/共196页第76页，共122页，周期解，周期解续上例续上例第76页/共196页第77页，共122页例例.求数学摆的数值解（无阻尼）求数学摆的数值解（无阻尼）其中eta初始值为0，theta初始值为pi/12.Fcn模块对应第二个方程右端形式为-9.8*sin(u)/1.仿真参数按如下设置：第77页/共196页第78页，共122页第78页/共196页第79页，共122页例例.求数学摆的数值解（有阻

47、尼）求数学摆的数值解（有阻尼）与上例相比只需相应加入 一项.第79页/共196页第80页，共122页第80页/共196页第81页，共122页求解求解*：Dx=-0.1*t*x第81页/共196页第82页，共122页时滞微分方程时滞微分方程自变量带延迟（例如：人口问题信号传输传染病）Dx=-5x(t)+5x(t-1)第82页/共196页第83页，共122页Dx=-5x(t)+5x(t-0.2*sin(x)第83页/共196页第84页，共122页作业作业II1.给出适当参数绘制出SIR模型的解并给出适当的解释;2.#熟悉ODE的解析解与数值解的Matlab求解方法；3.了解一维方程向量场绘图；自己

48、给出适当方程，绘图、观察积分曲线的走向；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中人

49、口（或种群）微分方程模型，并求解.第84页/共196页第85页，共122页哥尼斯堡七桥问题哥尼斯堡七桥问题 11.1 图论的基本概念图论的基本概念哥尼斯堡是一座城市哥尼斯堡是一座城市,位于普莱格尔河之上位于普莱格尔河之上,河中有两个岛屿河中有两个岛屿A与与D,为了沟通城市交通为了沟通城市交通,在岛与岛及岛与岸在岛与岛及岛与岸之间架设了之间架设了7座桥。座桥。的前提下的前提下,最后又返回到原出发地？最后又返回到原出发地？要从任何一地出发要从任何一地出发,能否在每座桥只允许通过一次能否在每座桥只允许通过一次哥尼斯堡哥尼斯堡(Konigsberg)七桥七桥 问：问：第85页/共196页第86页，共1

50、22页每块陆地每块陆地图中的点图中的点A,B,C,D:连接相应点之间的线连接相应点之间的线:连接陆地之间的连接陆地之间的7座桥座桥哥尼斯堡七桥问题哥尼斯堡七桥问题哥尼斯堡七桥问题哥尼斯堡七桥问题 过每一条边且仅能划一次，过每一条边且仅能划一次，从任意一点出发从任意一点出发,能否用笔划能否用笔划最后又回到原出发点？最后又回到原出发点？哥尼斯堡七桥问题不可能有解哥尼斯堡七桥问题不可能有解哥尼斯堡七桥问题不可能有解哥尼斯堡七桥问题不可能有解 与每个点相与每个点相连的边应为连的边应为偶数条偶数条每一次通每一次通过点的边过点的边总是两条总是两条进入和离进入和离开开哥尼斯堡七桥问题哥尼斯堡七桥问题 第86