还剩14页未读,继续阅读
本资源只提供10页预览,全部文档请下载后查看!喜欢就下载吧,查找使用更方便
文本内容:
数学与信息科学系实训报告常微分方程初值问题数值解法及其__TLAB实现姓名学号专业信息与计算科学年级2010级指导教师朱耀生完成时间2013年11月25—2013年12月6实验目的:1研究满足给定方程的可微函数的数值方法2培养__tlab编程和上机调试能力实验基本原理和内容:根据给定的初始条件,确定常微分方程惟一解的问题叫常微分方程初值问题大多数实际问题难以求得解析解,必须将微分问题离散化,用数值方法求其近似解 一阶常微分方程的初值问题的提法是,求出函数,使满足条件 1利用数值方法解问题1时,通常假定解存在且惟一,解函数及右端函数具有所需的光滑程度数值解法的基本思想是:先取自变量一系列离散点把微分问题1离散化,求出离散问题的数值解并以此作为微分问题解的近似例如取步长0以剖分区间[],令=+把微分方程离散化成一个差分方程以表微分方程初值问题的解以表差分问题的解,就是近似解的误差,称为全局误差因此,设计各种离散化模型求出近似解估计误差以及研究数值方法的稳定性和收敛性等构成了数值解法的基本内容
1.欧拉法和后退欧拉法原理:function[xy]=meulerdfxspany0h%用途改进欧拉公式解常微分方程y=fxyyx0=y0%格式[xy]=meulerdfaby0hdf为函数fxyxspan为求解%区间[x0xn]y0为初值yx0h为步长[xy]返回节点和数值解矩阵x=xspan1:h:xspan2;y1=y0;forn=1:lengthx-1k1=fevaldfxnyn;yn+1=yn+h*k1;k2=fevaldfxn+1yn+1;yn+1=yn+h*k1+k2/2;end例题给定的初值问题y′=-y+x+2,0=x=1Y0=-1取精确解yx=exp-x+x后退欧拉法,步长h=
0.003 h=
0.1求在节点k=1+
0.1k k=123……10处的数值解若h=
0.1; y=-1; x=1; for i=1:20 k1=h*oulei_wfxy; k2=h*oulei_wfx+hy+k1; y=y+
0.5*k1+
0.5*k2x=x+h; z=ouleij_qx t=y-z endz = 4z =
3.7000y = -
0.6150 z =
1.4329 z =
1.4329 t = -
2.0479 z =
3.7150 z =
3.4435 y = -
0.2571z =
1.5012 z =
1.5012 t =-
1.7583 z =
3.4571 z =
3.2___ y =
0.0763z =
1.5725 z =
1.5725 t =-
1.4962z =
3.2237 z =
3.0013 y =
0.3876z =
1.6466 z =
1.6466 t =-
1.2590z =
3.0124z =
2.8112y =
0.6788 z =
1.7231z =
1.7231t =-
1.0444 z =
2.8212z =
2.6391y =
0.9518 z =
1.8019z =
1.8019t =-
0.8501z =
2.6482z =
2.4834y =
1.2084z =
1.8827z =
1.8827 t =-
0.6743 z =
2.4916z =
2.3425 y =
1.4501 z =
1.9653z =
1.9653 t =-
0.5152实验结果的分析与研究
1.对于欧拉法,步长越小,精度越高,而产生的误差越小总对于欧拉法,步长越小,精度越高,而产生的误差越小对于欧拉法体来说,欧拉法的优点是形式简单计算方便,是形式简单,体来说,欧拉法的优点是形式简单,计算方便,缺点是总的运算精度比较低的增大,误差值也越来越大算精度比较低而且随着x的增大,误差值也越来越大根据欧拉公式的截断误差计算,欧拉法是一阶方法欧拉公式的截断误差计算,欧拉法是一阶方法
2.对于改进欧拉法,对于改进欧拉法,其基本特征与欧拉法相似,也是步长越小,对于改进欧拉法其基本特征与欧拉法相似,也是步长越小,精度越高,误差越小优点是精度相对欧拉法来说较高,精度越高,误差越小优点是精度相对欧拉法来说较高,截断误差为Oh^3,缺点是比欧拉法计算量大根据欧拉改进法,缺点是比欧拉法计算量大公式的截断误差计算,欧拉改进法是二阶方法公式的截断误差计算,欧拉改进法是二阶方法2.龙格-库塔方法原理对于龙哥库塔三阶利用Yn+1 = Yn + h/6K1 + 4K2 + K3 K1 = f(Xn Yn)K2 = f(Xn + 1/2*h Yn + h/2*K1)K3= f(Xn + h Yn - K1+2*h*k2)若 clear all h=
0.1; y=-1;x=1; for i=1:25;k1=oulei_wfxy; k2=oulei_wfx+h/2y+k1*h/2;k3=oulei_wfx+hy-k1*h+k2*2*h; y=y+k1+4*k2+k3*h/6; x=x+h; z=ouleij_qx; t=absy-z; A=[x y z t] endz = 4z =
3.8500z =
3.7300A =
1.1000 -
0.6145
1.4329
2.0474z =
3.7145z=
3.5788z =
3.4702z =
1.5012 A =
1.2000 -
0.2562
1.5012
1.7574 z =
3.4562z =
3.3334z =
3.2351z =
1.5725 A =
1.3000
0.0776
1.5725
1.4950z =
3.2224z =
3.1113z =
3.0224z =
1.6466 A =
1.4000
0.3__1
1.6466
1.2575 z =
3.0109z =
2.9104z =
2.8299z =
1.7231 A =
1.5000
0.6804
1.7231
1.0427z =
2.8196z =
2.7286z =
2.6558z =
1.8019 A =
1.6000
0.9536
1.8019
0.8483z =
2.6464z =
2.5641z =
2.4982z =
1.8827 A =
1.7000
1.2103
1.8827
0.6724 z =
2.4__7z =
2.4152z =
2.3556z =
1.9653 A =
1.8000
1.4521
1.9653
0.5132 z =
2.3479z =.2805z =
2.2266z =
2.0496 A =
1.9000
1.6803
2.0496
0.3692 例题1给定的初值问题y′=-y+x+2,0=x=1Y0=-1取精确解yx=exp-x+x龙格库塔法三阶方法 h=
0.1求在节点k=1+
0.1k k=123……10处的数值解例题2给定的初值问题y′=-y+x+2,0=x=1Y0=-1取精确解yx=exp-x+x龙格库塔四阶法 h=
0.1求在节点k=1+
0.1k k=123……10处的数值解原理function[xy]=m4rkodesdfxspany0h%用途:4阶经典龙格库塔公式解常微分方程组y=fxyyx0=y0%格式:[xy]=m4rkodesdfxspany0h%df为向量函数fxy表达式xspan为求解区间[x0xn]%y0为初值向量h为步长x返回节点y返回数值解向量x=xspan1:h:xspan2;y=zeroslengthy0lengthx;y:1=y0:;forn=1:lengthx-1k1=fevaldfxny:n;k2=fevaldfxn+h/2y:n+h/2*k1;k3=fevaldfxn+h/2y:n+h/2*k2;k4=fevaldfxn+1y:n+h*k3;y:n+1=y:n+h*k1+2*k2+2*k3+k4/6;endx=x;y=y;若h=
0.1; y=-1; x=1; for i=1:20 k1=h*oulei_wfxy; k2=oulei_wfx+h/2y+k1*h/2;k3=oulei_wfx+h/2y+k2*h/2;k4=oulei_wfx+hy+h*k3;y=y+k1+2*k2+2*k3+k4*h/6;x=x+h;z=ouleij_qx;t=absy-z;A=[xyzt]endz=4z=
3.8500z=
3.8575z=
3.7142z=
1.4329A=
1.1000-
0.
61451.
43292.0407z=
3.7145z=
3.5788z=
3.5856z=
3.4560z=
1.5012A=
1.2000-
0.
25621.
50121.7574z=
3.4562z=
3.3334z=
3.3395z=
3.2222z=
1.5725A=
1.
30000.
07751.
57251.4950z=
3.2225z=
3.1113z=
3.1169z=
3.0108z=
1.6466z=
1.
40000.3__
01.
64661.2576z=
3.0110z=
2.9104z=
2.9154z=
2.8194z=
1.7231A=
1.
50000.
68041.
72311.0427z=
2.8196z=
2.7286z=
2.7332z=
2.6463z=
1.8019A=
1.
60000.
95361.
80190.8483z=
2.6464z=
2.5641z=
2.5682z=
2.4__6z=
1.8827A=
1.
70001.
21021.
88270.6724z=
2.4__8z=
2.4153z=
2.4190z=
2.3479z=
1.9653A=
1.
80001.
45201.
96530.5133z=
2.3480z=
2.2806z=
2.2840z=
2.2196z=
2.0496A=
1.
90001.
68032.
04960.3693z=
2.2197z=
2.1587z=
2.1618z=
2.1035z=
2.1353A=
2.
00001.__
642.
13530.2390z=
2.1036z=
2.0485z=
2.0512z=
1.9985z=
2.2225A=
2.
10002.
10142.
22250.1211z=
1.9986z=
1.9487z=
1.9512z=
1.9035z=
2.3108A=
2.
20002.
29642.
31080.0144结果的分析与总结对于龙格—库塔法,优点是精度更高,同样的步长下精度比对于龙格—对于龙格库塔法,优点是精度更高,欧拉法高的多,欧拉法高的多,误差更小,截断误差为Oh^5,故龙格—库,故龙格—塔法是四阶方法缺点是每步都要计算四次微分值塔法是四阶方法缺点是每步都要计算四次微分值
3.亚当斯法%程序function[xy]=m4adamsdfxspany0h%用途4阶亚当斯预报-校正格式解常微分方程y=fxyyx0=y0%格式[xy]=m4adamsdfxspany0h%df为函数fxy表达式xspan为求解区间[x0xn]%y0为初值h为步长x返回节点y返回数值解x=xspan1:h:xspan2;[x1y]=m4rungekuttadf[x1x4]y0h;forn=4:lengthx-1p=yn+h/24*55*fevaldfxnyn-59*fevaldfxn-1yn-
1...+37*fevaldfxn-2yn-2-9*fevaldfxn-3yn-3;yn+1=yn+h/24*fevaldfxn-2yn-2-5*fevaldfxn-1yn-
1...+19*fevaldfxnyn+9*fevaldfxn+1p;end综合实例市场__模型对于纯粹的市场经济来说商品市场__取决于市场供需之间的关系市场__能促使商品的供给与需求相等这样的__称为静态均衡__.也就是说如果不考虑商品__形成的动态过程那么商品的市场__应能保证市场的供需平衡但是实际的市场__不会恰好等于均衡__而且__也不会是静态的应是随时间不断变化的动态过程.试建立描述市场__形成的动态过程的数学模型假设在某一时刻商品的__为它与该商品的均衡__间有差别此时存在供需差此供需差促使__变动.对新的__又有新的供需差如此不断调节就构成市场__形成的动态过程假设__的变化率与需求和供给之差成正比并记为需求函数为供给函数(为参数)于是其中为商品在时刻的__为正常数.若设则上式变为
①其中均为正常数其解为.下面对所得结果进行讨论:
(1)设为静态均衡__则其应满足即于是得从而__函数可写为令取极限得这说明市场__逐步趋于均衡__.又若初始__则动态__就维持在均衡__上整个动态过程就化为静态过程;
(2)由于所以当时单调下降向靠拢;当时单调增加向靠拢.这说明初始__高于均衡__时动态__就要逐步降低且逐步靠近均衡__;否则动态__就要逐步升高.因此式
①在一定程度上反映了__影响需求与供给而需求与供给反过来又影响__的动态过程并指出了动态__逐步向均衡__靠拢的变化趋势.总结通过这次与众不同的计算方法课程的学习,我和合作同学在此间都收获了很多从我们的角度将从以下几点谈论对这个学期计算方法学习做一个总结
1、首先从学习收获的角度从我们的交流过程中可以得到几点首先,我们对运用计算解决数学问题有了一定的了解,以前大家学习高数的时候整天研究做题,但是对高数中的积分微分等方法的原理没有一个基础的理解但是这次,我们能深刻的理解到计算的基本原理,这样对知识的掌握有一个新的高度其次,除了计算原理理解之外,在计算方法的学习和编程实现中也渐渐的明白了一个理论的公式计算(如同常微分中复杂的公式,高难度的运算)如何转化成程序去实现在这个过程中,明白了计算机运行,或者说是设计一个算法的基础是什么,简单的说是借用计算机告诉可重复的运算来弥补其中其离散的缺点这样虽然不能带到醉精确的数值但是可以得到一个你需要的数值在这个过程中我们从教材上,从课堂上去体会是怎么把一个问题去离散化,如何去用迭代等方式去逼近一个问题等等通过学习我们基本上认识到了怎么将一个现实的问题去用编程思想去考虑这不仅仅体现在我们算法的理解上面,更是我们对编程的一次新的理解
2、其次从授课方式角度曾老师改变以前的授课和考核方式,以小组的方式经行考核,这样对我们来说也是一个挑战但是从最终大家的结果来看收获颇丰第
一、让大家对教材的知识有了一个全新的认识采用阶段性的小组汇报,这样让大家有时间经行讨论,在大家写程序的时候对教材的知识有了一遍了解,在给大家讲解的时候有了一个更深的了解,然后到最后的报告总结,当大家能用自己的话把课本上的很多东西写出来的时候,或许大家不仅仅是看了教材上的内容,还会参考很多课外的资料,我想这个时候他们的理解会是非常深刻的第
二、培养了大家的团队合作能力,在这次的计算方法的实习中,我们从写程序到最后的报告总结始终贯彻着团队的精神,在写程序初期,我们两人一组,相互的补充相互的学习,使得这周的讨论效率更高,在讨论好大家之后,大家的分工明确在完成自己的工作得到同时主动帮助其他人完成任务于此同时大家相互的交换工作相互的知道让每个人有了很快的成长大家在团队中清楚了自己的角色,如何把一个比较大的工程分解到每个人身上,再怎么整合起来成为一个需要的产品这不仅仅在这次的学习中很重要,也渐渐的明白了团队的工作第
三、大家热情高涨,收获效率更上一层楼,像其他的组一样,开始的时候我们组也是存在部分,基础比较差编程能力相对来说比较弱,我们在按照老师的授课方式的基础上,在小组内实行两人一组等一系列的活动安排,调动我们的热情,渐渐的我们开始积极的展示自己的代码,对讨论的后的安排工作积极完成,一度出现了预约下一次任务的场景虽然大家做任务的水平参差不齐,但是还是本着大家自愿保留大家成果的原则经行小组汇报这样可能每次实习报告不是最好的,但是极大的调动了我们的积极性,收获肯定是最高的。