1、数值线性代数第二版徐树方高立张平文上机习题第四章实验报告第四章上机习题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=1step2令x=B*x0+gstep3若|
2、x-x0|2=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，其中0w= x0=x; x=B*x0+g; k=k+1; if k=N break endendend2 G-S迭代法编成的函数 x,k=GaussSeidel(A,b,c,N)function x,k=GaussSeidel(A,b,c,N)U=diag
3、(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)= x0=x; x=B*x0+g; k=k+1; if k=N break endendend 3 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
4、*(D-w*L)(-1)*b;x=B*x0+g;k=1;while norm(x-x0,2)= x0=x; x=B*x0+g; k=k+1; if k=N break endendend4 问题1求解 ex4_1clear;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);for i=2:n-1w A(i-1,i)=c+h; A(i,i-1)=c;endb=a*h2*ones(n-2,1);a*h2-(c+h);for i=1:n-1 x(i)=i*h; y(i)=(1-a)/(1-exp(-1/c
5、)*(1-exp(-x(i)/c)+a*x(i);endy1,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(迭代次
6、数为,num2str(n3);计算结果为（1）c=1时Jacobi迭代与精确解的差为迭代次数为11796G-S迭代与精确解的差为迭代次数为6227SOR迭代与精确解的差为迭代次数为15367（2）c=时Jacobi迭代与精确解的差为迭代次数为5353G-S迭代与精确解的差为迭代次数为2797SOR迭代与精确解的差为迭代次数为7300（3）c=时Jacobi迭代与精确解的差为迭代次数为532G-S迭代与精确解的差为迭代次数为318SOR迭代与精确解的差为迭代次数为834（4）c=时Jacobi迭代与精确解的差为迭代次数为116G-S迭代与精确解的差为迭代次数为108SOR迭代与精确解的差为迭代次
7、数为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|21 A(i-1,i)=-1; A(i,i-1)=-1; end if in(m)-1 A(i,i-n(m)+1)=-1;
8、 A(i-n(m)+1,i)=-1; end ii=ceil(i/(n(m)-1); if mod(i,n(m)-1)=0 jj=mod(i,n(m)-1); else jj=n(m)-1; end A(i,i)=4+exp(ii*jj*h2); b(i)=h3*(ii+jj); if ii=1|ii=n(m)-1 b(i)=b(i)+1; end if jj=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结果为n204080CPU time seconds seconds seconds迭代次数242627
