数值线性代数第二版徐树方高立张平文上机习题第四章实验报告

《数值线性代数第二版徐树方高立张平文上机习题第四章实验报告》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。

1、第四章上机习题1考虑两点边值问题d2ydy八t£H—61,0<6Z<1dx~dxy(O)=O,y(l)=l・容易知道它的精确解为—a-一y=(1-e£)+oxl-/7为了把微分方程离散化,把［0/L］区间n等分,令h二1/n,得到差分方程简化为£込戛如+沁jh2（£+〃）X+i一（2£+力）必+6L=。力2,从而离散化后得到的线性方程组的系数矩阵为-(2r+A)E+hE一(2e+/z)£+hE一(2e+/i)*•.*•.••・€+h£—(2£+/z)对£=1,Q=1/2/=100,分别用Jacobi迭代法,G-S迭代法和SOR迭代

2、法求线性方程组的解，要求有4位有效数字，然后比较与精确解得误差。对£=0.1,£=0.01,£=0.0001,考虑同样的问题。解（1）给岀算法：为解Ax=b,令A=D-L-U,其中A=[aij],D=diag（a{[,a22,...,ann）°-a]20一…~a23…一仏~Cl2,iu=••••••~an-,n0•••0利用Jacobi迭代法，G・S迭代法，SOR迭代法解线性方程组，均可以下步骤求解:stepl给定初始向量x0=(0,0,...,0),最大迭代次数N,精度要求c,令kJstep2令x=B*xO+gstep3若

3、

4、x

5、-xO

6、

7、2=N,算法停止，迭代失败，否则，令xO=x,转step2在Jacobi迭代法中，B=D1*(L+U)/g=D1*b在G-S迭代法中，B=D1*(L+U),g=D1*b在SOR迭代法中,B=(D-w*L)1*[(l-w)*D+w*U],g=w*(D-w*L)'1*b另外，在SOR迭代法屮，上面算法stepl中要给定松弛因子w,其中0

8、y(0)、y(l)的存在，使方程右端应减去y(0)、y(l)项，即b(l)=a*h2-c*y(0),b(n-l)=a*h2-(c+h)*y(n),b(i)=a*h2,i=2,.../n-l程序为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二tr订(A)-A;B二D"(-1)*(L+U);g二丁(T)*b;xO=zeros(length(A),1);x二B*xO+g;k=l;whilenorm

9、(x-xO,2)>=0.00001x0=x;x=B*xO+g;k=k+l;ifk>二Nbreakendendend2G-S迭代法编成的函数[x,k]二GaussSeidel(A,b,c,N)function[x,k]二GaussSeidel(A,b,c,N)U=diag(diag(A))-triu(A);xO=zeros(length(A),1);B=tril(A厂(-1)*U;g=tril(A厂(T)*b;x二B*xO+g;k=l;whilenorm(x-xO,2)>=0.00001x0=x;x=B*xO+g;k二k+1;ifk>

10、=Nbreakendendend1SOR迭代法编成的函数[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厂(T)*((1-w)*D+w*U);g=w*(D-w*L厂(T)*b;x二B*xO+g;k=l;whilenorm(x~x0,2)>=0.00001x0=x;x=B*x0+g;k二k+1;ifk>=Nbreakendendend4问题1求解ex

11、4_lclear;clc;%c=l;%c=0.1%c=0.01;c=0.0001;a=l/2;n=100;h=l/n;w二1/2;N二1000000;A二-(2*c+h)*eye(n-1);fori=2:n-lwA(i-1,i)二c+h;A(i,i-l)=c;endb二[a*h"2*ones(n-2,1);a*h"2-(c+h)];fori=l:n-lx(i)=i*h;y(i)=((l-a)/(l-exp(-1/c)))*(1-exp(-x(i)/c))+a*x(i);end[yl,nl]=Jacobi(A,b,c,N);[y2,n

12、2]=GaussSeidel(A,b,c,N);[y3,n3]=S0R(A,b,w,c,N);disp(['c二'，num2str(c),'时']);disp(['Jacobi迭代与精确解的差为'，num2str(norm(y,-yl