《MATALAB微分方程》PPT课件.ppt

《《MATALAB微分方程》PPT课件.ppt》由会员分享，可在线阅读，更多相关《《MATALAB微分方程》PPT课件.ppt（35页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、数学建模数学建模 微微 分分 方方 程程 在研究实际问题时，常常会联系到某些变量的变化在研究实际问题时，常常会联系到某些变量的变化率或导数，率或导数，这样所得到变量之间的关系式就是微分这样所得到变量之间的关系式就是微分方程模型。微分方程模型反映的是变量之间的间接方程模型。微分方程模型反映的是变量之间的间接关系，因此，要得到直接关系，就得求微分方程。关系，因此，要得到直接关系，就得求微分方程。求解微分方程有三种方法：求解微分方程有三种方法：1 1）求精确解；）求精确解；2 2）求数值解（近似解）；）求数值解（近似解）；3 3）定性）定性理论方法。理论方法。一、导弹追踪问题一、导弹追踪问题 设位于

2、坐标原点的甲舰向位于设位于坐标原点的甲舰向位于x轴上点轴上点A(1,0)处的乙舰处的乙舰发射导弹，导弹头始终对准乙舰发射导弹，导弹头始终对准乙舰.如果乙舰以最大的速度如果乙舰以最大的速度v0(是常数是常数)沿平行于沿平行于y轴的直线行驶，导弹的速度是轴的直线行驶，导弹的速度是5v0，求导求导弹运行的曲线方程弹运行的曲线方程.又乙舰行驶多远时，导弹将它击中又乙舰行驶多远时，导弹将它击中？解法一（解析法）由(1),(2)消去t整理得模型:二二 范范.梅格伦（梅格伦（Van Meegren）伪造名画案伪造名画案 第二次世界大战比利时解放后，荷兰保安机关开始搜捕纳粹分子的合作者，发现一名三流画家H.A

3、.Vanmeegren曾将17世纪荷兰著名画家Jan.Vermeer的一批名贵油画盗卖给德寇，于1945年5月29日通敌罪逮捕了此人。Vanmeegren被捕后宣称他从未出卖过荷兰的利益，所有的油画都是自己伪造的，为了证实这一切，在狱中开始伪造Vermeer的画耶稣在学者中间。当他的工作快完成时，又获悉他可能以伪造罪被判刑，于是拒绝将画老化，以免留下罪证。为了审理这一案件，法庭组织了一个由化学家、物理学家、艺术史学家等参加的国际专门小组，采用了当时最先进的科学方法，动用了X-光线透视等，对颜料成份进行分析，终于在几幅画中发现了现代物质诸如现代颜料钴蓝的痕迹。这样，伪造罪成立，Vanmeegre

4、n被判一年徒刑。1947年11月30日他在狱中心脏病发作而死去。但是，许多人还是不相信其余的名画是伪造的，因为，Vanmeegren在狱中作的画实在是质量太差，所找理由都不能使怀疑者满意。直到20年后，1967年，卡内基梅隆大学的科学家们用微分方程模型解决了这一问题。原理原理著名物理学家卢瑟夫（Rutherford）指出：物质的放射性正比于现存物质的原子数。设 时刻的原子数为 ，则有为物质的衰变常数。初始条件半衰期能测出或算出，只要知道 就可算出这正是问题的难处，下面是间接确定 的方法。断代。油画中的放射性物质油画中的放射性物质 白铅（铅的氧化物）是油画中的颜料之一，应用已有2000余年，白铅

5、中含有少量的铅(Pb210)和更少量的镭(Ra226)。白铅是由铅金属产生的，而铅金属是经过熔炼从铅矿中提取来出的。当白铅从处于放射性平衡状态的矿中提取出来时，Pb210的绝大多数来源被切断，因而要迅速蜕变，直到Pb210与少量的镭再度处于放射平衡，这时Pb210的蜕变正好等于镭蜕变所补足的为止。铀238镭226铅210钋210铅206（放射性）（无放射性）假设假设（1）镭的半衰期为1600年，我们只对17 世纪的油画感兴趣，时经300多年，白铅中镭至少还有原量的90%以上，所以每克白铅中每分钟镭的衰变数可视为常数，用 表示。（2）钋的半衰期为138天容易测定，铅210的半衰期为22年，对要鉴

6、别的300多年的颜料来说，每克白铅中每分钟钋的衰变数与铅210的衰变数可视为相等。建模建模设 时刻每克白铅中含铅210的数量为 ，为制造时刻 每克白铅中含铅210的数量。为铅210的衰变常数。则油画中铅210含量求解求解均可测出。可算出白铅中铅的衰变率 ，再于当时的矿物比较，以鉴别真伪。矿石中铀的最大含量可能 23%，若白铅中铅210每分钟衰变超过3 万个原子，则矿石中含铀量超过 4%。测定结果与分析测定结果与分析画名画名钋钋210衰变原子数衰变原子数镭镭226衰变原子数衰变原子数Emmaus的信徒们8.50.82洗足12.60.26读乐谱的妇人10.30.3弹曼陀林的妇人8.20.17做花边

