椭圆型微分方程.docx

1、椭圆型微分方程数学与计算科学学院实 验 报 告实验项目名称 椭圆型方程数值解 所属课程名称 微分方程数值解法 实 验 类 型 验证 实 验 日 期 班 级 信计0902 学 号 姓 名 成 绩 一、实验概述：【实验目的】1掌握椭圆型方程的五点差分方法2掌握求解线性方程组的Jacobi迭代法【实验原理】【实验环境】Microsoft Visual C+6.0二、实验内容：【实验方案】用五点差分格式的Jacobi迭代法求解单位正方形区域上的方程边值问题【实验过程】（实验步骤、记录、数据、分析）1，部分主要算法代码：(1)Dirichlet边值、系数矩阵A以及向量序列K的初步赋值for(i=0;i=

2、m;i+)/Dirichlet边界值 Ui0=f(xi,y0); Uim=f(xi,ym); U0i=f(x0,yi); Umi=f(xm,yi); for(i=0;in;i+) for(j=0;jn;j+) aij=0; Ki=0; for(i=0;in;i+)/系数矩阵 for(j=0;jn;j+) if(i=j) aij=4; else if(j-i=1&j%(m-1)!=0) aij=-1; aji=-1; else if(j-i=m-1) aij=-1; aji=-1; for(i=0;i0&in-m+1&in-1) Ki=Ui-(m-2)*(m-1)+1m; else if(i=n

3、-1) Ki=Um-1m+Umm-1; for(j=1;jm-2;j+) K(m-1)*j=U0j+1; K(m-1)*(j+1)-1=Umj+1; (2)Jacobi迭代法int Jacobi(int ann,double bn,double xn) int i,j; double Fan8;/定义向量无穷范数 double sum1,sum2;/定义累加器 int count=0;/计数器 double x0n,x1n; for(i=0;in;i+) x0i=1;/读入初值序列 for(i=0;in;i+) x1i=x0i;/初始化 do for(i=0;in;i+) x0i=x1i; f

4、or(i=1;i=n;i+) for(j=1,sum1=0;j=i-1;j+) sum1=sum1+ai-1j-1*x0j-1; for(j=i+1,sum2=0;j=n;j+) sum2=sum2+ai-1j-1*x0j-1; x1i-1=(bi-1-sum1-sum2)/ai-1i-1;/Jacobi迭代公式 count+; for(Fan8=fabs(x10-x00),i=0;in;i+) if(Fan8e); for(i=0;in;i+) xi=x1i;/将结果返回 printf(Jacobi方法迭代次数为 %d次!nn,count); return 0;(3)Gauss_seidel

5、迭代法仅将Jacobi迭代代码的sum1=sum1+ai-1j-1*x0j-1处改为sum1=sum1+ai-1j-1*x1j-1;(4)SOR迭代法int SOR(int ann,double bn,double xn) int i,j; int count=0;/计数器 double Fan8;/定义向量无穷范数 double sum1,sum2;/定义累加器 double x0n,x1n; /初值序列 for(i=0;in;i+) x0i=1; for(i=0;in;i+) x1i=x0i;/初始化 do for(i=0;in;i+) x0i=x1i; for(i=1;i=n;i+) f

6、or(j=1,sum1=0;j=i-1;j+) sum1=sum1+ai-1j-1*x1j-1; for(j=i,sum2=0;j=n;j+) sum2=sum2+ai-1j-1*x0j-1; x1i-1=x0i-1+w*(bi-1-sum1-sum2)/ai-1i-1;/SOR迭代公式 count+; for(Fan8=fabs(x10-x00),i=0;in;i+) if(Fan8e); for(i=0;in;i+) xi=x1i;/将结果返回 printf(SOR方法迭代次数为 %d次!nn,count); return 0;(5)边值函数（解析解）double f(double a,d

7、ouble b)/边值函数 double fun; fun=pow(a,2)-pow(b,2); return fun;2，结果分析：1，A为对称占优矩阵，且为对称正定矩阵，存在唯一解。2，迭代速度从低到高依次为：Jacobi迭代法、Gause_Seidel迭代法，SOR迭代法。【实验结论】（结果）【实验小结】（收获体会） 通过本次实验，我重新复习了Jacobi，Gause_Seidel以及SOR迭代法，也对课本上的正方形区域中的Laplace方程Dirichlet边值问题有了更深的了解，同时也提高了自己的上机编程能力，感受到了理论与实践相结合的乐趣。三、指导教师评语及成绩：评 语评语等级优良

8、中及格不及格1.实验报告按时完成,字迹清楚,文字叙述流畅,逻辑性强2.实验方案设计合理3.实验过程（实验步骤详细,记录完整,数据合理,分析透彻）4实验结论正确. 成 绩： 指导教师签名： 批阅日期：附录1：源 程 序#include#include#define m 5#define n 16#define e 1.0e-10#define w 1.25int Jacobi(int ann,double bn,double xn);/Jacobi迭代法int SOR(int ann,double bn,double xn); /SOR迭代法函数int Gauss_Seidel(int ann,

9、double bn,double xn); /Gauss-Seidel迭代法函数double f(double a,double b);/边值函数void main() printf(*椭圆型微分方程数值解问题(Dirichlet边值条件)*nn); double xm+1,ym+1; double Um+1m+1; double Kn; double u;/精确解 double xx1n,xx2n,xx3n; int ann;/系数矩阵 double h; h=1.0/m; int i,j; for(i=0;i=m;i+) xi=i*h; yi=i*h; for(i=0;i=m;i+)/Di

