五点差分格式

微分方程数值解实验报告 姓名 丁建伟 学号 200708020211 日期

2010.12.30 实验项目 五点差分格式 指导教师

徐强 一、上机实验的问题和要求（需求分析）：

实验内容：

(I) 分别在步长h=1/16,1/32,1/64,1/128情形下用五点差分格式计算?

??Ω?=Ω=-?on u in f u ;0; 其中Ω=（0，1）×（0，1）； 2=f πsin πxsin πy ；精确解为 u=sin πxsin πy ；

(I)给出系数矩阵A 的图像。

(II) 用G-S 迭代法求解AV=F ，并分析无穷范数和L2范数下的误差阶。

二、程序设计的基本思想，原理和算法描述：

基本思想及原理：

三、主要程序代码或命令：

#include

#include

#define MAX 200

#define M 3.1415926

void main()

{

int n,i,j;

float V[MAX][MAX]={{0}};

float F[MAX][MAX],U[MAX][MAX],e[MAX][MAX];

float h,x,y,err,temp;

printf("请输入分割数n的值:");

scanf("%d",&n);

h=1/float(n+1);

for(i=1,x=h;i<=n;i++,x+=h)

for(j=1,y=h;j<=n;j++,y+=h)

{F[i][j]=2*M*M*sin(M*x)*sin(M*y);

U[i][j]=sin(M*x)*sin(M*y);}

err=1.0;

while(err>1.0e-6)

{err=0.0;

for(j=1;j<=n;j++)

for(i=1;i<=n;i++)

{temp=0.25*(V[i+1][j]+V[i-1][j]+V[i][j-1]+V[i][j+1])+(h*h/4)*F[i][j]; err=(temp-V[i][j])*(temp-V[i][j])+err;

V[i][j]=temp;

}

err=sqrt(err);

}

float e1=0,e2=0;

for(i=1;i<=n;i++)

for(j=1;j<=n;j++)

{

e[i][j]=fabs(V[i][j]-U[i][j]);

e2+=h*e[i][j]*e[i][j];

if(e1

e1=e[i][j];

}

e2=sqrt(e2);

printf("请输出在无穷范数和二范数下的误差:e1=%f和e2=%f\n",e1,e2);

for(i=1;i<=n;i++)

for(j=1;j<=n;j++)

printf("请输出V[%d][%d]的值:%f\n",i,j,V[i][j]);

}

四、调试和运行程序过程中产生的问题及采取的措施：

1.编译时出错，若想在主函数和被调用函数都使用一些变量，必须把这些变量设为全局变量。

2.编译时，没有注意数据类型转换，如float h;h=1/n;是错误的，因为n是整形的，当n 值大于1时，h老为零。应进行模式转换，h=1/float（n）；这样才是正确的。

3.对浮点数求绝对值时，应使用fabs（）函数，而不是abs（）。

4.注意乘除的计算，不能直接写成2x等，必须加上符号，2*x。

五、运行输出结果及分析：

上述程序在Visual C++ 6.0环境下加以实现。经过多次测试，程序运行正确。例如：分别输入n值：20 ，40，运行结果如图所示，图中显示了每一步的值及端点误差。

请输入等分数n值：20

输入n值：40

由图可知：

a.四种步长下无穷范数分别为：h=1/20时为e1=0.013742，h=1/40时为e2=0.003477，h=1/80时为e3=0.000872，h=1/160时为e4=0.000218。计算可得差分解近似精确解在无穷范数下误差阶数为二阶。(e=max{u[i]})

b.四种步长下L2范数分别为：h=1/20时为e1'=0.009277，h=1/40时为e2'=0.002342，h=1/80时为e3'=0.000586，h=1/160时为e4'=0.000145。计算可得差分解近似精确解在L2范数下误差阶数为二阶。(2/111)('2∑-==n i i u

e )

通过这次课程设计:

1. 我又进一步巩固了C 语言的基础。

2. 做课程设计达到了理论与实践结合的目的，提高了自己的编程能力。

3. 对数值分析里求解线性方程组的追赶法复习了一下，加深了理解。

4. 掌握中心差分格式的程序实现和分析算法误差的方法。