编程作业.docVIP

 数值实验题三 3.用有限差分方法（五点差分格式）求解正方形域上的Poisson方程边值问题 用MATLAB语言编写求解线性方程组的算法程序，采用下列方法计算，比较计算结果和算法性能，对计算结果给出结论。 （1）用Jacobi迭代法求解线性方程组。 （2）用块Jacobi迭代法求解线性方程组。 （3）用（预条件）共轭斜向量法求解线性方程组。 解：由差分格式可得: 写成矩阵形式：Au=f 其中： 其中： （1）用Jacobi迭代法求解线性方程组Au=f： ①、方法的基本原理 对A作自然分解： A=D-L-U=D-(L+U) 其中D是由A的对角线元素组成的矩阵，L是由A的对角线以下的元素组成的矩阵，U是由A的对角线以上的元素组成的矩阵。 于是将M=D，N=L+U代入得到 Dx=(L+U)x+b 任取x的初值进行迭代。 ②、原程序及计算结果 第一步：编写文件名为Jacobi.m的M文件 function [U,k,er,t]=Jacobi(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; %tic用来保存当前时间 for k=1:10000 %迭代求解 er=0; for i=2:n+1 for j=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; if ere,break; %判断是否达到计算精度，如果达到则退出循环 end end t=toc; %记录程序完成时间 数值解见下图： （2）用块Jacobi迭代法求解线性方程组Au=f： ①、方法的基本原理 对A作自然分解： A=D-L-U=D-(L+U) 其中D是由A的对角线元素组成的矩阵，L是由A的对角线以下的元素组成的矩阵，U是由A的对角线以上的元素组成的矩阵。 于是将M=D，N=L+U代入得到 Dx=(L+U)x+b 任取x的初值进行迭代。 ②、原程序及运行结果 第一步：编写文件名为KJacobi.m的M文件 function [U,k,er,t]=KJacobi(n) %块Jacobi迭代法 % u 表示方程组的解h 表示步长；k 表示迭代次数；n 表示非边界点数； % f 表示线性方程组A*u=f的右端矩阵f；q 表示n+2维向量；a 表示方程组系数矩阵的下对角线元素 % b 表示方程组系数矩阵的主对角线元素；c 表示方程组系数矩阵的上对角线元素；d 表示追赶法所求方程的右端向量 % e 表示允许误差界；er 表示迭代误差；l 表示系数矩阵A所分解成的下三角阵L中的下对角线元素 l；z 表示系数矩阵A所分解成的上三角阵U中的主对角线元素 z () 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;tic用来保存当前时间for k=1:1000 % 迭代求解 er=0; for j=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+n