数值线性代数第二版徐树方高立张平文上机习题第四章实验报告.docx
《数值线性代数第二版徐树方高立张平文上机习题第四章实验报告.docx》由会员分享，可在线阅读，更多相关《数值线性代数第二版徐树方高立张平文上机习题第四章实验报告.docx（10页珍藏版）》请在冰豆网上搜索。
数值线性代数第二版徐树方高立张平文上机习题第四章实验报告
第四章上机习题
1考虑两点边值问题
容易知道它的精确解为
为了把微分方程离散化，把[0,1]区间n等分，令h=1/n，
得到差分方程
简化为
从而离散化后得到的线性方程组的系数矩阵为
对
分别用Jacobi迭代法，G-S迭代法和SOR迭代法求线性方程组的解，要求有4位有效数字，然后比较与精确解得误差。
对
考虑同样的问题。
解
（1）给出算法：
为解
，令
，其中
，
利用Jacobi迭代法，G-S迭代法，SOR迭代法解线性方程组，均可以下步骤求解：
step1给定初始向量x0=（0,0,...,0），最大迭代次数N，精度要求c，令k=1
step2令x=B*x0+g
step3若||x-x0||2step4若k>=N,算法停止，迭代失败，否则，令x0=x，转step2
在Jacobi迭代法中，B=D-1*（L+U）,g=D-1*b
在G-S迭代法中，B=D-1*（L+U）,g=D-1*b
在SOR迭代法中，B=（D-w*L）-1*[（1-w）*D+w*U],g=w*（D-w*L）-1*b
另外，在SOR迭代法中，上面算法step1中要给定松弛因子w，其中0为计算结果，规定w=。
（2）给定方程
在Ax=b中，矩阵A如题目所示，且为n-1阶矩阵
由于在第1、n-1个方程中，由于y（0）、y
（1）的存在，使方程右端应减去y（0）、y
（1）项，即b
（1）=a*h2-c*y（0），b（n-1）=a*h2-（c+h）*y（n），b（i）=a*h2,i=2,...,n-1
程序为
1Jacobi迭代法编成的函数[x,k]=Jacobi（A,b,c,N）
function[x,k]=Jacobi（A,b,c,N）
D=diag（diag（A））;
L=triu（A）-A;
U=tril（A）-A;
B=D^（-1）*（L+U）;
g=D^（-1）*b;
x0=zeros（length（A）,1）;
x=B*x0+g;
k=1;
whilenorm（x-x0,2）>=
x0=x;
x=B*x0+g;
k=k+1;
ifk>=N
break
end
end
end
2G-S迭代法编成的函数[x,k]=GaussSeidel（A,b,c,N）
function[x,k]=GaussSeidel（A,b,c,N）
U=diag（diag（A））-triu（A）;
x0=zeros（length（A）,1）;
B=tril（A）^（-1）*U;
g=tril（A）^（-1）*b;
x=B*x0+g;
k=1;
whilenorm（x-x0,2）>=
x0=x;
x=B*x0+g;
k=k+1;
ifk>=N
break
end
end
end
3SOR迭代法编成的函数[x,k]=SOR（A,b,w,c,N）
function[x,k]=SOR（A,b,w,c,N）
D=diag（diag（A））;
L=D-tril（A）;
U=D-triu（A）;
x0=zeros（length（A）,1）;
B=（D-w*L）^（-1）*（（1-w）*D+w*U）;
g=w*（D-w*L）^（-1）*b;
x=B*x0+g;
k=1;
whilenorm（x-x0,2）>=
x0=x;
x=B*x0+g;
k=k+1;
ifk>=N
break
end
end
end
4问题1求解ex4_1
clear;clc;
%c=1;
%c=
%c=;
c=;
a=1/2;n=100;h=1/n;
w=1/2;N=1000000;
A=-（2*c+h）*eye（n-1）;
fori=2:
n-1w
A（i-1,i）=c+h;
A（i,i-1）=c;
end
b=[a*h^2*ones（n-2,1）;a*h^2-（c+h）];
fori=1:
n-1
x（i）=i*h;
y（i）=（（1-a）/（1-exp（-1/c）））*（1-exp（-x（i）/c））+a*x（i）;
end
[y1,n1]=Jacobi（A,b,c,N）;
[y2,n2]=GaussSeidel（A,b,c,N）;
[y3,n3]=SOR（A,b,w,c,N）;
disp（['c=',num2str（c）,'时']）;
disp（['Jacobi迭代与精确解的差为',num2str（norm（y'-y1,inf））]）;
disp（['迭代次数为',num2str（n1）]）;
disp（['G-S迭代与精确解的差为',num2str（norm（y'-y2,inf））]）;
disp（['迭代次数为',num2str（n2）]）;
disp（['SOR迭代与精确解的差为',num2str（norm（y'-y3,inf））]）;
disp（['迭代次数为',num2str（n3）]）;
计算结果为
（1）
c=1时
Jacobi迭代与精确解的差为
迭代次数为11796
G-S迭代与精确解的差为
迭代次数为6227
SOR迭代与精确解的差为
迭代次数为15367
（2）
c=时
Jacobi迭代与精确解的差为
迭代次数为5353
G-S迭代与精确解的差为
迭代次数为2797
SOR迭代与精确解的差为
迭代次数为7300
（3）
c=时
Jacobi迭代与精确解的差为
迭代次数为532
G-S迭代与精确解的差为
迭代次数为318
SOR迭代与精确解的差为
迭代次数为834
（4）
c=时
Jacobi迭代与精确解的差为
迭代次数为116
G-S迭代与精确解的差为
迭代次数为108
SOR迭代与精确解的差为
迭代次数为267
结果分析
三种迭代法的误差基本相同，且G-S迭代法的收敛速度明显小于Jacobi迭代法，但SOR迭代法收敛速度较慢，原因是收敛因子非最佳。
2考虑偏微分方程
其中边界条件为u=1.沿x方向和y方向均匀剖分为N等份，令h=1/N，并设应用中心差分离散化后得到差分方程的代数方程组为
取g（x,y）和f（x,y）分别为exp（xy）和x+y，用G-S迭代法求解上述方程组，并请列表比较N=20,40,80时收敛所需要的迭代次数和所用的CPU时间。
迭代终止条件为||xk+1-xk||2<10-7.
解求解过程同问题1
给定方程，将
组成（n-1）2维列向量，记为
则对Ax=b，由于g（x,y）和f（x,y）分别为exp（xy）和x+y，可得A的形式为
由于边值
，可得，当
中的下标存在0,n时，h2（i*h+j*h）应相应的加上1或2（当有两个变量下标均含有0,n时）.
程序为
（1）Jacobi迭代法的程序同问题1
（2）问题2求解程序为ex4_2
clear;clc;
c=10^（-7）;
n=[20,40,80];
N=1000000;
form=1:
3
h=1/n（m）;
A=zeros（（n（m）-1）^2）;
b=zeros（（n（m）-1）^2,1）;
fori=1:
（n（m）-1）^2
ifi>1
A（i-1,i）=-1;A（i,i-1）=-1;
end
ifi>n（m）-1
A（i,i-n（m）+1）=-1;A（i-n（m）+1,i）=-1;
end
ii=ceil（i/（n（m）-1））;
ifmod（i,n（m）-1）~=0
jj=mod（i,n（m）-1）;
else
jj=n（m）-1;
end
A（i,i）=4+exp（ii*jj*h^2）;
b（i）=h^3*（ii+jj）;
ifii==1||ii==n（m）-1
b（i）=b（i）+1;
end
ifjj==1||jj==n（m）-1
b（i）=b（i）+1;
end
end
disp（['n=',num2str（n（m））]）
tic
[y,k]=GaussSeidel（A,b,c,N）;
toc
disp（['迭代次数为',num2str（k）]）;
end
结果为
n
20
40
80
CPUtime
seconds
seconds
seconds
迭代次数
24
26
27
升级会员 