椭圆型微分方程.docx
- 文档编号:6117210
- 上传时间:2023-01-03
- 格式:DOCX
- 页数:17
- 大小:120.66KB
椭圆型微分方程.docx
《椭圆型微分方程.docx》由会员分享,可在线阅读,更多相关《椭圆型微分方程.docx(17页珍藏版)》请在冰豆网上搜索。
椭圆型微分方程
数学与计算科学学院
实验报告
实验项目名称椭圆型方程数值解
所属课程名称微分方程数值解法
实验类型验证
实验日期
班级信计0902
学号
姓名
成绩
一、实验概述:
【实验目的】
1掌握椭圆型方程的五点差分方法
2掌握求解线性方程组的Jacobi迭代法
【实验原理】
【实验环境】
MicrosoftVisualC++6.0
二、实验内容:
【实验方案】
用五点差分格式的Jacobi迭代法求解单位正方形区域
上的
方程
边值问题
【实验过程】(实验步骤、记录、数据、分析)
1,部分主要算法代码:
(1)Dirichlet边值、系数矩阵A以及向量序列K的初步赋值
for(i=0;i<=m;i++)//Dirichlet边界值
{
U[i][0]=f(x[i],y[0]);
U[i][m]=f(x[i],y[m]);
U[0][i]=f(x[0],y[i]);
U[m][i]=f(x[m],y[i]);
}
for(i=0;i for(j=0;j { a[i][j]=0; K[i]=0; } for(i=0;i for(j=0;j { if(i==j)a[i][j]=4; elseif(j-i==1&&j%(m-1)! =0) { a[i][j]=-1; a[j][i]=-1; } elseif(j-i==m-1) { a[i][j]=-1; a[j][i]=-1; } } for(i=0;i { if(i==0)K[0]=U[1][0]+U[0][1]; elseif(i>0&&i elseif(i==m-2)K[i]=U[i+1][0]+U[i+2][1]; elseif(i==n-m+1)K[i]=U[0][m-1]+U[1][m]; elseif(i>n-m+1&&i elseif(i==n-1)K[i]=U[m-1][m]+U[m][m-1]; } for(j=1;j { K[(m-1)*j]=U[0][j+1]; K[(m-1)*(j+1)-1]=U[m][j+1]; } (2)Jacobi迭代法 intJacobi(inta[n][n],doubleb[n],doublex[n]) { inti,j; doubleFan8;//定义向量无穷范数 doublesum1,sum2;//定义累加器 intcount=0;//计数器 doublex0[n],x1[n]; for(i=0;i x0[i]=1;//读入初值序列 for(i=0;i x1[i]=x0[i];//初始化 do { for(i=0;i x0[i]=x1[i]; for(i=1;i<=n;i++) { for(j=1,sum1=0;j<=i-1;j++) sum1=sum1+a[i-1][j-1]*x0[j-1]; for(j=i+1,sum2=0;j<=n;j++) sum2=sum2+a[i-1][j-1]*x0[j-1]; x1[i-1]=(b[i-1]-sum1-sum2)/a[i-1][i-1];//Jacobi迭代公式 } count++; for(Fan8=fabs(x1[0]-x0[0]),i=0;i if(Fan8 Fan8=fabs(x1[i]-x0[i]);//进行无穷范数运算 }while(Fan8>e); for(i=0;i x[i]=x1[i];//将结果返回 printf("Jacobi方法迭代次数为%d次! \n\n",count); return0; } (3)Gauss_seidel迭代法 仅将Jacobi迭代代码的sum1=sum1+a[i-1][j-1]*x0[j-1]处改为sum1=sum1+a[i-1][j-1]*x1[j-1]; (4)SOR迭代法 intSOR(inta[n][n],doubleb[n],doublex[n]) { inti,j; intcount=0;//计数器 doubleFan8;//定义向量无穷范数 doublesum1,sum2;//定义累加器 doublex0[n],x1[n]; //初值序列 for(i=0;i x0[i]=1; for(i=0;i x1[i]=x0[i];//初始化 do { for(i=0;i x0[i]=x1[i]; for(i=1;i<=n;i++) { for(j=1,sum1=0;j<=i-1;j++) sum1=sum1+a[i-1][j-1]*x1[j-1]; for(j=i,sum2=0;j<=n;j++) sum2=sum2+a[i-1][j-1]*x0[j-1]; x1[i-1]=x0[i-1]+w*(b[i-1]-sum1-sum2)/a[i-1][i-1];//SOR迭代公式 } count++; for(Fan8=fabs(x1[0]-x0[0]),i=0;i if(Fan8 Fan8=fabs(x1[i]-x0[i]);//进行无穷范数运算 }while(Fan8>e); for(i=0;i x[i]=x1[i];//将结果返回 printf("SOR方法迭代次数为%d次! \n\n",count); return0; } (5)边值函数(解析解) doublef(doublea,doubleb)//边值函数 { doublefun; fun=pow(a,2)-pow(b,2); returnfun; } 2,结果分析: 1,A为对称占优矩阵,且为对称正定矩阵,存在唯一解。 2,迭代速度从低到高依次为: Jacobi迭代法、Gause_Seidel迭代法,SOR迭代法。 【实验结论】(结果) 【实验小结】(收获体会) 通过本次实验,我重新复习了Jacobi,Gause_Seidel以及SOR迭代法,也对课本上的正方形区域中的Laplace方程Dirichlet边值问题有了更深的了解,同时也提高了自己的上机编程能力,感受到了理论与实践相结合的乐趣。 三、指导教师评语及成绩: 评语 评语等级 优 良 中 及格 不及格 1.实验报告按时完成,字迹清楚,文字叙述流畅,逻辑性强 2.实验方案设计合理 3.实验过程(实验步骤详细,记录完整,数据合理,分析透彻) 4实验结论正确. 成绩: 指导教师签名: 批阅日期: 附录1: 源程序 #include #include #definem5 #definen16 #definee1.0e-10 #definew1.25 intJacobi(inta[n][n],doubleb[n],doublex[n]);//Jacobi迭代法 intSOR(inta[n][n],doubleb[n],doublex[n]);//SOR迭代法函数 intGauss_Seidel(inta[n][n],doubleb[n],doublex[n]);//Gauss-Seidel迭代法函数 doublef(doublea,doubleb);//边值函数 voidmain() { printf("***椭圆型微分方程数值解问题(Dirichlet边值条件)***\n\n"); doublex[m+1],y[m+1]; doubleU[m+1][m+1]; doubleK[n]; doubleu;//精确解 doublexx1[n],xx2[n],xx3[n]; inta[n][n];//系数矩阵 doubleh; h=1.0/m; inti,j; for(i=0;i<=m;i++) { x[i]=i*h; y[i]=i*h; } for(i=0;i<=m;i++)//Dirichlet边界值 { U[i][0]=f(x[i],y[0]); U[i][m]=f(x[i],y[m]); U[0][i]=f(x[0],y[i]); U[m][i]=f(x[m],y[i]); } for(i=0;i for(j=0;j { a[i][j]=0; K[i]=0; } for(i=0;i for(j=0;j { if(i==j)a[i][j]=4; elseif(j-i==1&&j%(m-1)! =0) { a[i][j]=-1; a[j][i]=-1; } elseif(j-i==m-1) { a[i][j]=-1; a[j][i]=-1; } } for(i=0;i { if(i==0)K[0]=U[1][0]+U[0][1]; elseif(i>0&&i elseif(i==m-2)K[i]=U[i+1][0]+U[i+2][1]; elseif(i==n-m+1)K[i]=U[0][m-1]+U[1][m]; elseif(i>n-m+1&&i elseif(i==n-1)K[i]=U[m-1][m]+U[m][m-1]; } for(j=1;j { K[(m-1)*j]=U[0][j+1]; K[(m-1)*(j+1)-1]=U[m][j+1]; } /* //测试结果正确性与否的已知数据如下(m=4,n=9): K[0]=100; K[1]=20; K[2]=20; K[3]=80; K[4]=0; K[5]=0; K[6]=260; K[7]=180; K[8]=180; */ printf("向量K序列值如下: \n"); for(i=0;i printf("K[%d]=%-10.2lf\n",i+1,K[i]); printf("系数矩阵A如下: \n"); for(i=0;i { for(j=0;j { printf("%-4d",a[i][j]); } printf("\n"); } Jacobi(a,K,xx1);//进行Jacobi迭代 Gauss_Seidel(a,K,xx2);//进行Gauss_Seidel迭代 SOR(a,K,xx3);//进行SOR迭代 printf("迭代求的解为: \n"); printf("U[i]精确解JacobiGauss_SeidelSOR\n"); 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.10lf\n",i,u,xx1[i-1],xx2[i-1],xx3[i-1]); } } intJacobi(inta[n][n],doubleb[n],doublex[n]) { inti,j; doubleFan8;//定义向量无穷范数 doublesum1,sum2;//定义累加器 intcount=0;//计数器 doublex0[n],x1[n]; for(i=0;i x0[i]=1;//读入初值序列 for(i=0;i x1[i]=x0[i];//初始化 do { for(i=0;i x0[i]=x1[i]; for(i=1;i<=n;i++) { for(j=1,sum1=0;j<=i-1;j++) sum1=sum1+a[i-1][j-1]*x0[j-1]; for(j=i+1,sum2=0;j<=n;j++) sum2=sum2+a[i-1][j-1]*x0[j-1]; x1[i-1]=(b[i-1]-sum1-sum2)/a[i-1][i-1];//Jacobi迭代公式 } count++; for(Fan8=fabs(x1[0]-x0[0]),i=0;i if(Fan8 Fan8=fabs(x1[i]-x0[i]);//进行无穷范数运算 }while(Fan8>e); for(i=0;i x[i]=x1[i];//将结果返回 printf("Jacobi方法迭代次数为%d次! \n\n",count); return0; } intSOR(inta[n][n],doubleb[n],doublex[n]) { inti,j; intcount=0;//计数器 doubleFan8;//定义向量无穷范数 doublesum1,sum2;//定义累加器 doublex0[n],x1[n]; //初值序列 for(i=0;i x0[i]=1; for(i=0;i x1[i]=x0[i];//初始化 do { for(i=0;i x0[i]=x1[i]; for(i=1;i<=n;i++) { for(j=1,sum1=0;j<=i-1;j++) sum1=sum1+a[i-1][j-1]*x1[j-1]; for(j=i,sum2=0;j<=n;j++) sum2=sum2+a[i-1][j-1]*x0[j-1]; x1[i-1]=x0[i-1]+w*(b[i-1]-sum1-sum2)/a[i-1][i-1];//SOR迭代公式 } count++; for(Fan8=fabs(x1[0]-x0[0]),i=0;i if(Fan8 Fan8=fabs(x1[i]-x0[i]);//进行无穷范数运算 }while(Fan8>e); for(i=0;i x[i]=x1[i];//将结果返回 printf("SOR方法迭代次数为%d次! \n\n",count); return0; } intGauss_Seidel(inta[n][n],doubleb[n],doublex[n]) { inti,j; intcount=0; doubleFan8;//定义向量无穷范数 doublesum1,sum2;//定义累加器 doublex0[n],x1[n]; //初值序列 for(i=0;i x0[i]=1; for(i=0;i x1[i]=x0[i];//初始化 do { for(i=0;i x0[i]=x1[i]; for(i=1;i<=n;i++) { for(j=1,sum1=0;j<=i-1;j++) sum1=sum1+a[i-1][j-1]*x1[j-1]; for(j=i+1,sum2=0;j<=n;j++) sum2=sum2+a[i-1][j-1]*x0[j-1]; x1[i-1]=(b[i-1]-sum1-sum2)/a[i-1][i-1];//Gauss-Seidel迭代公式 } count++; for(Fan8=fabs(x1[0]-x0[0]),i=0;i if(Fan8 Fan8=fabs(x1[i]-x0[i]);//进行无穷范数运算 }while(Fan8>e); for(i=0;i x[i]=x1[i];//将结果返回 printf("Gauss_Seidel方法迭代次数为%d次! \n\n",count); return0; } doublef(doublea,doubleb)//边值函数 { doublefun; fun=pow(a,2)-pow(b,2); returnfun; } 附录2: 实验报告填写说明 1.实验项目名称: 要求与实验教学大纲一致。 2.实验目的: 目的要明确,要抓住重点,符合实验教学大纲要求。 3.实验原理: 简要说明本实验项目所涉及的理论知识。 4.实验环境: 实验用的软、硬件环境。 5.实验方案(思路、步骤和方法等): 这是实验报告极其重要的内容。 概括整个实验过程。 对于验证性实验,要写明依据何种原理、操作方法进行实验,要写明需要经过哪几个步骤来实现其操作。 对于设计性和综合性实验,在上述内容基础上还应该画出流程图、设计思路和设计方法,再配以相应的文字说明。 对于创新性实验,还应注明其创新点、特色。 6.实验过程(实验中涉及的记录、数据、分析): 写明具体实验方案的具体实施步骤,包括实验过程中的记录、数据和相应的分析。 7.实验结论(结果): 根据实验过程中得到的结果,做出结论。 8.实验小结: 本次实验心得体会、思考和建议。 9.指导教师评语及成绩: 指导教师依据学生的实际报告内容,给出本次实验报告的评价。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 椭圆 微分方程
![提示](https://static.bdocx.com/images/bang_tan.gif)