10、richlet边界值 Ui0=f(xi,y0); Uim=f(xi,ym); U0i=f(x0,yi); Umi=f(xm,yi); for(i=0;in;i+) for(j=0;jn;j+) aij=0; Ki=0; for(i=0;in;i+)/系数矩阵 for(j=0;jn;j+) if(i=j) aij=4; else if(j-i=1&j%(m-1)!=0) aij=-1; aji=-1; else if(j-i=m-1) aij=-1; aji=-1; for(i=0;i0&in-m+1&in-1) Ki=Ui-(m-2)*(m-1)+1m; else if(i=n-1) Ki=U

11、m-1m+Umm-1; for(j=1;jm-2;j+) K(m-1)*j=U0j+1; K(m-1)*(j+1)-1=Umj+1; /* /测试结果正确性与否的已知数据如下(m=4,n=9)： K0=100; K1=20; K2=20; K3=80; K4=0; K5=0; K6=260; K7=180; K8=180; */ printf(向量K序列值如下：n); for(i=0;in;i+) printf(K%d=%-10.2lfn,i+1,Ki); printf(系数矩阵A如下：n); for(i=0;in;i+) for(j=0;jn;j+) printf(%-4d,aij); pr

12、intf(n); Jacobi(a,K,xx1);/进行Jacobi迭代 Gauss_Seidel(a,K,xx2);/进行Gauss_Seidel迭代 SOR(a,K,xx3);/进行SOR迭代 printf(迭代求的解为:n); printf(Ui 精确解 Jacobi Gauss_Seidel SORn); for(i=1;i=n;i+) /打印方程组的解 u=f(x(i-1)%(m-1)+1,y(i+m-2)/(m-1);/精确解 printf(U%-2d %-15.10lf %-15.10lf %-15.10lf %-15.10lfn,i,u,xx1i-1,xx2i-1,xx3i-1

13、); int Jacobi(int ann,double bn,double xn) int i,j; double Fan8;/定义向量无穷范数 double sum1,sum2;/定义累加器 int count=0;/计数器 double x0n,x1n; for(i=0;in;i+) x0i=1;/读入初值序列 for(i=0;in;i+) x1i=x0i;/初始化 do for(i=0;in;i+) x0i=x1i; for(i=1;i=n;i+) for(j=1,sum1=0;j=i-1;j+) sum1=sum1+ai-1j-1*x0j-1; for(j=i+1,sum2=0;j=

14、n;j+) sum2=sum2+ai-1j-1*x0j-1; x1i-1=(bi-1-sum1-sum2)/ai-1i-1;/Jacobi迭代公式 count+; for(Fan8=fabs(x10-x00),i=0;in;i+) if(Fan8e); for(i=0;in;i+) xi=x1i;/将结果返回 printf(Jacobi方法迭代次数为 %d次!nn,count); return 0; int SOR(int ann,double bn,double xn) int i,j; int count=0;/计数器 double Fan8;/定义向量无穷范数 double sum1,s

15、um2;/定义累加器 double x0n,x1n; /初值序列 for(i=0;in;i+) x0i=1; for(i=0;in;i+) x1i=x0i;/初始化 do for(i=0;in;i+) x0i=x1i; for(i=1;i=n;i+) for(j=1,sum1=0;j=i-1;j+) sum1=sum1+ai-1j-1*x1j-1; for(j=i,sum2=0;j=n;j+) sum2=sum2+ai-1j-1*x0j-1; x1i-1=x0i-1+w*(bi-1-sum1-sum2)/ai-1i-1;/SOR迭代公式 count+; for(Fan8=fabs(x10-x0

16、0),i=0;in;i+) if(Fan8e); for(i=0;in;i+) xi=x1i;/将结果返回 printf(SOR方法迭代次数为 %d次!nn,count); return 0;int Gauss_Seidel(int ann,double bn,double xn) int i,j; int count=0; double Fan8;/定义向量无穷范数 double sum1,sum2;/定义累加器 double x0n,x1n; /初值序列 for(i=0;in;i+) x0i=1; for(i=0;in;i+) x1i=x0i;/初始化 do for(i=0;in;i+)

17、x0i=x1i; for(i=1;i=n;i+) for(j=1,sum1=0;j=i-1;j+) sum1=sum1+ai-1j-1*x1j-1; for(j=i+1,sum2=0;j=n;j+) sum2=sum2+ai-1j-1*x0j-1; x1i-1=(bi-1-sum1-sum2)/ai-1i-1;/Gauss-Seidel迭代公式 count+; for(Fan8=fabs(x10-x00),i=0;in;i+) if(Fan8e); for(i=0;in;i+) xi=x1i;/将结果返回 printf(Gauss_Seidel方法迭代次数为%d次!nn,count); ret

18、urn 0;double f(double a,double b)/边值函数 double fun; fun=pow(a,2)-pow(b,2); return fun;附录2：实验报告填写说明 1实验项目名称：要求与实验教学大纲一致。2实验目的：目的要明确，要抓住重点，符合实验教学大纲要求。3实验原理：简要说明本实验项目所涉及的理论知识。4实验环境：实验用的软、硬件环境。5实验方案（思路、步骤和方法等）：这是实验报告极其重要的内容。概括整个实验过程。对于验证性实验，要写明依据何种原理、操作方法进行实验，要写明需要经过哪几个步骤来实现其操作。对于设计性和综合性实验，在上述内容基础上还应该画出流程图、设计思路和设计方法，再配以相应的文字说明。对于创新性实验，还应注明其创新点、特色。6实验过程（实验中涉及的记录、数据、分析）：写明具体实验方案的具体实施步骤，包括实验过程中的记录、数据和相应的分析。7实验结论（结果）：根据实验过程中得到的结果，做出结论。8实验小结：本次实验心得体会、思考和建议。9指导教师评语及成绩：指导教师依据学生的实际报告内容，给出本次实验报告的评价。