数值分析课程设计.doc

课 程 设 计 报 告课程名称 数值计算方法 课题名称 数值积分 专 业 信息与计算科学 班 级 学 号 姓 名 指导教师 年 月 日湖 南 工 程 学 院课 程 设 计 任 务 书课程名称 数值分析 课 题 数值积分 专业班级 学生姓名 学 号 指导老师 审 批 任务书下达日期 年 月 日任务完成日期 年 月 日一、设计内容与设计要求1. 设计内容：1、 分别用显式，隐式，改进的欧拉公式和四阶rungekutta，求初值问题 采用步长h=0.1。2、对计算结果进行误差分析并用matlab画出对应图形。2.设计要求：l 课程设计报告正文内容a. 问题的描述及算法设计；b. 算法的流程图（要求画出模块图）；c. 算法的理论依据及其推导；d. 相关的数值结果（通过程序调试），；e. 数值计算结果的分析；f. 附件（所有程序的原代码，要求对程序写出必要的注释）。l 书写格式a要求用a4纸打印成册b正文格式：一级标题用3号黑体,二级标题用四号宋体加粗,正文用小四号宋体;行距为22。c正文的内容:正文总字数要求在3000字左右（不含程序原代码）。d封面格式如下页。l 考核方式指导老师负责验收程序的运行结果，并结合学生的工作态度、实际动手能力、创新精神和设计报告等进行综合考评，并按优秀、良好、中等、及格和不及格五个等级给出每位同学的课程设计成绩。具体考核标准包含以下几个部分：a平时出勤 （占10%）b系统需求分析、功能设计、数据结构设计及程序总体结构合理与否（占10%）c程序能否完整、准确地运行，个人能否独立、熟练地调试程序（占40%）d设计报告（占30%）注意：不得抄袭他人的报告（或给他人抄袭），一旦发现，成绩为零分。e独立完成情况（占10%）。l 课程验收要求a判定算法设计的合理性，运行相关程序，获得正确的数值结果。b回答有关问题。c提交课程设计报告。d提交软盘（源程序、设计报告文档）。e依内容的创新程度，完善程序情况及对程序讲解情况打分。二、进度安排1、 班级： 信息与计算科学：1001、1002、10032、 主讲教师：聂存云3、 时间安排：第 15 周 星期一 8时：30分11时：30分 星期三 8时：30分11时：30分 星期五 8时：30分11时：30分第 16 周 星期四 8时：30分11时：30分 星期五 8时：30分11时：30分摘要 常微分方程的数值解是建立在离散化基础上的，由此导出了化导数为差商法、taylor展开式法，euler法，runge-kutta（r-k）等常用方法，是解决工程实例的常用的工具，建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题终归结出来的微分方程主要靠数值解法数值解法就是一个十分重要的手段，而欧拉方法和r-k方法又是数值解法最基础最常用的方法。 显式，隐式，及改进欧拉(euler)方法和四阶rk公式是解决常微分方程初值问题常用的数值解法，同时还提供了四阶adams方法，simpson数值积分，milne方法，hamming方法四种线性多步法解决常微分方程初值问题。本文通过编写java程序分别实现了这八种数值解法，并通过matlab作图形式对比其精度，加深对这八种方法的认识。关键词：显示欧拉方法；隐式欧拉方法；改进的欧拉方法；四阶rk公式；四阶adams方法；simpson数值积分；milne方法；hamming方法；eclipse实现一、问题的提出1、分别用显式、隐式、改进的欧拉公式、四阶rungekutta、四阶adams方法、simpson数值积分、milne方法和hamming方法求初值问题 采用步长h=0.1，并对计算结果进行误差分析。2、对计算结果进行误差分析并用matlab画出对应图形。二、问题的分析2. 1 显示欧拉和隐式欧拉 （1）显示euler公式 将求解区间a,b分成n等份，令 ， 在 点列出方程式 （1）将用一阶差商近似代替，即，令为的近似值，将上式代入（1）中可得 （3.1.1）此式称为显示欧拉（euler）公式 （2） 隐式euler公式 在点列出方程也用一阶差商逼近导数令的近似值，则可得 ， n =0,1,2,.,n-1 （3.1.2）此式称为隐式euler公式（后退euler公式）算法流程图：i =0 = i * h = i * h0i11n y y结束注：在0处分开后，按左边的运行为显示euler公式，右边的为隐式euler公式（3） 显示及隐式euler方法的误差估计，收敛性和稳定性分析 误差估计：用不同的方法求得的数值解相对于精确解的误差是不同的，误差的大小与方法的精确度有关。人们常在为准确的前提下，即在的前提下估计误差，也就是局部截断误差。同样，也有相关定义：如果一种方法的局部截断误差为则称该方法是p阶方法或具有p阶精度。通常可用泰勒展开式来估计欧拉方法的局部截断误差，对于显示euler（3.1.1）在的条件下，有而按泰勒展开式，有 （3.1.3）将上面两式相减，得到 （）这说明显示euler的局部截断误差为,显式euler也称为一阶方法。同理将隐式euler（3.1.2）和将其按泰勒展式展开后相减，有所以，隐式euler方法也具有一阶精度收敛性：局部截断误差是在假设的条件下，由迭代公式迭代一步所产生的误差，它属于方法本身的误差。事实上，每一步迭代的局部截断误差对后面的迭代结果都会产生影响。如果去掉的假设条件，直接考察精确解与迭代解的截断误差:则称之为迭代公式的整体截断误差。如果当（同时时）时，某种数值方法的整体截断误差，则称该数值方法是收敛的，否则称为不收敛。将式（3.1.3）减去（3.1.1），可得所以 因为已经假定充分光滑，所以它满足lipschitz条件，即存在正常数l使得 这样就得到 上式给出了第i+1步的整体截断误差与第i步的整体截断误差之间的关系，它对一切i都是成立的。倘若再注意到，不难得到利用式（）得 这说明在区间上用欧拉方法求初值问题（3.1.1）的数值解，在（i=1,2,.,n）各点上的整体截断误差为也就是说，当，因为h充分小时近似值能和充分接近，数值解是收敛的。稳定性：前面关于收敛性的讨论有个前提，必须假定递推公式的计算过程是准确的。实际情况并不是这样的，递推公式的求解还会有计算误差，因为常微分方程的初值问题数值解的每步计算都是在前一步计算的结果上进行的，所以必须考虑前面误差对后面计算的影响。误差的积累不会盖过真解，这就是数值方法的稳定性问题。令其中的近似值，若不大，则可以作为的近似值；很大时，就不能用了。因此，必须分析在什么情况下不会变的很大。为此，先给出下述定义。定义： 使用一个数值方法。求解微分方程：其中，是一个复常数，对给定步长，在计算时引入了误差，若这个误差在计算后面的（k=1,2，.）中所引的误差按绝对值均不增加，就说这个数值方法对于这个步长h和复数是绝对稳定的。为了保证方法的绝对稳定，步长h和都要受到一定的限制，它们的允许范围就称为该方法的绝对稳定区域。考察欧拉方法的稳定性。把欧拉法用到方程上，得到误差方程是于是前一步误差与后一步误差之比为要求误差不增加，就要比值小于或等于1，即 （1）不等式（1）表示的是复平面上以=-1为中心，1为半径的园，如下图所示。当式（1）满足时，欧拉法是绝对稳定的。 im（） -2 re（） 0欧拉法的稳定区域绝对稳定区域越小，这个方法的适应性就越强，也不会要求h太小。若绝对稳定区域包含复平面的整个左半平面，就称这个数值方法是a稳定的。 2.2 改进的euler公式（1）改进的euler公式实质为，显式、隐式欧拉二种算法的平均值改进的欧拉公式也称为梯形公式。以为出发点，按斜率跨一个步长，直线增长到。采用这种斜率计算是为了提高计算精度，梯形公式采用了梯形积分公式，它比采用左矩形、右矩形积分公式的euler格式的精度要高。（2）梯形euler方法的误差估计，收敛性和稳定性分析泰勒展开式推导改进欧拉法的局部截断误差：当时，因为 与二阶泰勒展开式进行比较，可知这说明改进的欧拉法是二阶方法。 从验证欧拉法收敛性的过程中不难得出，改进的欧拉法的整体截断误差是，所以是收敛的，比欧拉方法高一阶;但改进的欧拉公式每积分一步要计算两次函数值，这说明精度的提高是以计算量的增加为代价的。再考虑改进欧拉法的稳定性。因为对于模型方程，它误差方程为所以所以改进的欧拉法的稳定区域如下图所示。则可知改进的欧拉法的稳定区域大于显示欧拉法的稳定区域。im（） -20re（） 改进的欧拉法的稳定区域2.3 四阶rk公式（龙咯-库塔法）设y（x）的初值问题的准确解，根据微分中值定理为解曲线y（x）在区间上的平均斜率，只要能对平均斜率提供一种近似算法，就能得到一种对应的差分格式在区间上多预测几个点的斜率，用它们的加权平均值代替平均斜率，可能构造出精度更高的公式例如取m个点的斜率构造如下形式的公式 其中是与f，n无关的常数该公式称为m级龙格库塔（rungekutta）公式，简称r-k公式注：只需要将公式的局部截断误差在点进行taylor展开，令其前面尽可能多的项为0，便可导出所满足的方程组，即可从中求出这些系数。最常用的标准四阶r-k公式为： 由于龙格-库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算法而将步长h 取小。为了讨论龙格-库塔法的绝对稳定区域，把龙格-库塔法利用到模型方程上，可以得到因此 误差方程是于是龙格库塔公式的绝对稳定区域为，在复平面画出这个区域如下图：算法流程图：i =0 =i*hi11 n y y y结束2.4 线性多步法 （1）线性多步法基本思想：在逐步推进的求解过程中，充分利用前面多步的信息来预测，则可以期望获得较高的精度。 线性多步法公式可表示为 （2.4.1）其中为的近似，这里，为常数，不全为0，则称（2.4.1）为线性k步法，计算时需要先给出前面k个近似值再由（2.4.1）逐步求出，.（2） 线性多步法构造：基于数值积分的构造方法和基于taylor展开的的构造法：将在上积分，得到只要近似的计算出右边的积分，则可通过近似。四阶adams显示公式截断误差四阶adams隐式公式截断误差simpson公式局部截断误差milne公式milne公式是四阶四步显示方法，截断误差为：hamming公式其局部截断误差为应用线性多步法求解初值问题时，开始几点处的函数值要用别的方法先计算出来。一般选用与多步法同阶或更高阶的单步法，如rungekutta方法、taylor展开方法等对线性隐式多步法，初开始几点处的函数值需要单独计算外，还需要迭代求解或预估-校正方法求解三、程序调试1. 使用说明 1）输入的所求积分的上限、下限、和等分数，如下 2）下面的八个按钮分别是八种不同的计算方法2. 运行结果预览 注：程序代码见附录 4、 结果分析4.1 计算结果分析由计算结果可以看出，与我们理论推断的结果一样，显式euler和隐式euler都只有一阶精度，为一阶方法；改进的的euler公式具有二阶精度，四阶r-k公式具有四阶精度。显式、隐式euler精度比改进后的euler精度低，改进的euler的公式比四阶r-k公式精度低。四阶adams方法、simpson数值积分、milne方法、hamming这四种线性多步法的精度都比较高，其中hamming方法的精度最高。4.2 用matlab画出这四种种计算方法和准确解的图形程序代码： x=0:0.1:1; y1=1.0 . 1.7847708324979816 ; %由于值比较多，此处省略了部分值 y2=1.0 . 1.751188711082227;y3=1.0 . 1.7378674010354138;y4=1.0 . 1.7320563651655665;y5=1 . 1.7321;plot(x,y1,r,x,y2,b,x,y3,g,x,y4,y,x,y5,k);matlab画图代码同上（此处略）注：由于这四种方法的精确点比较高，所以看起来像只有一条五、总结经过一个学期对数值分析的学习，我学习到了基本的理论知识，知道了用一些计算机语言去解决一些数学问题，这些知识都为我的课程实践和进一步学习打下了坚实的基础。通过两星期的数值分析的课程设计，我从中受益匪浅，并且对数值分析这门课程有了更深一步的认识。在实习中，我们可以把这学期所学的理论知识和实践联系起来。虽然我们对这些理论知识还不是很熟悉，但是我相信我们也能在滴水穿石中成才起来。发现问题，提出问题，解决问题，使我们从不足之处出发，寻找新的学习方向。这次课程设计我通过编程实现了常微分方程初值问题数值解法中的欧拉方法及其改进的算法和四阶龙格-库塔法并比较了其数值解与精确解之间的误差。可以看出显示、隐式、梯形欧拉公式得到的数值解精确度较差，而四阶龙格-库塔法得到的结果则相对较好。经过这两周实习，我获得了许多在课堂上听不到的知识，在这里我要谢谢老师给我们安排的课程设计实习，感谢老师对我们热心的指导和帮助。这次课程设计让我认识到了自己的不足，认识到了动手能力的重要性。我会在以后的学习中更加努力锻炼自己，提高自己，让自己写出更完善的代码。六、参考文献1 张德丰matlab数值分析与应用m北京：国防工业出版社，20092 徐长发，王邦 数值分析m武汉：华中科技大学出版社，20093 薛莲，关治，江金生数值计算方法m北京：电子工业出版社，20074 张韵华，陈效群等数值计算方法和算法m武汉：科学出版社，2000附件：/*按钮监听类，实现所有按按钮的方法*/package com.yjb1217;import java.awt.event.mouseadapter; import java.awt.event.mouseevent;import java.io.bufferedwriter; import java.io.filewriter;import java.text.decimalformat; import java.util.date;import javax.swing.joptionpane;class mouselistener extends mouseadapter private int i; static int n;/ 等分数 static double a;static double b;/ 等分下限public boolean judgetextisnull() if (login.tf1.gettext().equals() | login.tf2.gettext().equals()| login.tf3.gettext().equals() joptionpane.showmessagedialog(null, 您好，请先输入积分下限、上限以及等分数);return false;return true;public mouselistener(int i) this.i = i;public void mouseclicked(mouseevent e) decimalformat df = new decimalformat(#0.000000);if(judgetextisnull()=false)return; a = double.parsedouble(login.tf1.gettext();b = double.parsedouble(login.tf2.gettext();n = integer.parseint(login.tf3.gettext();switch (i) case 0:euler_explicit1.y = new doublen + 1;euler_explicit1.euler();login.ta.append(显式欧拉法计算的值： + + 精确值+ + 误差 + rn);for (int i = 0; i euler_explicit1.y.length; i+) login.ta.append(y+ i+ =+ t+ df.format(euler_explicit1.yi)+ t+ + df.format(exactvalue.y0i)+ t+ + df.format(math.abs(exactvalue.y0i - euler_explicit1.yi) + rn);login.ta.append(-+ rn);break;case 1:euler_implicit1.y = new doublen + 1; euler_implicit1.euler();login.ta.append(隐式欧拉法计算的值： + + 精确值+ + 误差 + rn);for (int i = 0; i euler_implicit1.y.length; i+) login.ta.append(y+ i+ =+ t+ df.format(euler_implicit1.yi)+ t+ + df.format(exactvalue.y0i)+ t+ + df.format(math.abs(exactvalue.y0i -euler_implicit1.yi)+ rn);login.ta.append(-+ rn);break;case 2:euler_trapeziod1.y = new doublen + 1;euler_trapeziod1.euler();login.ta.append(梯形欧拉法计算的值： + + 精确值+ + 误差 + rn);for (int i = 0; i euler_trapeziod1.y.length; i+) login.ta.append(y+ i+ =+ t+ df.format(euler_trapeziod1.yi)+ t+ + df.format(exactvalue.y0i)+ t+ + df.format(math.abs(exactvalue.y0i - euler_trapeziod1.yi)+ rn);login.ta.append(-+ rn);break;case 3:rk.y = new doublen + 1;rk.euler();login.ta.append(按四阶r-k计算的值为： + + 精确值+ + 误差 + rn);for (int i = 0; i rk.y.length; i+) login.ta.append(y + i + = + t + df.format(rk.yi)+ t + + df.format(exactvalue.y0i) + t+ + df.format(math.abs(exactvalue.y0i - rk.yi) + rn);login.ta.append(-+ rn);break;case 4:try filewriter fw = new filewriter(result.txt, true);bufferedwriter bw = new bufferedwriter(fw);bw.write(保存时间： + new date();bw.newline();bw.write(login.ta.gettext();bw.close();joptionpane.showmessagedialog(null,文件保存在d:myspaceshuzhifenxiresult); catch (exception e1) throw new runtimeexception(写入文件失败);break;case 5:joptionpane.showmessagedialog(null, 该程序利用常微分的差分法计算y=y-2*x/y.+ n + 在输入积分上限、下限已及等分数之后, + n + 可以利用显式、隐式、梯形欧拉公式、milne方法+ n + 四阶r-k方法、四阶adams、simpson数值微分+ n + 以及hamming方法来计算此常微分方程);break;case 6:adams.y= new doublen + 1;adams.y1= new doublen + 1;adams.adams();login.ta.append(按四阶adams计算的值为： + + 精确值+ + 误差 + rn);for (int i = 0; i adams.y.length; i+) login.ta.append(y + i + = + t + df.format(adams.yi)+ t + + df.format(exactvalue.y0i) + t+ + df.format(math.abs(exactvalue.y0i - adams.yi) + rn);login.ta.append(-+ rn);break;case 7:milne.y = new doublen + 1;mne();login.ta.append(按milne法计算的值为： + + 精确值+ + 误差 + rn);for (int i = 0; i milne.y.length; i+) login.ta.append(y + i + = + t + df.format(milne.yi)+ t + + df.format(exactvalue.y0i) + t+ + df.format(math.abs(exactvalue.y0i - milne.yi) + rn);login.ta.append(-+ rn);break;case 8:simpson .y = new doublen + 1;simpson.simpson();login.ta.append(simpson数值积分计算值： + + 精确值+ + 误差 + rn);for (int i = 0; i simpson.y.length; i+) login.ta.append(y + i + = + t + df.format(simpson.yi)+ t + + df.format(exactvalue.y0i) + t+ + df.format(math.abs(exactvalue.y0i - simpson.yi) + rn);login.ta.append(-+ rn);break;case 9:hamming.y = new doublen + 1; hamming.hamming();login.ta.append(按hamming方法计算的值为： + + 精确值+ + 误差 + rn);for (int i = 0; i hamming.y.length; i+) login.ta.append(y + i + = + t + df.format(hamming.yi)+ t + + df.format(exactvalue.y0i) + t+ + df.format(math.abs(exactvalue.y0i - hamming.yi) + rn);login.ta.append(-+ rn);break;/*程序界面设置*/package com.yjb1217;import java.awt.button;import java.awt.color;import java.awt.font;import java.awt.label;import java.awt.textarea;import java.awt.textfield;import javax.swing.jframe;public class login static textfield tf1; static textfield tf2;static textfield tf3;static label l1; static label l2; static label l3;static button bt0; static button bt1; static button bt2;static button bt3; static button bt4; static button bt5;static button bt6; static button bt7; static button bt8;static button bt9; static textarea ta; static string str;public static void main(string args) throws exception lauchframe();public static void lauchframe() jframe frame = new jframe();frame.setlayout(null);l1 = new label(积分下限); l1.setbounds(25, 5, 70, 25);l1.setfont(new font(null, 1, 12);l2 = new label(积分上限); l2.setbounds(177, 5, 60, 25);l2.setfont(new font(null, 1, 12);l3 = new label(等分数); l3.setfont(new font(null, 1, 12);l3.setbounds(329, 5, 60, 25);tf1 = new textfield(5); tf1.setbounds(107, 5, 50, 25);tf1.setbackground(new color(153, 227, 55);tf2 = new textfield(5); tf2.setbackground(new color(153, 227, 55);tf2.setbounds(259, 5, 50, 25);tf3 = new textfield(5); tf3.setbackground(new color(153, 227, 55);tf3.setbounds(405, 5, 50, 25);bt4 = new button(保存文件); bt4.setbounds(50, 410, 80, 25);bt4.setbackground(new color(50,255,50);bt5 = new button(帮助); bt5.setbounds(350, 410, 80, 25);bt5.setbackground(new color(60,255,40);bt0 = new button(显示euler公式); bt0.setbackground(new color(80,255,10);bt0.setbounds(25, 40, 95, 25);bt6 = new button(四阶adams方法); bt6.setbounds(25, 72, 95, 25);bt6.setbackground(new color(80,255,100);bt1 = new button(隐式euler公式); bt1.setbounds(139, 40, 95, 25);bt1.setbackground(new color(40,255,20);bt7 = new button(milne方法); bt7.setbounds(139, 72, 95, 25);bt7.setbackground(new color(100,255,100);bt2 = new button(梯形euler公式); bt2.setbounds(253, 40, 95, 25);bt2.setbackground(new color(60,255,30); bt8 = new button(simpson积分); bt8.setbounds(253, 72, 95, 25);bt8.setbackground(new color(120,255,85);bt3 = new button(四阶r-k公式); bt3.setbounds(367, 40, 95, 25);bt3.setbackground(new color(80,255,40);bt9 = new button(hamming方法); bt9.setbounds(367, 72, 95, 25);bt9.setbackground(new color(120,255,90);bt0.addmouselistener(new mouselistener(0);bt1.addmouselistener(new mouselistener(1);bt2.addmouselistener(new mouselistener(2);bt3.addmouselistener(new mouselistener(3);bt4.addmouselistener(new mouselistener(4);bt5.addmouselistener(new mouselistener(5);bt6.addmouselistener(new mouselistener(6);bt7.addmouselistener(new mouselistener(7);bt8.addmouselistener(new mouselistener(8);bt9.addmouselistener(new mouselistener(9);ta = new textarea();ta.setbounds(0,105, 500, 300);ta.setbackground(new color(153, 227, 55);frame.add(l1); frame.add(tf1); frame.add(l2);frame.add(tf2); frame.add(l3); frame.add(tf3);frame.add(bt0); frame.add(bt1); frame.add(bt2);frame.add(bt3); frame.add(ta); frame.add(bt4);frame.add(bt5); frame.add(bt6); frame.add(bt7);frame.add(bt8); frame.add(bt9);frame.settitle(常微分方程的差分法+ +作者：尹家兵+ +指导老师：聂存云);frame.setbounds(50, 50, 500, 480);frame.setresizable(false); frame.setvisible(true);package com.yjb1217;/* 四阶adams */public class adams static double y1; static double y;public static void adams() exactvalue.exactvalue();for(int i=0;i=3;i+)y1i=exactvalue.formula(exactvalue.xi,exactvalue.y0i);yi = exactvalue.y0i;double x1 = exactvalue.y03;for(int i=3;iy.length-1;i+)yi+1=x1+exactvalue.h/24*(55*y1i-59*y1i-1+37*y1i-2-9*y1i-3);yi+1=x1+exactvalue.h/24*(9*exactvalue.formula(exactvalue.xi+1, yi+1)+19*y1i-5*y1i-1+y1i-2);x1 = yi+1;y1i+1 = exactvalue.formula(exactvalue.xi+1, yi+1);package com.yjb1217;/* 欧拉显示法 */public class euler_explicit1 static double y ;/* euler()方法用来计算步长h，xi用来存放x的值 */public static void euler() y0 = 1; exactvalue.exactvalue();for (int i = 0; i y.length - 1; i+) yi+1=yi+exactvalue.h*exactvalue.formula(exactvalue.xi, yi) ;package com.yjb1217;/* 欧拉隐式法*/public class euler_implicit1 static double y;/* euler()方法用来计算步长h，xi用来存放x的值*/public static void euler() double y1 = new doublemouselistener.n+1;y0 = 1; y10 = 1;exactvalue.exactvalue();for (int i = 0; i y.length - 1; i+) y1i+1= y1i + exactvalue.h *exactvalue.formula(exactvalue. xi,y1i);yi + 1 = yi + exactvalue.h *exactvalue.formula(exactvalue. xi+1,y1i+1);package com.yjb