微分方程数值解法实验2.docx
《微分方程数值解法实验2.docx》由会员分享，可在线阅读，更多相关《微分方程数值解法实验2.docx（14页珍藏版）》请在冰豆网上搜索。
微分方程数值解法实验2
一、实验概述：
【实验目的】
熟练掌握求解
方程第一边值问题的五点差分格式。
【实验原理】
1.取定沿X轴和Y轴的步长h1和h2，
作两族与坐标轴平行的直线：
，两族直线的交点称为网点或节点，记为
，假定
为正则内点。
沿x,y方向分别用二阶中心差商代替
则得
，由于差分方程中只出现u在（i,j）及其四个邻点上的值，故称为五点差分格式。
2.
迭代法
设线性方程组
（1）
的系数矩阵
可逆且主对角元素
均不为零，令
，并将
分解成
（2）
从而
（1）可以写成
，令
，其中
（3）
以
为迭代矩阵的迭代法
（4）
称为
迭代法，用向量的分量来表示，（4）为
（5）
其中
为初始向量。
3.
迭代法
把矩阵
分解成
（6）
其中
，
分别为
的主对角元除外的下三角和上三角部分，于是方程组
（1）可以写成
，即
，其中
（7）
以
为迭代矩阵构成的迭代法
（8）
称为
迭代法，用向量表示的形式为
4.
迭代
【实验环境】
MATLABR2014a
二、实验内容：
【实验方案】
1.求解边值问题
取步长h=k=1/64,1/128,做五点差分格式。
用jacobi迭代，Gauss-Seidel迭代和SOR迭代（取
）求解差分方程，以前后两次重合到小数点后四位的迭代值作为解的近似。
比较三种解法的迭代次数以及差分解
与精确解
的精度。
2.在
软件上编写代码；
3.运行得出结论。
【实验过程】（实验步骤、记录、数据、分析）
1.分析问题。
2.在
上编写代码，调试。
3.运行代码，发现问题后改进代码。
4.分别输入n值为64和128，得出结论。
【实验结论】（结果）
n=64时:
1.Jacobi迭代
迭代次数
i=7811
2.Guass-Seidel迭代
迭代次数
i=4211
3.SOR超松弛法
迭代次数
i=121
n=128时：
1.Jacobi迭代
迭代次数
i=20001
2.Guass-Seidel迭代
迭代次数
i=15660
3.SOR超松弛法
迭代次数
i=242
【实验小结】（收获体会）
通过本次实验，掌握了求解
方程第一边值问题的五点差分格式，比较了三种解法的迭代次数以及差分解与精确解的精度，更加熟悉的掌握了利用MATLAB求解数学问题的方法。
三、指导教师评语及成绩：
评语
评语等级
优
良
中
及格
不及格
1.实验报告按时完成,字迹清楚,文字叙述流畅,逻辑性强
2.实验方案设计合理
3.实验过程（实验步骤详细,记录完整,数据合理,分析透彻）
4实验结论正确.
成绩：
指导教师签名：
批阅日期：
附录1：
源程序
clc;
n=input（'请输入等分数n='）;
h=1/n;
A=zeros
（1）;
A=sparse（A）;
b=zeros（（n-1）^2,1）;
b=sparse（b）;
u=zeros（（n-1）^2,1）;
u=sparse（u）;
I=speye（（n-1）^2）;
D=zeros（（n-1）^2,1）;
D=sparse（D）;
fori=1:
（n-1）
forj=1:
（n-1）
xi=i*h;
yj=j*h;
b（j+（i-1）*（n-1）,1）=2*pi^2*exp（pi*（xi+yj））*（sin（pi*xi）*cos（pi*yj）+cos（pi*xi）*sin（pi*yj））;
u（j+（i-1）*（n-1）,1）=exp（pi*（xi+yj））*sin（pi*xi）*sin（pi*yj）;%精确解
end
end
%---------------------建立系数矩阵A--------------------------
j=0;
fori=1:
（n-1）^2
if（i>=n）
j=j+1;
A（i,j）=1/h^2;
A（j,i）=1/h^2;
end
end
fori=1:
（n-1）^2
A（i,i）=-4/h^2;
end
j=0;
fori=2:
（n-1）^2
j=j+1;
if（mod（j,n-1）~=0）
A（i,j）=1/h^2;
A（j,i）=1/h^2;
end
end
A;
D=diag（diag（A））;
L=tril（A,-1）;
R=triu（A,1）;
L=sparse（L）;
R=sparse（R）;
Z=speye
（1）;
x=1/n:
1/n:
（n-1）/n;
y=1/n:
1/n:
（n-1）/n;
[X,Y]=meshgrid（x,y）;
Z1=exp（pi*（X+Y））.*sin（pi*X）.*sin（pi*Y）;
subplot（2,2,1）;
mesh（X,Y,Z1）;%精确解图像
title（'精确解图像'）;
xlabel（'x'）;
ylabel（'y'）;
zlabel（'u'）;
%------------Jacboi迭代使用------------------
G=-inv（D）*（L+R）;
H=inv（D）*b;
%------------Guass-Seidel迭代使用----------------
M=inv（D+L）;
M=sparse（M）;
R1=-M*R;
R1=sparse（R1）;
M=M*b;
%-----------SOR超松弛法使用--------------------
Aw=zeros
（1）;
Aw=sparse（Aw）;
B=zeros
（1）;
B=sparse（B）;
B=I-inv（D）*A;
q=eig（B）;
q=q
（1）;
w=2/（1+sqrt（1-q^2））;
P=inv（D-w*（-L））;
Aw=P*（w*（-R）+（1-w）*D）;
P=w*P*b;
%-----------------------------------------------------
fork=1:
3
switchk
%--------------------------Jacobi迭代--------------
case1
U=speye（（n-1）^2,1）;
U0=zeros（（n-1）^2,1）;
U0=sparse（U0）;
disp（'1.Jacobi迭代'）;
i=0;
whilenorm（U-U0）>1e-4
U0=U;
U=G*U0+H;
i=i+1;
if（i>20000）
break
end
end
disp（'迭代次数'）;
i
fora=1:
n-1
forb=1:
n-1
Z（a,b）=U（（n-1）*（a-1）+b,1）;
end
end
subplot（2,2,2）;
mesh（X,Y,Z）;
title（'Jacobi迭代图像'）;
xlabel（'x'）;
ylabel（'y'）;
zlabel（'u'）;
%----------------Guass-Seidel迭代------------------
case2
U=speye（（n-1）^2,1）;
U0=zeros（（n-1）^2,1）;
U0=sparse（U0）;
disp（'2.Guass-Seidel迭代'）;
i=0;
whilenorm（U-U0）>1e-4
U0=U;
U=R1*U0+M;
i=i+1;
if（i>20000）
break
end
end
disp（'迭代次数'）;
i
fora=1:
n-1
forb=1:
n-1
Z（a,b）=U（（n-1）*（a-1）+b,1）;
end
end
subplot（2,2,3）;
mesh（X,Y,Z）;
title（'Guass-Seidel迭代图像'）;
xlabel（'x'）;
ylabel（'y'）;
zlabel（'u'）;
%---------------SOR超松弛法------------------------
case3
U=speye（（n-1）^2,1）;
U0=zeros（（n-1）^2,1）;
U0=sparse（U0）;
disp（'3.SOR超松弛法'）;
i=0;
whilenorm（U-U0）>1e-1
U0=U;
U=Aw*U0+P;
i=i+1;
if（i>20000）
break
end
end
disp（'迭代次数'）;
i
fora=1:
n-1
forb=1:
n-1
Z（a,b）=U（（n-1）*（a-1）+b,1）;
end
end
subplot（2,2,4）;
mesh（X,Y,Z）;
title（'SOR超松弛法图像'）;
xlabel（'x'）;
ylabel（'y'）;
zlabel（'u'）;
end
end
附录2：
实验报告填写说明
1．实验项目名称：
要求与实验教学大纲一致。
2．实验目的：
目的要明确，要抓住重点，符合实验教学大纲要求。
3．实验原理：
简要说明本实验项目所涉及的理论知识。
4．实验环境：
实验用的软、硬件环境。
5．实验方案（思路、步骤和方法等）：
这是实验报告极其重要的内容。
概括整个实验过程。
对于验证性实验，要写明依据何种原理、操作方法进行实验，要写明需要经过哪几个步骤来实现其操作。
对于设计性和综合性实验，在上述内容基础上还应该画出流程图、设计思路和设计方法，再配以相应的文字说明。
对于创新性实验，还应注明其创新点、特色。
6．实验过程（实验中涉及的记录、数据、分析）：
写明具体实验方案的具体实施步骤，包括实验过程中的记录、数据和相应的分析。
7．实验结论（结果）：
根据实验过程中得到的结果，做出结论。
8．实验小结：
本次实验心得体会、思考和建议。
9．指导教师评语及成绩：
指导教师依据学生的实际报告内容，给出本次实验报告的评价。
升级会员 