椭圆型微分方程.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;ifor（j=0;j{
a[i][j]=0;
K[i]=0;
}
for（i=0;ifor（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&&ielseif（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&&ielseif（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;ix0[i]=1;//读入初值序列
for（i=0;ix1[i]=x0[i];//初始化
do
{
for（i=0;ix0[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;iif（Fan8Fan8=fabs（x1[i]-x0[i]）;//进行无穷范数运算
}while（Fan8>e）;
for（i=0;ix[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;ix0[i]=1;
for（i=0;ix1[i]=x0[i];//初始化
do
{
for（i=0;ix0[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;iif（Fan8Fan8=fabs（x1[i]-x0[i]）;//进行无穷范数运算
}while（Fan8>e）;
for（i=0;ix[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;ifor（j=0;j{
a[i][j]=0;
K[i]=0;
}
for（i=0;ifor（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&&ielseif（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&&ielseif（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;iprintf（"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;ix0[i]=1;//读入初值序列
for（i=0;ix1[i]=x0[i];//初始化
do
{
for（i=0;ix0[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;iif（Fan8Fan8=fabs（x1[i]-x0[i]）;//进行无穷范数运算
}while（Fan8>e）;
for（i=0;ix[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;ix0[i]=1;
for（i=0;ix1[i]=x0[i];//初始化
do
{
for（i=0;ix0[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;iif（Fan8Fan8=fabs（x1[i]-x0[i]）;//进行无穷范数运算
}while（Fan8>e）;
for（i=0;ix[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;ix0[i]=1;
for（i=0;ix1[i]=x0[i];//初始化
do
{
for（i=0;ix0[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;iif（Fan8Fan8=fabs（x1[i]-x0[i]）;//进行无穷范数运算
}while（Fan8>e）;
for（i=0;ix[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．指导教师评语及成绩：
指导教师依据学生的实际报告内容，给出本次实验报告的评价。
升级会员 