




版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领
文档简介
1、会计学1常微分方程的数值解法常微分方程的数值解法2. 把高阶方程转换成一阶微分方程组把高阶方程转换成一阶微分方程组1. 列出微分方列出微分方程程初始条件初始条件令令)0()0(210yyYyyyy21,),(),(),(2121YtfYtfYtfyyY(2.1) (2.2)(2.3),(tyyfy 00)0(,)0(yyyy第1页/共35页0)1(2yyyy 令令 yyyy21,降为一阶降为一阶1221221) 1(yyyyyyY21yyY初始条件初始条件2010210)0()0(yyyyY第2页/共35页把把t,Y作为输入宗量,把作为输入宗量,把 作为输出宗量作为输出宗量Y %M funct
2、ion file name: dYdt.m function Yd = f (t, Y) Yd = f (t,Y) 的展开式的展开式例例Van der Pol方程方程 %M function file name: dYdt.m function Yd = f (t, Y) Yd=zeros(size(Y););1 () 2 (*) 12).1 () 2 ();2 () 1 (YYYYdYYd第3页/共35页4. 使编写好的使编写好的ODE函数文件和初值函数文件和初值 供微分供微分方程解算指令(方程解算指令(solver)调用)调用0YSolver解算指令的使用格式解算指令的使用格式t, Y=s
3、olver (ODE函数文件名函数文件名, t0, tN, Y0, tol);ode45)()()()()()(222112110201tytytytytytyY输出宗量形式输出宗量形式)() 1(:,1tyY)()2(:,2tyY说明:说明:t0:初始时刻;:初始时刻;tN:终点时刻:终点时刻Y0:初值;:初值; tol:计算精度:计算精度第4页/共35页0)1(2yyyy % 主程序主程序 (程序名:程序名:VanderPol _ex1.m) t0 = 0; tN = 20; tol = 1e-6; Y0 = 0.25; 0.0; t, Y=ode45 (dYdt, t0, tN, Y0,
4、 tol); subplot (121), plot (t, Y) subplot (122), plot (Y( :, 1), Y( :, 2)解法解法1:采:采用用ODE命令命令第5页/共35页Van der Pol方程方程% 子程序子程序 (程序名:程序名: dYdt.m ) function Ydot = dYdt (t, Y)Ydot=Y(2);-Y(2)*(Y(1)2-1)-Y(1);或写为或写为function Ydot = dYdt (t, Y)Ydot=zeros(size(Y);Ydot(1)=Y(2);Ydot(2)=-Y(2)*(Y(1).2-1)-Y(1);第6页/共
5、35页第7页/共35页解法指令解题类型特 点适合场合ode45非刚性非刚性采用采用4、5阶阶RungeKutta法法大多数场合的首选算法大多数场合的首选算法ode23非刚性非刚性采用采用Adams算法算法较低精度(较低精度(103)场合)场合ode113非刚性非刚性多步法;采用多步法;采用Adams算法;高算法;高低精度均可(低精度均可(103106)ode45计算时间太长时计算时间太长时取代取代ode45ode23t适度刚适度刚性性采用梯形法则算法采用梯形法则算法适度刚性适度刚性ode15s刚性刚性多步法;采用多步法;采用2阶阶Rosenbrock算式,精度中等算式,精度中等当当ode45失
6、败时使用;失败时使用;或存在质量矩阵时或存在质量矩阵时ode23s刚性刚性一步法;采用一步法;采用2阶阶Rosenbrock算式,低精度算式,低精度低精度时,比低精度时,比ode15s有有效;或存在质量矩阵时效;或存在质量矩阵时ode23tb刚性刚性采用梯形法则反向数值微分采用梯形法则反向数值微分两阶段算法,低精度两阶段算法,低精度低精度时,比低精度时,比ode15s有有效;或存在质量矩阵时效;或存在质量矩阵时各种各种solver 解算指令的特点解算指令的特点第8页/共35页二二. 四四 阶阶 Runge-Kutta 法法00)(),(ytybtaytfdtdybtttaN110对对 I=a,
7、b作分割作分割1, 1 , 0,1Nihhhiii步长第9页/共35页初值问题的数值解法分为两大类初值问题的数值解法分为两大类单步法-Runge-Kutta 方法多步法-Admas方法计算 的近似值 时只用到 ,是自开始方法 )(1ntynnyt ,1nyuRunge-Kutta法是常微分方程的一种经典解法常微分方程的一种经典解法uMATLAB 对应命令:对应命令:ode45第10页/共35页),()2,21()2,21(),()22(6342312143211hkyhtfkkhyhtfkkhyhtfkytfkkkkkhyynnnnnnnnnn第11页/共35页四四 阶阶 Runge-Kutt
8、a 法计算流程图法计算流程图开始开始Next ifor i = 1 : N Plot)22(643211KKKKhyynn11,nnnnyytt初始条件:初始条件: ; 积分步长:积分步长: 迭代次数:迭代次数:0t0yhN输出结果httnn10ttn子程序计算iKEnd第12页/共35页第13页/共35页第14页/共35页第15页/共35页子程序:子程序:RK_sub.m function ydot = vdpol (t, y) ydot=zeros(size(y); ydot(1) = y(2); ydot(2) = -y(2)*(y(1)2-1)-y(1); 或写为:或写为: ydot
9、= y(1) ;-y(2)*(y(1)2-1)-y(1);第16页/共35页调用格式调用格式: t, y=ode23 (函数文件名函数文件名, t0, tN, y0, tol) t, y=ode45 (函数文件名函数文件名, t0, tN, y0, tol)默认精度:默认精度: ode231e-3 ode451e-6说明:说明:t0:初始时刻;:初始时刻;tN:终点时刻:终点时刻y0:初值;:初值; tol:计算精度:计算精度第17页/共35页3月月15日作业日作业: 1.Van der Pol 方程的两种解法方程的两种解法:1)采采用用ode45命令命令 2)Runge-Kutta方法方法2
10、.Duffing 方程的求解方程的求解(Runge-Kutta方法,计算步长方法,计算步长h=0.005,计算时间,计算时间t0=0.0,tN=100)要求:要求:写出程序体,打印所绘图形,图形标题用个写出程序体,打印所绘图形,图形标题用个人的名字。人的名字。)2.1cos(28.03.03tyyyy Duffing 方程方程0.0)0(,01.0)0(yy第18页/共35页第19页/共35页2. 二阶方程转成一阶方程二阶方程转成一阶方程00)0(,)0()(),()()()(XXXXXPKXXCXM nRttttt(1) 0 () 0 () 0 (,)(,)()()(2XXYYXXYnRtt
11、tt令令:)()()(tttPAYY(2)第20页/共35页12112211)()(,nnnnnnnnRttRPM0PCMKMI0A其中:其中:)()()()()(1111tttttnnnnnPM0XXCMKMI0XX 即:即:00)0(,)0()(),()()()(XXXXXPKXXCXM nRttttt)()()(tttPAYY(2)第21页/共35页)()()(tttPAYY第22页/共35页function ydot = f (t, Y,P) M=, K=, C= %输入结构参数输入结构参数 P1=zeros(3,1);inv(M)*P; A=zeros(0,0), eye(n,n)
12、; -M-1K, -M-1C ydot =AY+P1)()()(tttPAYY第23页/共35页m1m2m3k1k2k3x1x2x3k4P0sin(wt)0)sin()(02121111tPxxkxkxm 0)()(32321222xxkxxkxm 0)(3432333xkxxkxm 第24页/共35页00)sin(0000000003214333322221321321tPxxxkkkkkkkkkkxxxmmm 矩阵表示矩阵表示其中:其中:321xxxX321000000mmmM433332222100kkkkkkkkkkK00)sin()(0tPtP第25页/共35页动力学方程:动力学方程
13、:)sin(22221. 0)sin(30363. 2)sin(121088. 0)()()(000321twkPtwkPtwkPtxtxtx解析解:已知参数已知参数:m1=m2=m3=1, k1=2, k2=2, K3=1, K4=2, P0=1, 要求:要求:采用四阶龙格库塔法编程计算三个质量的响应时程.计算时间 0 50例如:例如:)()()(tttPKXXM 3第26页/共35页05101520253035404550-0.4-0.3-0.2-0.100.10.20.30.405101520253035404550-0.4-0.3-0.2-0.100.10.20.30.44阶龙格库塔法
14、的结果阶龙格库塔法的结果ode45 的结果的结果第一个质量的位移响应时程第一个质量的位移响应时程结果完全一致结果完全一致MATLAB程序程序(1)4阶阶RK方法:方法: (2)采用)采用ode45: m_chap2_ex2_1.m,m_chap2_ex2_1_sub.m 第27页/共35页例题例题3: 蹦极跳系统的动态仿真蹦极跳系统的动态仿真蹦极者系着一根弹性绳从高处的桥梁(或山崖等)向下跳。在下落的过程中,蹦极者几乎处于失重状态。按照牛顿运动规律,自由下落的物体由下式确定蹦极者系着一根弹性绳从高处的桥梁(或山崖等)向下跳。在下落的过程中,蹦极者几乎处于失重状态。按照牛顿运动规律,自由下落的物
15、体由下式确定:xxaxamgxm 210 , 00,)(xxkxxb其中,其中,m 为人体的质量,为人体的质量,g 为重力加速度,为重力加速度,x 为物体的位置,第二项和第三项表示空气的为物体的位置,第二项和第三项表示空气的阻力。其中位置阻力。其中位置 x 的基准为桥梁的基准面(的基准为桥梁的基准面(即选择桥梁作为位置的起点即选择桥梁作为位置的起点 x0),低于桥),低于桥梁的位置为正值,高于桥梁的位置为负值。梁的位置为正值,高于桥梁的位置为负值。如果人体系在一个弹性常数为如果人体系在一个弹性常数为 k 的弹性绳索的弹性绳索上,定义绳索下端的初始位置为上,定义绳索下端的初始位置为 0,则其对落
16、,则其对落体位置的影响为:体位置的影响为:地面地面x桥梁基准面桥梁基准面0梯梯子子h2h1空气的阻力空气的阻力第28页/共35页xxaxaxbmgxm 21)(整个蹦极系统的数学模型为:整个蹦极系统的数学模型为:设桥梁距离地面为设桥梁距离地面为 50 m,即,即 h2=50,蹦极者的起始位置为绳索的长度,蹦极者的起始位置为绳索的长度 30 m,即,即 h1=30,蹦极者起始速度为,蹦极者起始速度为 0,其余的参数分别为,其余的参数分别为 k20, a2a11;m70 kg,g10 m/s2。地面地面x桥梁基准面桥梁基准面0梯梯子子h2h1初始条件:初始条件:0)0(;30)0(xx已知参数:已
17、知参数:10,70, 1,20,50,302121gmaakhh0 , 00,)(xxkxxb第29页/共35页xxY0,00100,0102121xgYxmamaYxgYxmamamkY令令: 0300)0()0(xxY初始条件变为:初始条件变为:xxaxaxbmgxm 21)(0 , 00,)(xxkxxb第30页/共35页y0=-30; 0; % 初始位移和初始速度初始位移和初始速度t,y=ode45(bengji_sub, 0:0.01:100, y0);x1=50. - y(:,1); % x1代表蹦极者与地面之间的距离代表蹦极者与地面之间的距离plot(t,x1); gridplo
18、t(t,y(:,1); grid % y(:,1)代表位移代表位移主程序主程序 (程序名:bengji.m)Matlab程序程序第31页/共35页function ydot=f(t,y)m=70; k=20; a1=1; a2=1; g=10; x=y(1); % x代表蹦极者的位移代表蹦极者的位移x_dot=y(2); % x_dot 代表代表 x 的速度的速度if x0 ydot=0,1;-k/m,-a1/m-(a2/m)*abs(x_dot)*y+0;g;else ydot=0,1;0,-a1/m-(a2/m)*abs(x_dot)*y+0;g;end子程序子程序 (程序名:(程序名:bengji_sub.m )0,00100,0102121xgYxmamaYxgYxmamamkY第32页/共35页y(:,1)x1结果分析:结果分析: 右上图为蹦极者与地面之间的距离。从结果可看出,对
温馨提示
- 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
- 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
最新文档
- 消防维修采购合同范本
- 水泵大修合同范本
- N-N-didesmethyl-AH-7921-生命科学试剂-MCE
- 眼部肿瘤患者的营养支持与饮食调整
- 1-2-Distearoyl-3-butyryl-glycerol-TG-18-0-18-0-4-0-生命科学试剂-MCE
- 现代办公楼宇与生态环境的融合设计案例
- 电子器件与人工智能的交叉应用研究
- 科技白水泥品牌的国际市场拓展策略
- 房屋权属转移合同范本
- 购买枣树合同范本
- 2025年海域使用权租赁合同
- 四年级希望杯历年数学竞赛试题与答案1-13届+奥数分类专项练习集等
- 《走近世界民间美术》 课件 2024-2025学年人美版(2024)初中美术七年级下册
- (2025春)人教版三年级数学下册全册教案
- 2025年江苏省高职单招《职测》高频必练考试题库400题(含答案)
- 2025云南红河州个旧市大红屯粮食购销限公司招聘及人员高频重点模拟试卷提升(共500题附带答案详解)
- X证书失智老年人照护讲解
- 传感器技术-武汉大学
- 2024新版有限空间作业安全大培训
- (正式版)JBT 14449-2024 起重机械焊接工艺评定
- 骆驼祥子选择题100道及答案
评论
0/150
提交评论