书签 分享 收藏 举报 版权申诉 / 50
上传文档赚钱

类型Matlab在常微分方程求解中应用课件.ppt

  • 上传人(卖家):ziliao2023
  • 文档编号:5961708
  • 上传时间:2023-05-19
  • 格式:PPT
  • 页数:50
  • 大小:367.50KB
  • 【下载声明】
    1. 本站全部试题类文档,若标题没写含答案,则无答案;标题注明含答案的文档,主观题也可能无答案。请谨慎下单,一旦售出,不予退换。
    2. 本站全部PPT文档均不含视频和音频,PPT中出现的音频或视频标识(或文字)仅表示流程,实际无音频或视频文件。请谨慎下单,一旦售出,不予退换。
    3. 本页资料《Matlab在常微分方程求解中应用课件.ppt》由用户(ziliao2023)主动上传,其收益全归该用户。163文库仅提供信息存储空间,仅对该用户上传内容的表现方式做保护处理,对上传内容本身不做任何修改或编辑。 若此文所含内容侵犯了您的版权或隐私,请立即通知163文库(点击联系客服),我们立即给予删除!
    4. 请根据预览情况,自愿下载本文。本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
    5. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007及以上版本和PDF阅读器,压缩文件请下载最新的WinRAR软件解压。
    配套讲稿:

    如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。

    特殊限制:

    部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。

    关 键  词:
    Matlab 微分方程 求解 应用 课件
    资源描述:

    1、Matlab在常微分方程在常微分方程求解中的应用求解中的应用实验目的实验目的 (1 1)学会用学会用MatlabMatlab软件求解微分方程的初值问题软件求解微分方程的初值问题 (2 2)了解微分方程数值解思想,掌握基本的微分方程数值解方法)了解微分方程数值解思想,掌握基本的微分方程数值解方法 (3 3)学会根据实际问题建立简单微分方程数学模型)学会根据实际问题建立简单微分方程数学模型 (4 4)了解计算机数据仿真、数据模拟的基本方法)了解计算机数据仿真、数据模拟的基本方法 n17世纪:初等解法世纪:初等解法n18世纪:初等解法和无穷级数方法世纪:初等解法和无穷级数方法n19世纪:解的存在性、

    2、奇点理论、定性世纪:解的存在性、奇点理论、定性理论、稳定性理论理论、稳定性理论l 包含一个自变量和它的未知函数以及未知函数的导数的包含一个自变量和它的未知函数以及未知函数的导数的等式等式l 形成和发展与力学、天文学、物理学及其他自然科学技形成和发展与力学、天文学、物理学及其他自然科学技术的发展互相促进和推动术的发展互相促进和推动 常微分方程常微分方程yxf,ybxaD,:1,yx2,yx2121,yyLyxfyxf xy定理 设函数设函数在区域在区域上连续,且在区域上连续,且在区域D内满足李普希兹内满足李普希兹(Lipschitz)条件,即存在条件,即存在正数正数L,使得对于,使得对于R内任意

    3、两点内任意两点与与,恒有恒有则初值问题则初值问题(1)的解的解存在并且唯一。存在并且唯一。常微分方程的解析解常微分方程的解析解 求微分方程(组)的解析解命令:dsolve(方程方程1,方程方程2,方程方程n,初始条件初始条件,自变量自变量)结 果:u=tg(t-c)记号:在表达微分方程时,用字母D表示求微分,D2、D3等表示求高阶微分。D后所跟字母为因变量,自变量可以指定或由系统规则选定为缺省。例如:微分方程 可以表示为D2y=0.022dxyd 解解 输入命令:y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x)结 果 为:y=3e-2xsin(5x)例例

    4、3 求微分方程组的通解.zyxdtdzzyxdtdyzyxdtdx244354332解解 输入命令:x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z,t);x=simple(x)%将将x化简化简 y=simple(y)z=simple(z)结 果 为:x=(c1-c2+c3+c2e-3t-c3e-3t)e2t y=(c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t z=(-c1e-4t+c2e-4t+c1-c2+c3)e2t 虽然说解析解是最精确的,但是实际问虽然说解析解是最精确的,但是实际问题中

    5、常要求研究常微分方程的数值解题中常要求研究常微分方程的数值解常微分方程的数值解常微分方程的数值解 常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)。另外,有些初值问题虽然有初等解,但由于形式太复杂不便于应用。因此,有必要探讨常微分方程初值问题的数值解法。以下主要介绍一阶常微分方程初值问题几种经典数值解方法:欧拉法及改进的欧拉法。其它方法:龙格-库塔法、阿达姆斯方法;一阶微分方程组与高阶方程初值问题的数值解法;二阶常微分方程值问题的差分方法等。求解常微分方程初值问题的数值解的整体思路:00)(),(yxyyxfdxdy (1))(xynxxxx210,210nyyyy ny1iii

    6、xxh寻求准确解在一系列离散节点:上的近似值 称为问题的数值解,数值解所满足的离散方程统称为差称为步长,实用中常取定步长。分格式,建立数值解法的一些途径建立数值解法的一些途径001i)y(xy)f(x,y ,1,2,1,0 ,xynihxi解微分方程:可用以下离散化方法求设1、用差商代替导数、用差商代替导数 若步长h较小,则有hxyhxyxy)()()(故有公式:1-n,0,1,2,i )(),(001xyyyxhfyyiiii此即欧拉法欧拉法Euler法的几何意义:2、使用数值积分、使用数值积分对方程y=f(x,y),两边由xi到xi+1积分,并利用梯形公式,有:)(,()(,(2)(,()

    7、()(11111iiiiiixxiixyxfxyxfxxdttytfxyxyii实际应用时,与欧拉公式结合使用:,2,1,0 ),(),(2),()(11)1(1)0(1kyxfyxfhyyyxhfyykiiiiikiiiii的计算。然后继续下一步,取时,当满足,对于已给的精确度)(y y 2i111i)(1)1(1kikikiyyy此即改进的欧拉法改进的欧拉法故有公式:)(),(),(200111xyyyxfyxfhyyiiiiii使得中值定理可知,存在可微的条件下,由微分在函数,10)(xy)()()(1hxyhxyxynnn上的平均斜率在区间为函数称,)()()(11nnnnxxxyhx

    8、yxyhxyxyxxyxfxyxfxxxxyxnnnnnnnnnnnnn)()(,),(),(,)(111-111上的平均斜率区间的平均值代替了曲线在处的预测斜率和在处的斜率线在改进的欧拉公式:用曲上的平均斜率在区间了函数的斜率(导数值)代替欧拉公式:用曲线在的精度更高的计算公式计算由此可构造出由上的平均斜率,在区间的加权平均值作为曲线后将这些点处的斜率值值,然内多预测几个点的斜率间的公式;如果设法在区一种计算算法,就可以得到上的平均斜率提供一种间启发:只要对曲线在区11111,nnnnnnnnnyyxxxxyxx3、使用泰勒公式、使用泰勒公式 以此方法为基础,有龙格龙格-库塔法库塔法、线性多

    9、步法线性多步法等方法。4、数值公式的精度、数值公式的精度 当一个数值公式的截断误差可表示为O(hk+1)时(k为正整数,h为步长),称它是一个k阶公式阶公式。k越大,则数值公式的精度越高。欧拉法是一阶公式,改进的欧拉法是二阶公式。龙格-库塔法有二阶公式和四阶公式。线性多步法有四阶阿达姆斯外插公式和内插公式。ODE 指令列表 nMATLAB 用於求解常微分方程式的指令:指令方法应用ODE类别ode45Explicit Runge-Kutta(4,5)pair of Dormand-PrinceNonstiff ODEode23Explicit Runge-Kutta(2,3)pair of Bo

    10、gacki and ShampineNonstiff ODEode113Variable order Adams-Bashforth-Moulton PECE solverNonstiff ODEode15sNumerical differentiation formulas(NDFS)Stiff ODEode23sModified Rosenbrock formula of order 2Stiff ODEode23t Trapezoidal rule with a“free”interpolantStiff ODEode23tbImplicit Runge-Kutta formula wi

    11、th a backward differentiation formula of order twoStiff ODE 适用于 Nonstiff 系統 n速率(即微分值)差异相常大 n使用一般的ode45、ode23或ode113來求解,可能会使得积分的步长(Step Sizes)变得很小,以便降低积分误差至可容忍范围以內,会导致计算时间过长 n专门对付Stiff系统的指令,例如ode15s、ode23s、ode23t及ode23tb 适用于 Nonstiff 系統 一般的常微分方程式都是 Nonstiff 系統 直接采用 ode45、ode23 或 ode113 來求解 提示 n使用 Sim

    12、ulink 來求解常微分方程式Simulink是和MATLAB共同使用的一套软件 可使用拖拉的方式來建立动力系统 可直接产生C语言代码或进行动画演示 功能非常强大 n使用ODE指令时,必须先将要求解的ODE表示 成一个函数 输入为t(时间)及y(状态函数:State Variables)输出则为dy(状态变量的微分值)n若ODE函数的文件为odeFile.m,则调用ODE指令 的格式如下:t,y=solver(odeFile,t0,t1,y0)t0,t1 是积分的时间区间 y0代表起始条件(Initial Conditions)solver是前表所列的各种ODE指令 t是输出的时间向量 y是对

    13、应的状态变量向量 t,y=solver(f,ts,y0,options)ode45 ode23 ode113ode15sode23s由待解方程写成的m-文件名ts=t0,tf,t0、tf为自变量的初值和终值函数的初值ode23:组合的2/3阶龙格-库塔-芬尔格算法ode45:运用组合的4/5阶龙格-库塔-芬尔格算法自变量值函数值用于设定误差限(缺省时设定相对误差10-3,绝对误差10-6),命令为:options=odeset(reltol,rt,abstol,at),rt,at:分别为设定的相对误差和绝对误差.1、在解n个未知函数的方程组时,y0和y均为n维向量,m-文件中的待解方程组应以y

    14、的分量形式写成.2、使用Matlab软件求数值解时,高阶微分方程必须等价地变换成一阶微分方程组.注意注意:Van der Pol微分方程n1928年荷兰的范德波耳(Van der Pol)为描述LC回路的电子管振荡器建立了著名的vanderPol方程.它在自激振荡理论中有着重要的意义,一直作为数学物理方程中的一个基本方程.这是一个具有可变非线性阻尼的微分方程,代表了一类极为典型的非线性问题.和其他非线性微分方程在数学上无法精确求解一样,人们一直在努力寻找求解这类方程近似解析解的方法,并乐于用Van der Pol方程来检验求解方法的有效性.Van der Pol微分方程其方程形式为:n令 0)

    15、1(2yyyy121yyyy将Van der Pol微分方程化成标准形式 1221221)1(yyyyyy),(tyFyy写成向量的形式:为一个向量,代表状态变量n设 =1,ODE文件(vdp1.m)显示如下:type vdp1.m type vdp1.m function dy=vdp1(t,y)mu=1;dy=y(2);mu*(1-y(1)2)*y(2)-y(1);建立了vdp1.m后,即可选用前述ODE指令来求解 n在 =1时,van der Pol方程并非 Stiff系统,所以使用ode45来显示结果 odeBasic01.mode45(vdp1,0 25,3 3);n 0 25代表积

    16、分的时间区间,3 3 则代表起始条件n因为没有输出变量,所以上述程序执行结束后,MATLAB 只画出状态变量对时间的图形0510152025-3-2-101234 odeGetData01.mt,y=ode45(vdp1,0 25,3 3);plot(t,y(:,1),t,y(:,2),:);xlabel(Time t);ylabel(Solution y(t)and y(t);legend(y(t),y(t);0510152025-3-2-101234Time tSolution y(t)and y(t)y(t)y(t)n也可以画出 及 在相位平面(Phase Plane)的运动情况 ode

    17、PhastPlot01.m)(ty)(tyt,y=ode45(vdp1,0 25,3 3);plot(y(:,1),y(:,2),-o);xlabel(y(t);ylabel(y(t);-3-2-101234-3-2-10123y(t)y(t)n当 值越来越大时,Van der Pol方程就渐变为一个Stiff系统,此时若要求解该类问题,就必须用对应于Stiff系统的指令 n将 值改成 1000,ODE 文件改写成(vdp2.m):type vdp2.m function dy=vdp2(t,y)mu=1000;dy=y(2);mu*(1-y(1)2)*y(2)-y(1);n选用专门用来解决

    18、Stiff 系统问题的 ODE 指令,例如 ode15s,來求解此系统并作图显示 ode15s01.mt,y=ode15s(vdp2,0 3000,2 1);subplot(2,1,1);plot(t,y(:,1),-o);xlabel(Time t);ylabel(y(t);subplot(2,1,2);plot(t,y(:,2),-o);xlabel(Time t);ylabel(y(t);n由上图可知,的变化相当剧烈(超过 ),这也就是 Stiff 系统的特点 050010001500200025003000-4-2024Time ty(t)05001000150020002500300

    19、0-2000-1000010002000Time ty(t)(ty1000n绘制二维平面相位图:ode15s02.mn若要产生在某些特定时间点的状态变量值,则调用 ODE 指令的格式可改成:t,y=solver(odeFile,t0,t1,tn,y0)其中t0,t1,tn即是特定时间点所形成的向量 t,y=ode15s(vdp2,0 3000,2 1);subplot(1,1,1);plot(y(:,1),y(:,2),-o);xlabel(y(t);ylabel(y(t)-2.5-2-1.5-1-0.500.511.522.5-1500-1000-500050010001500y(t)y(t

    20、)ODE 指令的选项nODE 指令可以接受第四个输入变量,代表积分过程用到的各种选项(Options),此种ODE指令的格式为:t,y=solver(odeFile,t0,tn,y0,options)其中 options 是由 odeset 指令來控制的结构变量 结构变量即包含了积分过程用到的各种选项 odeset的一般格式如下:options=odeset(name1,value1,name2,value2,)nname1的值为value1,name2的值为value2 n也可以只改变一个现有的options中某栏位的值,格式如下:newOptions=odeset(options,name

    21、,value);n若要获取某栏位的值,可用 odeget,格式如下:value=odeget(otpions,name);n当使用 odeset 指令時,若无输入任何变量名,则 odeset 会显示所有变量名及其值,并以大括号代表预设值 odeset AbsTol:positive scalar or vector 1e-6 RelTol:positive scalar 1e-3 NormControl:on|off NonNegative:vector of integers OutputFcn:function_handle OutputSel:vector of integers Ref

    22、ine:positive integer Stats:on|off InitialStep:positive scalar MaxStep:positive scalar BDF:on|off MaxOrder:1|2|3|4|5 Jacobian:matrix|function_handle JPattern:sparse matrix Vectorized:on|off Mass:matrix|function_handle MStateDependence:none|weak|strong MvPattern:sparse matrix MassSingular:yes|no|maybe

    23、 InitialSlope:vector Events:function_handle 由 odeset 产生的 ODE 选项 310610类别栏位名称资料形态预设值说明誤差容忍度之相关栏位RelTol正数量相对误差容忍度AbsTo1正数量或向量绝对误差容忍度积分输出之相关栏位OutPutFcn字串odeplot输出函式(若 ODE 指令无输出变量,則在数值积分执行完毕后,MATLAB会呼叫此输出函数)OutputSel索引向量全部ODE 指令之输出变量的索引值,以決定那些输出变量之元素将被送到输出函式Refine正整数1或4(for ode45)Refine=2 可产生两倍数量的输出点,Re

    24、fine=3 可产生三倍数量的输出点,依此类推Statson 或 offoffStats=on 产生计算过程的各种统计资料常用到的档位说明 nf在积分误差容忍度方面,每一次积分所产生的局部误差e(i),必须满足下列方程式:max(RelTol*,AbsTol(i)其中i代表第i个状态 降低RelTol及AbsTol来求得更精确的积分结果)(ie)(iyodeRelTol01.msubplot(2,1,1);ode45(vdp1,0 25,3 3);title(RelTol=0.01);options=odeset(RelTol,0.00001);subplot(2,1,2);ode45(vdp

    25、1,0 25,3 3,options);title(RelTol=0.0001);n第一个图所使用的相对误差值是0.01(预设值),第二个图所使用的相对误差值是0.00001,因此我们得到较细密的点,但所花的计算时间也会比较长 0510152025-4-2024RelTol=0.010510152025-4-2024RelTol=0.0001积分输出说明n积分输出的相关处理 选用一个 OutputFcn OutputFcn 的预设值是“odeplot”,其功能为画出所有的状态变量的图像 其它可用的函式 nodephas2:画出 2-D 的相位平面(Phase Plane)nodephas3:画出 3-D 的相位平面nodeprint:随时在指令视窗显示计算结果options=odeset(OutputFcn,odephas2);%使用 odephas2 进行绘图ode45(vdp1,0 25,3 3,options);-3-2-101234-3-2-10123

    展开阅读全文
    提示  163文库所有资源均是用户自行上传分享,仅供网友学习交流,未经上传用户书面授权,请勿作他用。
    关于本文
    本文标题:Matlab在常微分方程求解中应用课件.ppt
    链接地址:https://www.163wenku.com/p-5961708.html

    Copyright@ 2017-2037 Www.163WenKu.Com  网站版权所有  |  资源地图   
    IPC备案号:蜀ICP备2021032737号  | 川公网安备 51099002000191号


    侵权投诉QQ:3464097650  资料上传QQ:3464097650
       


    【声明】本站为“文档C2C交易模式”,即用户上传的文档直接卖给(下载)用户,本站只是网络空间服务平台,本站所有原创文档下载所得归上传人所有,如您发现上传作品侵犯了您的版权,请立刻联系我们并提供证据,我们将在3个工作日内予以改正。

    163文库