第三章編程作業(yè)

1、數(shù)值實(shí)驗(yàn)題三3.用有限差分方法（五點(diǎn)差分格式）求解正方形域上的Poisson方程邊值問題用MATLAB語言編寫求解線性方程組的算法程序，采用下列方法計算，比較計算結(jié)果和算法性能，對計算結(jié)果給出結(jié)論。（1）用Jacobi迭代法求解線性方程組。（2）用塊Jacobi迭代法求解線性方程組。（3）用（預(yù)條件）共軛斜向量法求解線性方程組。解：由差分格式可得:寫成矩陣形式：Au=f其中：其中：（1）Jacobi迭代法：function U,k,er,t=Jcb(n)% Jacobi迭代法% U 表示方程組的解；h 表示步長； A 表示迭代矩陣；k 表示迭代次數(shù)；n 表示非邊界點(diǎn)數(shù)% f 表示線性方程組A*

2、U=f的右端矩陣f ；e 表示允許誤差界；er 表示迭代誤差% t 表示計算時間h=1/(n+1);f(2:n+1,2:n+1)= 2*h2; %初始化fA=zeros(n+2); %初始化A為n+2階零矩陣U=zeros(n+2); %初始化U為n+2階零矩陣e=0.; %設(shè)置誤差界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); %估計當(dāng)前誤差 end end U=A; %將矩

3、陣A的值賦予U er=er/n2; if er U,k,er,t=Jcb(9)k = 304er = 9.7771e-010t =0.0077U如下：對U作圖： （2）塊Jacobi方法：function U,k,er,t=KJcb(n)%塊Jacobi迭代法% u 表示方程組的解h 表示步長；k 表示迭代次數(shù)；n 表示非邊界點(diǎn)數(shù)； % f 表示線性方程組A*u=f的右端矩陣f；q 表示n+2維向量；a 表示方程組系數(shù)矩陣的下對角線元素 % b 表示方程組系數(shù)矩陣的主對角線元素；c 表示方程組系數(shù)矩陣的上對角線元素；d 表示追趕法所求方程的右端向量% e 表示允許誤差界；er 表示迭代誤差；l

4、 表示系數(shù)矩陣A所分解成的下三角陣L中的下對角線元素 l(i)；z 表示系數(shù)矩陣A所分解成的上三角陣U中的主對角線元素 z(i)h=1/(n+1);f(2:n+1,2:n+1)=2*h2; %初始化fA=zeros(n+2); %初始化A為n+2階零矩陣U=zeros(n+2); %初始化U為n+2階零矩陣e=0.; %設(shè)置誤差界 a=-1*ones(1,n);b=4*ones(1,n);c=-1*ones(1,n); 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+

5、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/n2; if ere ,break;end %判斷是否達(dá)到計算精度，如果達(dá)到則退出循環(huán)endt=toc;追趕法：function x=zg(a,b,c,d)n=length(b);u(1)=b(1);y(1)=d(1);for i=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

6、=d end x(n)=y(n)/u(n); % 追趕法求解之趕過程 求解Uz=y for m=n-1:-1:1 if u(m)=0 ,D=0,break; end x(m)=(y(m)-c(m)*x(m+1)/u(m); end運(yùn)行結(jié)果：U,k,er,t=KJcb(9)其中U如下：k = 163er = 9.5903e-010t = 0.1329（3）共軛斜量法function A=xishu(N) %存儲離散化后非邊界點(diǎn)構(gòu)成的系數(shù)矩陣A=zeros(N2);for i=1:N2 A(i,i)=4; if i+N1e-10 w=a*p; t=q0/(p*w); x=x+t*p; r=r-t*

7、w; q=r*r; s=q/q0; p=r+s*p; q0=q;end運(yùn)行： a=xishu(9);for i=1:92b(i,1)=0.02;x(i,1)=0;end x=cg3(a,b,x)reshape(x,9,9)結(jié)果：下面的數(shù)據(jù)均為非邊界點(diǎn)的值。4.用有限差分方法（五點(diǎn)差分格式）求解正方形域上的Poisson方程邊界問題用Matlab語言編寫求解線性方程組的算法程序，采用下列方法計算，并比較方法的計算速度。（1）用SOR迭代法求解線性方程組，用試算法確定最佳松弛因子。（2）用塊SOR迭代法求解線性方程組，用試算法確定最佳松弛因子。（3）用（預(yù)條件）共軛斜向量法求解差分方程組。解：（1

8、）SORfunction U,k,er,t=SOR(n,w)% SOR迭代法% U 表示方程組的解；h 表示步長； k 表示迭代次數(shù)；n 表示非邊界點(diǎn)數(shù)% f 表示線性方程組A*U=f的右端矩陣f ；e 表示允許誤差界；er 表示迭代誤差% t 表示計算時間h=1/(n+1);f(2:n+1,2:n+1)= 3*h2; %初始化fU=zeros(n+2); %初始化U矩陣for l=1:n+2 U(1,l)=(l-1)*h*(1-(l-1)*h); U(n+2,l)=(l-1)*h*(1-(l-1)*h);ende=0.; %設(shè)置誤差界 tic;for k=1:10000 %迭代求解 er=0

9、; for i=2:n+1 for j=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); %估計當(dāng)前誤差 end end er=er/n2; if ere,break; %判斷是否達(dá)到計算精度，如果達(dá)到則退出循環(huán) endendt=toc; %獲得計算時間最佳松弛因子：function goodw,goodk=zuijiaw(n)%計算最佳w%精度為0.01U,k,er,t=SOR(n,0.1);goodk=k;goodw=0.01

10、;for w=0.02:0.01:2 U,k,er,t=SOR(n,w); if k=goodk goodk=k; goodw=w; endend（2）塊SORfunction U,k,er,t=KSOR(n,w)%塊SOR迭代法% u 表示方程組的解h 表示步長；k 表示迭代次數(shù)；n 表示非邊界點(diǎn)數(shù)； % f 表示線性方程組A*u=f的右端矩陣f；q 表示n+2維向量；a 表示方程組系數(shù)矩陣的下對角線元素 % b 表示方程組系數(shù)矩陣的主對角線元素；c 表示方程組系數(shù)矩陣的上對角線元素；d 表示追趕法所求方程的右端向量% e 表示允許誤差界；er 表示迭代誤差；l 表示系數(shù)矩陣A所分解成的下三

11、角陣L中的下對角線元素 l(i)；z% 表示系數(shù)矩陣A所分解成的上三角陣U中的主對角線元素 z(i)h=1/(n+1);f(2:n+1,2:n+1)=3*h2; %初始化fU=zeros(n+2); %初始化U矩陣for l=1:n+2 U(1,l)=(l-1)*h*(1-(l-1)*h); U(n+2,l)=(l-1)*h*(1-(l-1)*h);ende=0.; %設(shè)置誤差界 a=-1*ones(1,n);b=4*ones(1,n);c=-1*ones(1,n); tic;for k=1:1000 % 迭代求解 er=0; for j=2:n+1 Ub=U(:,j); q=zeros(1,

12、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); for i=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/n2; if ere ,break;end %判斷是否達(dá)到計算精度，如果達(dá)到則退出循環(huán)endt=toc;最佳松弛因子：function goodw,goodk=zuijiaw(n)%計算最佳w%精度為0.01U,k,er,t=KSOR(n,0.01);goodk=k;goodw=0.01;for w=0.02:0.01:2 U,k,er,t=KSOR(n,w); if k=goodk goodk=k; goodw=w; endend（3）共軛斜量法