数值分析上机作业

CHINAUNIVERSITYOFPETROLEUM数值分析上机作业姓名：黄亚航学号：2013213004班级：5班院系：化学工程学院任课老师：张明完成日期:2013年11月11日1数值分析上机作业一一、数值求解正方形域上的Poisson方程边值问题2222(,)1,0,1(0,)(1,)(1),01(,0)(,1)(1),01uufxyxyxyuyuyyyyuxuxxxx二、由椭圆型第一边值问题的五点差分格式得到线性方程组为2,1,1,,1,10,1,,0,14,1,(1),(1),01,0,1(1),(1),01,0,1ijijijijijijjNjiiNuuuuuhfijNuyyuyyyijNuxxuxxxijN写成矩阵形式Au=f。其中：1122NNAIIAAIIA1122NNvbvbufvb4114114iiA11,12,1,121,22,2,21,2,,2211,12,1,121,22,2,221,2,,,(,,...,),(,,...,),......,(,,...,)(,,...,),(,,...,),......,(,,...,)1,999,0.10.0111TTNNTNNNNNTTNNTNNNNNijvuuuvuuuvuuubhfffbhfffbhfffhNhNf,,1,2,...,ijN三、编写求解线性方程组Au=f的算法程序，用下列方法编程计算，并比较计算速度。1.用Jacobi迭代法求解线性方程组Au=f。2.用块Jacobi迭代法求解线性方程组Au=f。3.用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。4.用块SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。5.用（预条件）共轭斜量法求解差分方程组Au=b。四、上机报告要求1.简述方法的基本原理。2.程序中要加注释。3.对程序中的主要变量给出说明。4.附原程序及计算结果。5.迭代法效率分析：比较迭代法的效率，记录所需的迭代次数，所花费的计算时间。2（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)=h^2;%初始化fA=zeros(n+2);%初始化A为n+2阶零矩阵U=zeros(n+2);%初始化U为n+2阶零矩阵forx=1:n+2%给U赋上边界条件U(1,x)=(x-1)*h*(1-(x-1)*h);U(n+2,x)=(x-1)*h*(1-(x-1)*h);U(x,1)=(x-1)*h*(1-(x-1)*h);U(x,n+2)=(x-1)*h*(1-(x-1)*h);ende=0.000000001;%设置误差界tic;fork=1:1000%迭代求解er=0;fori=2:n+1forj=2:n+1A(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));%估计当前误差endendU(2:n+1,2:n+1)=A(2:n+1,2:n+1);%将矩阵A的值赋予Uer=er/n^2;ifere,break;end%判断是否达到计算精度，如果达到则退出循环endt=toc;%获得计算时间运行结果：[U,k,er,t]=Jcb(9)k=322er=9.9052e-010t=0.0032U如下：3对U作图：（2）块Jacobi迭代法：function[U,k,er,t]=KJcb(n)%块Jacobi迭代法%U表示方程组的解；h表示步长；A表示迭代矩阵；k表示迭代次数；n表示非边界点数；%f表示线性方程组A*u=f的右端矩阵f；r表示余矩阵；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)=h^2;%初始化fU=zeros(n+2);%初始化U为n+2阶零矩阵forx=1:n+2%给U赋上边界条件U(1,x)=(x-1)*h*(1-(x-1)*h);U(n+2,x)=(x-1)*h*(1-(x-1)*h);U(x,1)=(x-1)*h*(1-(x-1)*h);U(x,n+2)=(x-1)*h*(1-(x-1)*h);endA=U;%初始化Ar=zeros(n);%初始化余矩阵r(1,:)=U(1,2:n+1);r(n,:)=U(n+2,2:n+1);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+1d=f(2:n+1,j)+U(2:n+1,j-1)+U(2:n+1,j+1)+r(:,j-1);%块Jacobi迭代d=d';A(2:n+1,j)=zg(a,b,c,d);er=er+norm(A(:,j)-U(:,j),1);%估计误差4endU=A;er=er/n^2;ifere,break;end%判断是否达到计算精度，如果达到则退出循环endt=toc;%获得计算时间追赶法程序：functionx=zg(a,b,c,d)n=length(b);u(1)=b(1);%LU分解。fork=2:nifu(k-1)==0,D=0,return;endl(k)=a(k)/u(k-1);u(k)=b(k)-l(k)*c(k-1);endy(1)=d(1);%追赶法求解之追过程，求解Ly=d。fork=2:ny(k)=d(k)-l(k)*y(k-1);endifu(n)==0,D=0,return;endx(n)=y(n)/u(n);%追赶法求解之赶过程，求解Ux=yfork=n-1:-1:1x(k)=(y(k)-c(k)*x(k+1))/u(k);end运行结果：[U,k,er,t]=KJcb(9)k=172er=9.8209e-010t=0.0787U如下：5对U作图：（3）SOR迭代法：function[U,k,er,t]=SOR(n,w)%SOR迭代法%U表示方程组的解；h表示步长；A表示迭代矩阵；k表示迭代次数；n表示非边界点数%f表示线性方程组A*U=f的右端矩阵f；e表示允许误差界；er表示迭代误差%t表示计算时间h=1/(n+1);f(2:n+1,2:n+1)=h^2;%初始化fU=zeros(n+2);%初始化U为n+2阶零矩阵forx=1:n+2%给U赋上边界条件U(1,x)=(x-1)*h*(1-(x-1)*h);U(n+2,x)=(x-1)*h*(1-(x-1)*h);U(x,1)=(x-1)*h*(1-(x-1)*h);U(x,n+2)=(x-1)*h*(1-(x-1)*h);ende=0.000000001;%设置误差界tic;fork=1:1000%迭代求解er=0;fori=2:n+1forj=2:n+1Ub=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));%估计当前误差endender=er/n^2;ifere,break;end%判断是否达到计算精度，如果达到则退出循环endt=toc;%获得计算时间最佳松弛因子：function[goodw,goodk]=zuijiaw(n)%计算最佳w，精度为0.01[U,k,er,t]=SOR(n,0.1);goodk=k;goodw=0.01;6forw=0.02:0.01:2[U,k,er,t]=SOR(n,w);ifk=goodkgoodk=k;goodw=w;endend运行结果：[goodw,goodk]=zuijiaw(9)goodw=1.5400goodk=35[U,k,er,t]=SOR(9,1.54)k=35er=7.4839e-010t=9.3024e-004U如下：对U作图：7（4）块SOR迭代法function[U,k,er,t]=KSor(n,w)%块SOR迭代法%U表示方程组的解；h表示步长；k表示迭代次数；n表示非边界点数；r表示余矩阵%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)=h^2;%初始化fU=zeros(n+2);%初始化U为n+2阶零矩阵forx=1:n+2U(1,x)=(x-1)*h*(1-(x-1)*h);U(n+2,x)=(x-1)*h*(1-(x-1)*h);U(x,1)=(x-1)*h*(1-(x-1)*h);U(x,n+2)=(x-1)*h*(1-(x-1)*h);endr=zeros(n);%初始化余矩阵r(1,:)=U(1,2:n+1);r(n,:)=U(n+2,2:n+1);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+1Ub=U(:,j);d=f(2:n+1,j)+U(2:n+1,j-1)+U(2:n+1,j+1)+r(:,j-1);%块SOR迭代x=zg(a,b,c,d');fori=1:nU(i+1,j)=w*x(i)+(1-w)*U(i+1,j);ender=er+norm(Ub-U(:,j),1);%估计误差ender=er/n^2;ifere,break;end%判断是否达到计算精度，如果达到则退出循环endt=toc;%获得计算时间追赶法程序：functionx=zg(a,b,c,d)n=length(b);u(1)=b(1);%LU分解。fork=2:nifu(k-1)==0,D=0,return;endl(k)=a(k)/u(k-1);u(k)=b(k)-l(k)*c(k-1);8endy(1)=d(1);%追赶法求解之追过程，求解Ly=d。fork=2:ny(k)=d(k)-l(k)*y(k-1);endifu(n)==0,D=0,return;endx(n)=y(n)/u(n);%追赶法求解之赶过程，求解Ux=yfork=n-1:-1:1x(k)=(y(k)-c(k)*x(k+1))/u(k);end最佳松弛因子：function[goodw,goodk]=zuijiaw(n)%计算最佳w，精度为0.01[U,k,er,t]=KSor(n,0.1);goodk=k;goodw=0.01;forw=0.02:0.01:2[U,k,er,t]=SOR(n,w);ifk=goodkgoodk=k;goodw=w;endend运行结果：[goodw,goodk]=zuijiaw(9)goodw=1.4200goodk