7、的人1.51.4欢笑的女孩5.26.0若第一幅画是真品，铅210每分钟每克衰变不合理，为赝品。同理可检验第2，3，4幅画亦为赝品，而后两幅画为真品。微分方程数值解微分方程数值解 1欧拉方法欧拉方法欧拉方法的基本原理欧拉方法的基本原理:消除导数项（离散化）。消除导数项（离散化）。欧拉法的一般步骤欧拉法的一般步骤:欧拉方法的特点欧拉方法的特点:易于理解，计算量小，精度低。易于理解，计算量小，精度低。2梯形法梯形法梯形法的一般步骤梯形法的一般步骤:梯形法的特点梯形法的特点:，计算量大，精度高。，计算量大，精度高。龙格库塔法龙格库塔法 龙格库塔法的基本思想：龙格库塔法的基本思想：由微分中值定理由微分中

8、值定理：微分方程的解析解微分方程的解析解 求单个微分方程的解析解命令求单个微分方程的解析解命令:dsolve(方程方程1,初始条件初始条件,自变量自变量)结 果：u=tg(t-c)求微分方程组的解析解求微分方程组的解析解:x1,x2,xndsolve(方程方程1,方程方程2,方程方程n,初始条件初始条件,自变量自变量)例例 2求微分方程求微分方程 的通解，并验证。的通解，并验证。y=dsolve(Dy+2*x*y=x*exp(-x2),x)syms x;diff(y)+2*x*y-x*exp(-x2)例例 3求微分方程求微分方程 在初值条件在初值条件 下的特解，并画出解函数的图形。下的特解，并

9、画出解函数的图形。y=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1),x);ezplot(y)例例4 求微分方程组求微分方程组 在初值条件在初值条件 下的特解，并画出解函数的图形。下的特解，并画出解函数的图形。x,y=dsolve(Dx+5*x+y=exp(t),Dy-x-3*y=0,x(0)=1,y(0)=0,t)ezplot(x,y,0,1.3)微分方程的数值解微分方程的数值解常微分方程数值解的定义常微分方程数值解的定义 在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解。而在实际上对初值问题，一般是要求得到解在若干个点上满足规定精确度的近似值，或者得到一

10、个满足精确度要求的便于计算的表达式。因此，研究常微分方程的数值解法是十分必要的因此，研究常微分方程的数值解法是十分必要的。用用Matlab软件求常微分方程的数值解软件求常微分方程的数值解t，y=solver（f,ts,y0）ode45 ode23 ode113ode15sode23s由待解方程写成的m-文件名ts=t0，tf，t0、tf为自变量的初值和终值函数的初值ode23：组合的2/3阶龙格-库塔-芬尔格算法ode45：运用组合的4/5阶龙格-库塔-芬尔格算法自变量值函数值Matlab提供的ODE求解器求解器 ODE类型特点说明ode45非刚性非刚性单步法；单步法；4，5 阶阶 R-K 方

11、法；方法；累计截断误差为累计截断误差为(x)3大部分场合的大部分场合的首选方法首选方法ode23非刚性非刚性单步法；单步法；2，3 阶阶 R-K 方法；方法；累计截断误差为累计截断误差为(x)3使用于精度较低的情形使用于精度较低的情形ode113非刚性非刚性多步法；多步法；Adams算法；高低精算法；高低精度均可到度均可到 10-310-6计算时间比计算时间比 ode45 短短ode23t适度刚性适度刚性 采用梯形算法采用梯形算法适度刚性情形适度刚性情形ode15s刚性刚性多步法；多步法；Gears 反向数值微反向数值微分；精度中等分；精度中等若若 ode45 失效时，可失效时，可尝试使用尝试

12、使用ode23s刚性刚性单步法；单步法；2 阶阶Rosebrock 算算法；低精度法；低精度当精度较低时，计算时当精度较低时，计算时间比间比 ode15s 短短ode23tb刚性刚性梯形算法；低精度梯形算法；低精度当精度较低时，计算时当精度较低时，计算时间比间比ode15s短短 1、在解n个未知函数的方程组时，x0和x均为n维向量，m-文件中的待解方程组应以x的分量形式写成.2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.注意注意:x,y=ode23(fun,0,0.5,1)fun=inline(-2*y+2*x2+2*x,x,y);x,y=ode23(fun,

13、0,0.5,1)求初值问题求初值问题 的数值解，求解范围的数值解，求解范围为为 0,0.5例例 1图中，图中，y1的图形为实线，的图形为实线，y2的图形为的图形为“*”线，线，y3的图形为的图形为“+”线线.function f=cxd1(t,y)f=y(2)*y(3);-y(1)*y(3);-0.51*y(1)*y(2)T,Y=ode45(cxd1,0 1 2,0 1 1)plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+)解解 令令 y1=x，y2=y11、建立M文件function f=cxd(t,y)f=y(2);1000*(1-y(1)2)*y(2)-y(1)2、取、取t0=0，tf=3000，输入命令：输入命令：T,Y=ode15s(cxd,0 3000,2 0);plot(T,Y(:,1),-)