数值分析课程设计--初值问题、边值问题.doc

课 程 设 计 报 告 课程名称 课题名称 专 业 班 级 学 号 姓 名 指导教师 2009年 12 月 15 日课 程 设 计 任 务 书课程名称 数值分析 课 题 初值问题、边值问题 专业班级 学生姓名 学 号 指导老师 审 批 任务书下达日期 2011 年 12 月 5 日任务完成日期 2011 年 12 月 15日 目 录一、设计内容与要求.4二、问题提出的背景.6三、问题的提出.6四、算法的理论依据及推导.7五、程序模块 1、函数模块.13 2、模块流程图.14六、实验结果与调试.24七、心得体会.27八、附件.28一、设计内容与设计要求1设计内容： 对课程计算方法中的常见算法进行综合设计或应用（具体课题题目见后面的供选题目）。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依内容的创新程度，完善程序情况及对程序讲解情况打分。2、进度安排1、 班级：2、 主讲教师：辅导教师： 上机时间安排：第 12 周 星期一 8时：30分11时：30分 星期三 8时：30分11时：30分 星期五 8时：30分11时：30分第 13 周 星期三 8时：30分11时：30分 星期五 8时：30分11时：30分 二、问题提出的背景 人类社会进入20世纪以来，在科学技术和生产力飞速发展的同时世界各国的人口也以空前的规模增长。我们知道人口是反应国情、国力基本情况的重要指标，例如评价一个国家或一个地区的发展潜力时离不开现在与今后各类人口数量比例指数和年龄分布。故人口的预测是制定和顺利实现社会经济各项战略设想的基础和出发点，制定正确人口政策的科学依据。三、问题的提出常用的用于人口预测的模型是阻滞增长模型，即Logistic模型。该模型的微分方程为： （其中：为固有增长率（即人口平均增长率），表示平均人口死亡率，表示在自然资源和环境条件允许的条件下所能容纳的最大人口数量，称为人口容量。） 利用分离变量法可以求解得： 1. 对该微分方程，分别用显示欧拉法，隐式欧拉法，改进欧拉法，梯形公式，经典龙格库塔法，理查德森外推法近似求解，得出估算人口；2. 有限差分法求解方程: 四、算法的理论基础及其推导 1、显示欧拉公式： 推导：把区间a,b分成n等分，记分点为 那么由上式有： 用向前差商法近似代替导数得： 即： 若用近似值代替，可得显示欧拉公式。 2、隐式欧拉公式： 推导：把区间a,b分成n等分，记分点为 用向后差商法近似代替导数得： 即： 亦即： 若用近似值代替，可得显示欧拉公式。3、梯形公式： 推导：把区间a,b分成n等分，记分点为 分别由显示欧拉公式和隐式欧拉公式，求其平均值可得。4、 改进欧拉方法：推导：把区间a,b分成n等分，记分点为 由显示欧拉公式和梯形公式，可得改进欧拉公式。5、 经典龙格库塔公式： 推导：由微分中值定理知 从而得到：若简单地取左端点处的斜率： 作为平均斜率的近似值，则得 若用点处的斜率近似值与右端点处的斜率近似值 的算术平均值作为斜率的近似值，则得 由上式可知，利用两个点上斜率近似值与的加权平均值作为平均斜率的近似值，利用欧拉公式预估，可构造形如： 的计算公式（称为二阶龙格库塔法），为确定系数与，将在处泰勒展开，并注意到，就有 另一方面：由上两式知： 利用类似的推导法即可求得四阶龙格库塔法.6、 理查德森外推法公式（欧拉步长自适应法）： 推导：从结点出发，先以h为步长，经一步计算出的近似值,如果使用的是p阶的方法，则： .（1）在一般情况下，式中系数c既依赖于h，又依赖于，但当h较小时，可近似地看作常数。 然后将步长折半，即以为步长，仍从出发，经两步计算求得的另一个近似值，其中每一步的截断误差约为，故：.（2）以乘以（1），并与（2）相减，得： 即： .（3）若取： .（4）作为的近似值，则其精度显然比与都要高，这种修正的思想实际上与龙贝格数值积分法思想是一致的。 若将（4）式改写为： 则由立即可得： 由此可见，若以作为的近似值，则其误差可用前后两次计算结果之差来表示，即可用： 来判断所选取色步长是否适当。具体的做法是：（1） 若（由精度要求确定），则反复将步长折半进行计算，直至，并取最后一次步长所得值作为；（2） 若，则反复将步长加倍进行计算，直至，并取上一次步长所得值作为.7、 有限差分法：将区间a,bN等分，记分点为 在内结点 处有 由上两式知： 故： 当h较小时，略去式中余项，则得确定近似值的线性方程组： 将它稍加整理后写成矩阵形式，即：其中： 方程组的，系数矩阵为对角占优的三对角线矩阵，可用追赶法解之。五、程序模块1、 函数模块1、 void GJEulerMethord(double N,int a,int b,double y0,double xm,double b2,double l) /改进欧拉公式2、 void ExplicitEulerMethord(double N,int a,int b,double y0,double xm,double b2,double l) /显示欧拉公式3、 void ImplicitEulerMethord(double N,int a,int b,double y0,double xm,double b2,double l) /隐式欧拉公式4、 void TxingMethord(double N,int a,int b,double y0,double xm,double b2,double l) /梯形公式5、 void R_KMethord(double N,int a,int b,double y0,double xm,double b2,double l) /经典龙格库塔公式6、 void JSuMethord(double N,int a,int b,double y0,double xm,double b2,double l) /理查德森外推公式7、 void FiniteDifferenceMethod(double f,double a1,double b1,double c,int max) /有限差分法2、 模块流程图（1）改进欧拉公式：开始输入：N,a,b,y00h=(b-a)/N; x=a;y=y0; m=xm/y-1;1 = kt1=y+h*f(x,y);x=x+h;y1=xm/(1+m*exp(-0.02557*x);t2=y+h*f(x,t1);y=(t1+t2)/2;输出：x,y,y1k+N结束否是（2）显示欧拉法：开始输入：N,a,b,y00h=(b-a)/N; x=a;y=y0; m=xm/y-1;1 = iy=y+h*f(x,y);x=x+h;y1=xm/(1+m*exp(-0.02557*x);i+N输出：x,y,y1结束否是开始输入：N,a,b,y00h=(b-a)/N; x=a;y=y0; m=xm/y-1;1 = iyi+1=yi+h*f(x,yi);k=0;yy=yi+h*f(x+h,yi+1)-yi+1;yp=h*(b-2*l*yi+1)-1;y1=xm/(1+m*exp(-0.02557*(x+h);fabs(yp)0.0000001输出：X+h;yi+1y1;yy=yi+1-yy/yp;是fabs(yy-yi+1)esp否是（3）隐式欧拉法：输出：X+h;yi+1y1;否yi+1=yy;k+100是x=x+h;否i+ iyi+1=yi+h*f(x,yi);k=0;yy=yi+h*f(x+h,yi+1)-yi+1;yp=h*(b-2*l*yi+1)-1;y1=xm/(1+m*exp(-0.02557*(x+h);fabs(yp)0.0000001输出：X+h;gi+1y1;yy=yi+1-yy/yp;是fabs(yy-yi+1)esp否是gi+1=yi+h/2*(f(x+h,yi)+f(x,yi+1);输出：X+h;yi+1y1;yi+1=yy;k+100是x=x+h;否i+ kk1=f(x,y);k2=f(x+h/2,y+(h*k1)/2);k3=f(x+h/2,y+(h*k2)/2);k4=f(x+h,y+h*k3);y=y+h*(k1+2*k2+2*k3+k4)/6;x=x+h;y1=xm/(1+m*exp(-0.02557*x);输出：x,y,y1k+N结束否是（6）理查德森外推法：开始输入：N,a,b,y0h=(b-a)/N; x=a;y=y0; m=xm/y-1;2 = ih=(b-a)/N;yi+1=yi+h*f(x,yi,b2,l);pi+1=yi+h/2*f(x,yi,b2,l);x=(i+1)*h;y1=xm/(1+m*exp(-0.02557*x);fabs(yi+1-pi+1)=0.01输出：x1=c;-1=c;C=1?是否否是fabs(yi+1-pi+1)=0.01输出：yi+11,y1h=h/2.0;yi+1=yi+h*f(x,yi,b2,l);pi+1=yi+h/2*f(x,yi,b2,l);yi=yi+1;h=h*2.0;yi+1=yi+h*f(x,yi,b2,l);pi+1=yi+h/2*f(x,yi,b2,l);yi=yi+1;fabs(yi+1-pi+1)=0.01输出：yi+11,y1是否i+ ini=ci/(b1i-a1i*ni-1);i+ iyi=(fi-1-a1i*yi-1)/(b1i-a1i*ni-1);i+ ixi=yi-ni*xi+1; 输出：xi结束i-=1是否六、实验结果与调试1、 选择1，改进欧拉法实现；2、 选择2，显示欧拉法实现；3、 选择3，隐式欧拉法实现；4、 选择4，梯形公式实现；5、 选择5，经典龙格库塔法实现；6、 选择6，理查德森外推法实现；7、 选择7，有限差分法解方程组；七、心得体会两周的课程设计很快就结束了，回顾这两周的课程设计，真是布满了艰辛与坎坷。由于聂老师是一个非常细致，严谨的人，对我们的要求也很严格，这使得我备受压力，深怕自己不能按时完成老师指定的任务，但是在老师细心的指导下，我逐步找到了问题的突破口，不再举手无措。在实验过程中，我都是严格按着老师的要求，先理清问题的解题思路，然后一步一步的翻译成计算机语言程序，上机操作，不断调试。最初，在调试的过程中出现了很多的错误，例如：语法错误，算法错误，层出不穷。但是在老师的指导下，经过几天的实践操作后，出现的错误逐渐减少了。特别是在写到理查德森外推法和隐式欧拉法时，花费了我很多很多的时间，我一遍一遍的翻阅教科书，在图书馆查阅相关的资料，同时不断的征求老师的意见，与周围的同学不断讨论，不断的修正。最终，问题解决了通过这次课程设计使我懂得了理论与实际相结合是很重要的，使自己的实际动手能力和独立思考的能力能得到很大的提高。虽说困难重重，但是只要你不断的坚持，只要你有足够的耐心，只要你努力，胜利就一定是属于你的。 八、附件（源代码）#include#include#includeusing namespace std;#define esp 0.01double q(double x)return 1;double r(double x)return x;double f(double x,double y,double b2,double l) return (b2*y)-(l*y*y);void p()printf(以下为近十年来的人口计算情况:nn);void GJEulerMethord(double N,int a,int b,double y0,double xm,double b2,double l)/检查正确int k;double h,x,t1,t2,y,y1,m;h=(b-a)/N;x=a;y=y0;m=xm/y-1;for(k=0;kN;k+)t1=y+h*f(x,y,b2,l);x=x+h;y1=xm/(1+m*exp(-(b2)*x);t2=y+h*f(x,t1,b2,l);y=(t1+t2)/2;printf(年份:%.0f,x);printf();printf(估计人口(千万):%f,y);printf();printf(计算人口(千万):%fn,y1);void ExplicitEulerMethord(double N,int a,int b,double y0,double xm,double b2,double l)/检查正确int i;double h,x,y,y1,m;h=(b-a)/N;x=a;y=y0;m=xm/y-1;for(i=0;i10;i+)y=y+h*f(x,y,b2,l);x=x+h;y1=xm/(1+m*exp(-(b2)*x);printf(年份:%.0f,x);printf();printf(估计人口(千万):%f,y);printf();printf(计算人口(千万):%fn,y1);void ImplicitEulerMethord(double N,int a,int b,double y0,double xm,double b2,double l)/检查正确int k,i;double h,x,y20,yy,yp,m,y1;h=(b-a)/N;y0=y0;x=a;m=xm/y0-1;for(i=0;i10;i+)yi+1=yi+h*f(x,yi,b2,l);k=0;while(k+100005)yy=yi+h*f(x+h,yi+1,b2,l)-yi+1;yp=h*(b-2*l*yi+1)-1;y1=xm/(1+m*exp(-(b2)*(x+h);if(fabs(yp)0.0000001)printf(年份:%.0f,x+h); printf(); printf(估计人口(千万):%f,yi+1);printf();printf(计算人口(千万):%fn,y1);break;yy=yi+1-yy/yp;if(fabs(yy-yi+1)esp)printf(年份:%.0f,x+h);printf();printf(估计人口(千万):%f,yi+1);printf();printf(计算人口(千万):%fn,y1);break;elseyi+1=yy;x=x+h;void TxingMethord(double N,int a,int b,double y0,double xm,double b2,double l)/检查正确int k,i;double h,x,y20,yy,yp,m,y1,g20;h=(b-a)/N;y0=y0;x=a;m=xm/y0-1;for(i=0;i10;i+)yi+1=yi+h*f(x,yi,b2,l);k=0;while(k+100)yy=yi+h*f(x+h,yi+1,b2,l)-yi+1;yp=h*(b-2*l*yi+1)-1;y1=xm/(1+m*exp(-b2)*(x+h);if(fabs(yp)0.0000001)gi+1=yi+h/2*(f(x+h,yi,b2,l)+f(x,yi+1,b2,l);printf(年份:%.0f,x+h); printf(); printf(估计人口(千万):%f,gi+1);printf();printf(计算人口(千万):%fn,y1);break;yy=yi+1-yy/yp;if(fabs(yy-yi+1)esp)gi+1=yi+h/2*(f(x+h,yi,b2,l)+f(x,yi+1,b2,l);printf(年份:%.0f,x+h);printf();printf(估计人口(千万):%f,gi+1);printf();printf(计算人口(千万):%fn,y1);break;elseyi+1=yy;x=x+h;void R_KMethord(double N,int a,int b,double y0,double xm,double b2,double l)/检查正确int k;double h,x,k1,k2,k3,k4,y,y1,m;h=(b-a)/N;x=a;y=y0;m=xm/y-1;for(k=0;kN;k+)k1=f(x,y,b2,l);k2=f(x+h/2,y+(h*k1)/2,b2,l);k3=f(x+h/2,y+(h*k2)/2,b2,l);k4=f(x+h,y+h*k3,b2,l);y=y+h*(k1+2*k2+2*k3+k4)/6;x=x+h;y1=xm/(1+m*exp(-(b2)*x);printf(年份:%.0f,x);printf();printf(估计人口(千万):%f,y);printf();printf(计算人口(千万):%fn,y1);void JSuMethord(double N,int a,int b,double y0,double xm,double b2,double l)/算法待定int i,c;double h,x,y20,p20,y1,m;h=(b-a)/N;y0=y0;x=a;m=xm/y0-1;for(i=0;i=0.01)c=1;elsec=-1;switch(c)case 1: while(fabs(yi+1-pi+1)=0.01)h=h/2.0;yi+1=yi+h*f(x,yi,b2,l);pi+1=yi+h/2*f(x,yi,b2,l);yi=yi+1;if(fabs(yi+1-pi+1)0.01)printf(估计人口(千万):%f,yi+1);printf();printf(计算人口(千万):%fn,y1);break;break;case -1: while(fabs(yi+1-pi+1)=0.01) printf(估计人口(千万):%f,yi+1);printf();printf(计算人口(千万):%fn,y1); break; ;break;void FiniteDifferenceMethod(double f,double a1,double b1,double c,int max)/有限差分法int i,k;double n10,y10,x10,y0,y1,h,x1,k1,k2;cout请输入y0y0;cout请输入y1y1;cout请输入k1k1;cout请输入k2k2;cout请输入kk;h=(k1-k2)/k;for(i=1;i=max-1;i+)/cout请输入cici;/coutendl;ci=1;for(i=1;i=max;i+)/cout请输入b1ib1i;/coutendl;x1=h*i;b1i=-2-q(x1)*h*h;for(i=2;i=max;i+)/cout请输入a1ia1i;/coutendl;a1i=1;for(i=0;imax;i+)/cout请输入fifi;/coutendl;x1=h*i;fi=r(x1)*h*h;f0=f0-y0;fmax-1=fmax-1-y1;n1=c1/b11;for(i=2;i=max-1;i+)ni=ci/(b1i-a1i*ni-1); y1=f0/b11;for(i=2;i=1;i-)xi=yi-ni*xi+1;printf();printf(x%d=%fn,i-1,xi);void main()int a=0,b=100,N=10,max,d;double l,b2;doublexm;double y0,f9;double c100,a1100,b1100;char e;printf(对logistic人口预测模型运用数值方法进行人口预测:nnn);printf(*n);printf(*1:改进欧拉法 *n);printf(*2:显示欧拉法 *n);printf(*3:隐式欧拉法 *n);printf(*4:梯形公式 *n);printf(*5:经典龙格-库塔法 *n);printf(*6:理查森外推法 *n);printf(*7:线性有限差分法 *n);printf(*nnn);doprintf(请选择进行人口预测所使用的算法（1、2、3、4、5、6、7）：n);cind;if(d=1) p();printf(请输入该国的平均人口出生率：n); cinb2;printf(请输入该国十年前的人口数量(千万)：n);ciny0; pr