第三章编程作业.docx
- 文档编号:3508615
- 上传时间:2022-11-23
- 格式:DOCX
- 页数:13
- 大小:160.61KB
第三章编程作业.docx
《第三章编程作业.docx》由会员分享,可在线阅读,更多相关《第三章编程作业.docx(13页珍藏版)》请在冰豆网上搜索。
第三章编程作业
数值实验题三
3.用有限差分方法(五点差分格式)求解正方形域上的Poisson方程边值问题
用MATLAB语言编写求解线性方程组
的算法程序,采用下列方法计算,比较计算结果和算法性能,对计算结果给出结论。
(1)用Jacobi迭代法求解线性方程组
。
(2)用块Jacobi迭代法求解线性方程组
。
(3)用(预条件)共轭斜向量法求解线性方程组
。
解:
由差分格式可得:
写成矩阵形式:
Au=f
其中:
其中:
(1)Jacobi迭代法:
function[U,k,er,t]=Jcb(n)
%Jacobi迭代法
%U表示方程组的解;h表示步长;A表示迭代矩阵;k表示迭代次数;n表示非边界点数
%f表示线性方程组A*U=f的右端矩阵f;e表示允许误差界;er表示迭代误差
%t表示计算时间
h=1/(n+1);
f(2:
n+1,2:
n+1)=2*h^2;%初始化f
A=zeros(n+2);%初始化A为n+2阶零矩阵
U=zeros(n+2);%初始化U为n+2阶零矩阵
e=0.000000001;%设置误差界
tic;
fork=1:
10000%迭代求解
er=0;
fori=2:
n+1
forj=2:
n+1
A(i,j)=(U(i-1,j)+U(i+1,j)+U(i,j+1)+U(i,j-1)+f(i,j))/4;
er=er+abs(A(i,j)-U(i,j));%估计当前误差
end
end
U=A;%将矩阵A的值赋予U
er=er/n^2;
ifer end end t=toc;%获得计算时间 >>[U,k,er,t]=Jcb(9) k= 304 er= 9.7771e-010 t= 0.0077 U如下: 对U作图: (2)块Jacobi方法: function[U,k,er,t]=KJcb(n) %块Jacobi迭代法 %u表示方程组的解h表示步长;k表示迭代次数;n表示非边界点数; %f表示线性方程组A*u=f的右端矩阵f;q表示n+2维向量;a表示方程组系数矩阵的下对角线元素 %b表示方程组系数矩阵的主对角线元素;c表示方程组系数矩阵的上对角线元素;d表示追赶法所求方程的右端向量 %e表示允许误差界;er表示迭代误差;l表示系数矩阵A所分解成的下三角阵L中的下对角线元素l(i);z表示系数矩阵A所分解成的上三角阵U中的主对角线元素z(i) h=1/(n+1); f(2: n+1,2: n+1)=2*h^2;%初始化f A=zeros(n+2);%初始化A为n+2阶零矩阵 U=zeros(n+2);%初始化U为n+2阶零矩阵 e=0.000000001;%设置误差界 a=-1*ones(1,n); b=4*ones(1,n); c=-1*ones(1,n); tic; fork=1: 1000%迭代求解 er=0; forj=2: n+1 q=zeros(1,n); d=f(2: n+1,j)+U(2: n+1,j-1)+U(2: n+1,j+1);%块Jacobi迭代 d=d'; A(2: n+1,j)=zg(a,b,c,d); er=er+norm(A(: j)-U(: j),1);%估计误差 end U=A; er=er/n^2; ifer end t=toc; 追赶法: functionx=zg(a,b,c,d) n=length(b); u (1)=b (1);y (1)=d (1); fori=2: n l(i)=a(i)/u(i-1); u(i)=b(i)-l(i)*c(i-1); y(i)=d(i)-l(i)*y(i-1);%追赶法求解之追过程求解Ly=d end x(n)=y(n)/u(n);%追赶法求解之赶过程求解Uz=y form=n-1: -1: 1 ifu(m)==0,D=0,break;end x(m)=(y(m)-c(m)*x(m+1))/u(m); end 运行结果: [U,k,er,t]=KJcb(9) 其中U如下: k= 163 er= 9.5903e-010 t= 0.1329 (3)共轭斜量法 functionA=xishu(N)%存储离散化后非边界点构成的系数矩阵 A=zeros(N^2); fori=1: N^2 A(i,i)=4; ifi+N A(i,i+N)=-1; A(i+N,i)=A(i,i+N); end ifmod(i,N)~=0 A(i,i+1)=-1; A(i+1,i)=A(i,i+1); end end functionx=cg(a,b,x) %共轭向量法求解线性方程组 %初始输入变量a: 系数矩阵b: 方程右端向量x输入的初始解 n=length(b); r=b-a*x; p=r; q0=r'*r; whileq0>1e-10 w=a*p; t=q0/(p'*w); x=x+t*p; r=r-t*w; q=r'*r; s=q/q0; p=r+s*p; q0=q; end 运行: >>a=xishu(9); fori=1: 9^2 b(i,1)=0.02; x(i,1)=0; end >>x=cg3(a,b,x) reshape(x,9,9) 结果: 下面的数据均为非边界点的值。 4.用有限差分方法(五点差分格式)求解正方形域上的Poisson方程边界问题 用Matlab语言编写求解线性方程组 的算法程序,采用下列方法计算,并比较方法的计算速度。 (1)用SOR迭代法求解线性方程组 ,用试算法确定最佳松弛因子。 (2)用块SOR迭代法求解线性方程组 ,用试算法确定最佳松弛因子。 (3)用(预条件)共轭斜向量法求解差分方程组 。 解: (1)SOR function[U,k,er,t]=SOR(n,w) %SOR迭代法 %U表示方程组的解;h表示步长;k表示迭代次数;n表示非边界点数 %f表示线性方程组A*U=f的右端矩阵f;e表示允许误差界;er表示迭代误差 %t表示计算时间 h=1/(n+1); f(2: n+1,2: n+1)=3*h^2;%初始化f U=zeros(n+2);%初始化U矩阵 forl=1: n+2 U(1,l)=(l-1)*h*(1-(l-1)*h); U(n+2,l)=(l-1)*h*(1-(l-1)*h); end e=0.000000001;%设置误差界 tic; fork=1: 10000%迭代求解 er=0; fori=2: n+1 forj=2: n+1 Ub=U(i,j); U(i,j)=w*((U(i-1,j)+U(i+1,j)+U(i,j+1)+U(i,j-1)+f(i,j))/4)+(1-w)*U(i,j); er=er+abs(Ub-U(i,j));%估计当前误差 end end er=er/n^2; ifer end end t=toc;%获得计算时间 最佳松弛因子: function[goodw,goodk]=zuijiaw(n) %计算最佳w %精度为0.01 [U,k,er,t]=SOR(n,0.1); goodk=k;goodw=0.01; forw=0.02: 0.01: 2 [U,k,er,t]=SOR(n,w); ifk<=goodk goodk=k; goodw=w; end end (2)块SOR function[U,k,er,t]=KSOR(n,w) %块SOR迭代法 %u表示方程组的解h表示步长;k表示迭代次数;n表示非边界点数; %f表示线性方程组A*u=f的右端矩阵f;q表示n+2维向量;a表示方程组系数矩阵的下对角线元素 %b表示方程组系数矩阵的主对角线元素;c表示方程组系数矩阵的上对角线元素;d表示追赶法所求方程的右端向量 %e表示允许误差界;er表示迭代误差;l表示系数矩阵A所分解成的下三角阵L中的下对角线元素l(i);z %表示系数矩阵A所分解成的上三角阵U中的主对角线元素z(i) h=1/(n+1); f(2: n+1,2: n+1)=3*h^2;%初始化f U=zeros(n+2);%初始化U矩阵 forl=1: n+2 U(1,l)=(l-1)*h*(1-(l-1)*h); U(n+2,l)=(l-1)*h*(1-(l-1)*h); end e=0.000000001;%设置误差界 a=-1*ones(1,n); b=4*ones(1,n); c=-1*ones(1,n); tic; fork=1: 1000%迭代求解 er=0; forj=2: n+1 Ub=U(: j); q=zeros(1,n); d=f(2: n+1,j)+U(2: n+1,j-1)+U(2: n+1,j+1);%块SOR迭代 d=d'; x=zg(a,b,c,d); fori=1: n U(i+1,j)=w*x(i)+(1-w)*U(i+1,j); end er=er+norm(Ub-U(: j),1);%估计误差 end er=er/n^2; ifer end t=toc; 最佳松弛因子: function[goodw,goodk]=zuijiaw(n) %计算最佳w %精度为0.01 [U,k,er,t]=KSOR(n,0.01); goodk=k;goodw=0.01; forw=0.02: 0.01: 2 [U,k,er,t]=KSOR(n,w); ifk<=goodk goodk=k; goodw=w; end end (3)共轭斜量法 functionA=xishu(N)%存储离散化后非边界点构成的系数矩阵 A=zeros(N^2); fori=1: N^2 A(i,i)=4; ifi+N A(i,i+N)=-1; A(i+N,i)=A(i,i+N); end ifmod(i,N)~=0 A(i,i+1)=-1; A(i+1,i)=A(i,i+1); end end functionU=cg(n) %共轭向量法求解线性方程组 h=1/(n+1); U=zeros(n+2);%初始化U矩阵 forl=1: n+2 U(1,l)=(l-1)*h*(1-(l-1)*h); U(n+2,l)=(l-1)*h*(1-(l-1)*h); end fori=1: n^2 b(i,1)=0.03; x(i,1)=0; end a=xishu(n); n=length(b); r=b-a*x; p=r; q0=r'*r; whileq0>1e-10 w=a*p; t=q0/(p'*w); x=x+t*p; r=r-t*w; q=r'*r; s=q/q0; p=r+s*p; q0=q; end N=n^0.5 xx=reshape(x,N,N); U(2: N+1,2: N+1)=xx;
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第三 编程 作业