实验四

内江师范学院数值分析实验报告册编制张莉审定牟廉明专业：班级：级班学号：姓名：数学与信息科学学院2013年9月说明一、学生在做实验之前必须要准备实验，主要包括预习与本次实验相关的理论知识，熟练与本次实验相关的软件操作，收集整理相关的实验参考资料，要求学生在做实验时能带上充足的参考资料；若准备不充分，则学生不得参加本次实验，不得书写实验报告；二、要求学生要认真做实验，主要是指不得迟到、早退和旷课，在做实验过程中要严格遵守实验室规章制度，认真完成实验内容，极积主动地向实验教师提问等；三、各个实验按照学生水平分别设置了A、B、C、D四个等级，其中对应的难度系数为1、0.8、0.7、0.6，也可根据实际完成情况制定相应地的难度系数，但总体保证难度排序为A级难度最大，B级次之，C级较易，D级最简单。四、学生可以根据自己对各个实验涉及到的知识点掌握的程度自由选取A、B、C、D等级的实验题目。五、学生要认真工整地书写实验报告，实验报告的内容要紧扣实验的要求和目的，不得抄袭他人的实验报告；四、根据实验准备、实验态度、实验报告的书写、实验报告的内容进行综合评定，并给出实验成绩评定分。实验名称：实验四常微分方程的数值算法指导教师：吴开腾张莉实验时数：2实验设备：安装了Matlab、C++、VF软件的计算机实验日期：2013年11月13日实验地点：第五教学楼北902实验目的：1.掌握常微分方程的各种差分方法的基本思想和基本步骤。2.理解各类差分方法的优缺点，并能自行编程求解，体会参数和步长对问题解的影响。实验准备：1.在开始本实验之前，请回顾教科书的相关内容；2.需要一台准备安装WindowsXPProfessional操作系统和装有数学软件的计算机。实验内容及要求A题求解常微分方程初值问题d1,01d(0)1yyxxxy其中5050。（1）对参数取不同的值，取步长0.01h，编写程序用4阶经典R-K法计算，将计算结果画图比较，并分析相应的初值问题的性态。（2）取参数为一个绝对值不大的负值和两个不同的步长，一个步长使参数,h在经典R-K法的稳定域内，另一步长在经典R-K法的稳定域外，分别用R-K法编写程序计算并比较计算结果，取全域等距的10个点上的计算值。（3）用MATLAB中的内部函数dsolve求此常微分方程初值问题的解并与（1）与（2）中的结果进行比较。B题求解微分方程组d2000999.751000.25(0)0,(0)2uuvtvuvtuvddd任选一显式方法，编写程序取不同步长求解上述方程。C题给定初值问题2d2,12d(0)1xyyxexxxy其精确解为2()xyxee。试用改进的欧拉法，步长0.05,0.1hh，编写程序求在节点10.1(1,2,10)kxkk处的数值解及误差，并用MATLAB中的内部函数dsolve求此常微分方程初值问题的解并与上述结果进行比较。D题给定初值问题2d2,12d(0)1xyyxexxxy其精确解为2()xyxee。试用欧拉法，步长0.025,0.1hh，编写程序求在节点10.1(1,2,10)kxkk处的数值解及误差，并用MATLAB中的内部函数dsolve求此常微分方程初值问题的解并与上述结果进行比较。说明：实验过程应包括对问题的简要分析、求解方法、求解步骤、程序及其必要的图表等内容。实验过程：A题的实验过程1、实验中问题的重述问题1：对参数取不同的值，取步长0.01h，编写程序用4阶经典R-K法计算，将计算结果画图比较，并分析相应的初值问题的性态。问题2：取参数为一个绝对值不大的负值和两个不同的步长，一个步长使参数,h在经典R-K法的稳定域内，另一步长在经典R-K法的稳定域外，分别用R-K法编写程序计算并比较计算结果，取全域等距的10个点上的计算值。问题3：用MATLAB中的内部函数dsolve求此常微分方程初值问题的解并与（1）与（2）中的结果进行比较。2、问题一的求解过程2.1问题一的分析首先，需要将参数赋不同值，于是利用for循环，取1为步长。并且将分为大于0和小于0两类，分别进行讨论。然后，采用4阶经典R-K法计算此常微分。最后，利用plot函数画出图形。2.2问题一的求解程序（见附录）2.3问题一的求解结果0000.10.20.30.40.50.60.70.80.9100.511.522.500.10.20.30.40.50.60.70.80.910123456x10213、问题二的求解过程3.1问题二的分析首先，取参数为一个绝对值不大的负值和两个不同的步长，此处取为-40，步长分别为0.01和0.2；然后，利用R-K法编写程序计算；最后，取等距的10个点上的计算值。3.2问题二的求解程序（见附录）3.3问题二的求解结果h=0.01且a=-40h=0.2且a=-4000.10.20.30.40.50.60.70.80.910.10.20.30.40.50.60.70.80.9100.10.20.30.40.50.60.70.80.91024681012141618x1094、问题三的求解过程4.1问题三的分析必须利用matlab内部函数dsolve求解此常微分方程，此处取a=-40。4.2问题三的求解程序（见附录）4.3问题三的求解结果00.10.20.30.40.50.60.70.80.910123456x1021B题的实验过程1、实验中问题的重述求解微分方程组d2000999.751000.25(0)0,(0)2ddduuvtvuvtuv任选一显式方法，编写程序取不同步长求解上述方程。2、实验分析分析该题题意可得到，该题是一个已知初值与未知量所满足的两个微分方程所组成的一阶微分方程组。需要我们做的就是去应用或构造一种Euler显式方法对其微分方程组进行求解。分析该一阶微分方程组的求解方法，引进,0nhttnn=1,2,···，以un和vn表示节点tn上的近似值，此处应用Euler方法中的梯形格式形式进行求解，即将显式Euler格式与隐式Euler格式的结合，是其两者的算术平均。这样既在一定程度上提高了计算精度，也缩减了计算过程的复杂性，计算简单可行。3、实验求解过程首先，确定初值。由已知已经得到u(0）=0;v(0)=-2。除此之外，对向量t，u，v赋初值。假定所需得到的结果个数为n，则将t，u，v均赋初值为一个含有n+1个元素均为0的向量。然后，确定步长。假定所求方程组的值的范围是在[a,b]内，所需得到的结果个数为n，则可得到步长h为(b-a)/n。其次，引用梯形格式。令f(u,v)=-2000u+999.75v+1000.25;g(u,v)=u-v。引进初值问题：00dy(,)()fxydxyxy，其中n=1,2,···，以un和vn表示节点tn上的近似值，则其梯形公式具有以下形式：)],(),([2111nnnnkkvufvufhuu；)],(),([2111nnnnkkvugvughvv4、实验求解结果通过对梯形公式的运用，在Matlab的命令窗口输入a,b,n的具体取值，则可运算出u,v的具体数值值，即可得到向量u和v，见下表：表1参数值C题的实验过程1、实验中问题的重述给定初值问题2d2,12d(0)1xyyxexxxy其精确解为2()xyxee。试用改进的欧拉法，步长0.05,0.1hh，编写程序求在节点10.1(1,2,10)kxkk处的数值解及误差，并用MATLAB中的内部函数dsolve求此常微分方程初值问题的解并与上述结果进行比较。2、实验分析该问题给定初值求解常微分问题,并用改进的欧拉法计算,且需利用不同的步长进行处理比较,其主要目的是研究步长对求解常微分问题的影响程度.tuv00-20.10.05-1.89750.20.22371875-1.7889390630.3-11.54774764-2.2656935540.41129.9956653.769986720.5-109081.7181-5346.5046340.610531931.71516063.32040.7-1016864424-49826364.480.89817888273148107658690.9-9.47923E+12-4.64482E+1119.15226E+144.48461E+133、实验求解过程(代码见附录)利用MATLAB编程实现求解过程，程序见附录。4、实验求解结果步长和数值解的比较如下：x1.051.11.151.21.251.31.351.41.45y0.150.340.580.861.201.602.072.613.241.51.551.61.651.71.751.81.851.91.9523.964.785.716.767.959.2810.7712.4414.3016.3618.65x1.11.21.31.41.51.61.71.81.92y0.340.861.592.603.945.687.9110.7214.2418.58步长不同对误差的影响：单位（E-03）x1.051.11.151.21.251.31.351.41.45y3.0032.9052.8082.7142.6222.5322.4452.3602.2791.51.551.61.651.71.751.81.851.91.952.2002.1242.0501.9791.9101.8441.7791.7181.6581.600x1.11.21.31.41.51.61.71.81.9y10.249.6099.0008.4197.8697.3516.8646.4075.979由以上数据可知步长越小其所对应的误差越小。D题的实验过程1、实验中问题的重述给定初值问题2d2,12d(0)1xyyxexxxy其精确解为2()xyxee。试用欧拉法，步长0.025,0.1hh，编写程序求在节点10.1(1,2,10)kxkk处的数值解及误差，并用MATLAB中的内部函数dsolve求此常微分方程初值问题的解并与上述结果进行比较。2、实验分析这道题主要考察用欧拉法编写程序求此常微分方程初值问题，对于误差分析，主要采用局部截断误差，因为实验中的累积误差较复杂。方程中含有导数项，这是微分方程的本质特征，导数是极限过程的结果，而计算过程则总是有限的，所以数值解法的第一步就是消除导数项dy/dx,即y’,这个过程也叫离散化。算法原理：对于初值问题00dy(,)()fxydxyxy,(),nnnnxyyx假设第步在点的值计算没有误差即1,ny由单步法计算出则111()nnnTyxy称为点1nx上的局部截断误差。欧拉法：若用差商hxyxynn)()(1替代其中的导数项y’(nx),则有近似关系式：))(()()(,1nnnnxyxhfxyxy若用ny替代)(nxy，这样设计出的计算公式为：1ny=ny+),(nnyxhf则该初值问题变为：1ny=ny+),(nnyxhf;)(00xyy局部截断误差为：23211()()(,(),)()()()26yxhyxhxyxhhyxhyxOh所以欧拉方法是一阶方法，其局部截断误差为)(2hO。3、实验求解过程分别令h=0.025和h=0.01,带入欧拉公式，当h=0.025和h=0.01时，求出x=(1:h:2)对应的ny的值，并画出图像。打开matlab软件，新建一个空白M文件，输入如下程序并保存（见附录）；然后，在commandwindow中输入以下程序，并画出图像如下：函数真实值与欧拉步长为h=0.025和h=0.01的误差：4、实验求解结果结果分析：由以上结果可知：采取欧拉法，步长取得越