数值线性代数第二版徐树方高立张平文上机习题第四章实验报告.docx
- 文档编号:24490999
- 上传时间:2023-05-28
- 格式:DOCX
- 页数:10
- 大小:72.58KB
数值线性代数第二版徐树方高立张平文上机习题第四章实验报告.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||2 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)给定方程 在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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 线性代数 第二 版徐树方高立张平文 上机 习题 第四 实验 报告