计算实习任务1.doc

用有限差分法求解正方形域上的Poisson方程边值问题 解：用五点差分格式对椭圆型第一边值问题得到线性方程组为 写成矩阵形式Au=f，其中 用Jacobi迭代法求解线性方程组Au=f。 Jacobi迭代法 对 k=0,1,2… 即 其中 迭代矩阵 右端向量 分量形式 （i=1,2,…n） 程序： %Jacobi迭代法 function [U,k,er,t]=jacobi(n) %变量说明： % 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; %初始化f U=zeros(n+2); %初始化U为n+2阶零矩阵 e=0.000000001; %设置误差界 for j=0:n+1 %设置矩阵A，U边界点的值 U(1,j+1)=j*h*(1-j*h); U(n+2,j+1)=U(1,j+1); U(j+1,1)=j*h*(1-j*h); U(j+1,n+2)=U(j+1,1); end A=U tic; for k=1:10000 %迭代求解 er=0; max=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=abs(A(i,j)-U(i,j)); %估计当前误差 if ermax max=er; end end end U=A; %将矩阵A的值赋予U if maxe,break; %判断是否达到计算精度，如果达到则退出循环 end end t=toc; %获得计算时间 end 命令窗口下输入： [U,k,er,t]=jacobi(9) （2）用块Jacobi迭代法求解线性方程组Au=f。 先把矩阵中的矩阵Aii看成Jacobi迭代法中的点 对 k=0,1,2… 即 其中 迭代矩阵 右端向量 分量形式 （i=1,2,…n） 第二步由于是三对角阵所以对每一个 （i=1,2,…n）用追赶法。追赶法原理： 追赶法： ①计算分解因子阵 对k=2,3,…,n ②求解LY=d， ③求解Ux=Y, 程序： %块Jacobi迭代法 function [U,k,er,t]=bjacobi(n) %变量说明： % u --- 方程组的解 % h --- 步长 % k --- 迭代次数 % n --- 非边界点数 % p --- n+2维向量 % q --- n+2维向量 % f --- 线性方程组A*u=f的右端矩阵f % 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; %初始化f 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); for j=0:n+1 %设置矩阵A，U边界点的值 U(1,j+1)=j*h*(1-j*h); U