数值线性代数第二版上机习题第四章实验报告.doc

<p>第四章上机习题
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||2</p><c,算法停止,输出解与迭代次数k,否则,转step4 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 1="" w="0、5。" ax="b中,矩阵A如题目所示,且为n-1阶矩阵" i="2,、、、,n-1" function="" 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;" while="">=0、00001
&nbsp; &nbsp;x0=x;
&nbsp; &nbsp;x=B*x0+g;
&nbsp; &nbsp;k=k+1;
&nbsp; &nbsp;if k&gt;=N
&nbsp; &nbsp; &nbsp; &nbsp;break
&nbsp; &nbsp;end
end
end
2 G-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;
while norm(x-x0,2)&gt;=0、00001
&nbsp; &nbsp;x0=x;
&nbsp; &nbsp;x=B*x0+g;
&nbsp; &nbsp;k=k+1;
&nbsp; &nbsp;if k&gt;=N
&nbsp; &nbsp; &nbsp; &nbsp;break
&nbsp; &nbsp;end
end
end
3 &nbsp;SOR迭代法编成得函数 [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;
while norm(x-x0,2)&gt;=0、00001
&nbsp; &nbsp;x0=x;
&nbsp; &nbsp;x=B*x0+g;
&nbsp; &nbsp;k=k+1;
&nbsp; &nbsp;if k&gt;=N
&nbsp; &nbsp; &nbsp; &nbsp;break
&nbsp; &nbsp;end
end
end
4 问题1求解 ex4_1
clear;clc;
%c=1;
%c=0、1
%c=0、01;
c=0、0001;
a=1/2;n=100;h=1/n;
w=1/2;N=1000000;
A=-(2*c+h)*eye(n-1);
for i=2:n-1w
&nbsp; &nbsp;A(i-1,i)=c+h;
&nbsp; &nbsp;A(i,i-1)=c;
end
b=[a*h^2*ones(n-2,1);a*h^2-(c+h)];
for i=1:n-1
&nbsp; &nbsp;x(i)=i*h;
&nbsp; &nbsp;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([&#39;c=&#39;,num2str(c),&#39;时&#39;]);
disp([&#39;Jacobi迭代与精确解得差为&#39;,num2str(norm(y&#39;-y1,inf))]);
disp([&#39;迭代次数为&#39;,num2str(n1)]);
disp([&#39;G-S迭代与精确解得差为&#39;,num2str(norm(y&#39;-y2,inf))]);
disp([&#39;迭代次数为&#39;,num2str(n2)]);
disp([&#39;SOR迭代与精确解得差为&#39;,num2str(norm(y&#39;-y3,inf))]);
disp([&#39;迭代次数为&#39;,num2str(n3)]);
计算结果为
(1)
c=1时
Jacobi迭代与精确解得差为0、0021999
迭代次数为11796
G-S迭代与精确解得差为0、0017027
迭代次数为6227
SOR迭代与精确解得差为0、004511
迭代次数为15367
(2)
c=0、1时
Jacobi迭代与精确解得差为0、0094349
迭代次数为5353
G-S迭代与精确解得差为0、0093007
迭代次数为2797
SOR迭代与精确解得差为0、010279
迭代次数为7300
(3)
c=0、01时
Jacobi迭代与精确解得差为0、066098
迭代次数为532
G-S迭代与精确解得差为0、066089
迭代次数为318
SOR迭代与精确解得差为0、06615
迭代次数为834
(4)
c=0、0001时
Jacobi迭代与精确解得差为0、0049526
迭代次数为116
G-S迭代与精确解得差为0、0049507
迭代次数为108
SOR迭代与精确解得差为0、0049789
迭代次数为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&lt;10-7、 1=&quot;&quot; ax=&quot;b,由于g(x,y)与f(x,y)分别为exp(xy)与x+y,可得A得形式为&quot; ex4_2=&quot;&quot; c=&quot;10^(-7);&quot; n=&quot;1000000;&quot; for=&quot;&quot; m=&quot;1:3&quot; h=&quot;1/n(m);&quot; a=&quot;zeros((n(m)-1)^2);&quot; b=&quot;zeros((n(m)-1)^2,1);&quot; i=&quot;&quot; if=&quot;&quot;&gt;1
&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;A(i-1,i)=-1; &nbsp;A(i,i-1)=-1;
&nbsp; &nbsp; &nbsp; &nbsp;end
&nbsp; &nbsp; &nbsp; &nbsp;if i&gt;n(m)-1
&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;A(i,i-n(m)+1)=-1; &nbsp;A(i-n(m)+1,i)=-1;
&nbsp; &nbsp; &nbsp; &nbsp;end
&nbsp; &nbsp; &nbsp; &nbsp;ii=ceil(i/(n(m)-1));
&nbsp; &nbsp; &nbsp; &nbsp;if mod(i,n(m)-1)~=0
&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;jj=mod(i,n(m)-1);
&nbsp; &nbsp; &nbsp; &nbsp;else
&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;jj=n(m)-1;
&nbsp; &nbsp; &nbsp; &nbsp;end
&nbsp; &nbsp; &nbsp; &nbsp;A(i,i)=4+exp(ii*jj*h^2);
&nbsp; &nbsp; &nbsp; &nbsp;b(i)=h^3*(ii+jj);
&nbsp; &nbsp; &nbsp; &nbsp;if ii==1||ii==n(m)-1
&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;b(i)=b(i)+1;
&nbsp; &nbsp; &nbsp; &nbsp;end
&nbsp; &nbsp; &nbsp; &nbsp;if jj==1||jj==n(m)-1
&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;b(i)=b(i)+1;
&nbsp; &nbsp; &nbsp; &nbsp;end
&nbsp; &nbsp;end
&nbsp; &nbsp;disp([&#39;n=&#39;,num2str(n(m))])
&nbsp; &nbsp;tic
&nbsp; &nbsp;[y,k]=GaussSeidel(A,b,c,N);
&nbsp; &nbsp;toc
&nbsp; &nbsp;disp([&#39;迭代次数为&#39;,num2str(k)]);
end
结果为
n
20
40
80
CPU time
0、080063 seconds
1、941715 seconds
112、580604 seconds
迭代次数
24
26
27

<!--10-7、--></w<2></c,算法停止,输出解与迭代次数k,否则,转step4>