1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1.0000 1.0011 1.0022 1.0031 1.0039 1.0044 1.0047 1.0045
1.0000 1.0022 1.0042 1.0061 1.0076 1.0087 1.0091 1.0088
1.0000 1.0031 1.0061 1.0088 1.0110 1.0126 1.0133 1.0128
1.0000 1.0039 1.0076 1.0110 1.0138 1.0159 1.0168 1.0162
1.0000 1.0044 1.0087 1.0126 1.0159 1.0189
1.0000 1.0047 1.0091 1.0133 1.0168 1.0203
1.0000 1.0045 1.0088 1.0128 1.0162 1.0201
1.0000 1.0037 1.0073 1.0107 1.0136 1.0175
1.0000 1.0023 1.0045 1.0066 1.0084 1.0113
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Columns 9 through 11
1.0000 1.0000 1.0000 1.0037 1.0023 1.0000 1.0073 1.0045 1.0000 1.0107 1.0066 1.0000 1.0136 1.0084 1.0000 1.0160 1.0100 1.0000 1.0174 1.0110 1.0000 1.0175 1.0113 1.0000 1.0155 1.0103 1.0000 1.0103 1.0072 1.0000 1.0000 1.0000 1.0000
>> contourf (u, 'DisplayName', 'u'); figure(gcf)
1.0183 1.0194 1.0194 1.0208 1.0189 1.0203 1.0160 1.0174 1.0100 1.0110 1.0000 1.0000
11109876543211234567891011
图二 块Gauss-Sediel迭代法
三 （预条件）共轭斜量法求解线性方程组Au=f。 1.基本原理
（1）预条件共轭斜量法原理 ?1?TT?1Au?b?CACCu?Cb?Au?b
?1?T?1T其中A?CAC,b?Cb,u?Cu，
令M=CC称为预优矩阵，当M接近A时，A接近单位阵，cond(A)2接近1，对Au?f用共轭斜量法求解，可达到加速的目的。T（2）预优矩阵的选取
?A212??IA???????IA33????I????I??AN?1,N?1??
M取为块Jacobi迭代的块对角阵，?A22?M???????????An?1,n?1??A33?
2. 算法
(1)取初值u(0)?R(n),r(0)(0)?b?Au(0)(0),(0)(2)解方程组Mz(3)k?0,1,...,(0)?r,求出z，令p?z(0)?k?ur(k?1)(k?1)(z(k),r(k)))(k)(k)(k?1)(Ap?u?r(k)(k),p(k)??kp(k)??kAp(k?1)(k?1)(k)解方程组Mz?r)?k?(4)直到r(k?1)(z(k?1)(k),r,r(z)(k?1)p。(k?1)?z(k?1)??kp(k)??,输出u3.程序
%用J-PCG求解正方形域上的Poisson方程边值问题
%输入n,对x、y轴进行n等分；先确定场u的边界及场源b； %用k和u分别返回迭代次数和精确解 function [k,u]=J_CG(n)
h=1/n; %h步长
u=zeros(n+1,n+1); %对u赋初值
u(1,1:n+1)=1;u(n+1,1:n+1)=1;u(1:n+1,1)=1;u(1:n+1,n+1)=1; b=zeros(n+1,n+1); %计算场源b for i=2:n+1 for j=2:n+1
b(i,j)=(i-1)*(j-1)*h^2; end end b=h^2*b; for i=2:n
b(2,i)=b(2,i)+u(1,i);b(i,n)=b(i,n)+u(i,n+1);
b(n,i)= b(n,i)+u(n+1,i);b(i,2)= b(i,2)+u(i,1);
end
z=zeros(n-1,n-1);r=zeros(n-1,n+1);p=zeros(n-1,n-1);bb=zeros(n-1,n-1); A=zeros(n-1,n-1); %定义矩阵的子块A for i=1:n-1
if i>1, A(i,i-1)=-1; end if i end for j=2:n if j==2, bb(:,1)=A*u(2:n,2)-u(2:n,3); end ?=A*u if j==n, bb(:,n-1)=A*u(2:n,n)-u(2:n,n-1); end if j>2&j r=b(2:n,2:n)-bb; %计算r0 for i=1:n-1 %计算z0 z(:,i)=pinv(A)*r(:,i); end p=z; %计算p0 tic; e=0.000000001; %e是精度 for k=1:1000 error=0;a=0;aa=0;cc=0;ub=u; for i=1:n-1 aa=aa+dot(z(:,i),r(:,i)); if i==1, cc=cc+dot((A*p(:,i)-p(:,i+1)),p(:,i));end if i>1&i a=aa/cc; %a是(z(k),r(k))/(Ap(k),p(k))，确定搜索步长 u(2:n,2:n)= u(2:n,2:n)+a*p; %对u进行更新计算 rr=r;zz=z; %rr和zz存储第k-1次的r和z for i=1:n-1 if i==1, r(:,i)= r(:,i)-a*(A*p(:,i)-p(:,i+1));end if i>1&i z(:,1:n-1)=pinv(A)*r(:,1:n-1); %Mz kk=0;mm=0; for i=1:n-1 kk=kk+dot(z(:,i),r(:,i)); mm=mm+dot(zz(:,i),rr(:,i)); end beita=kk/mm; ?ita是(z(k?1)(k?1)?r(k?1) ,r(k?1))/(z(k),r(k)) p=z+beita*p; %对p进行更新，确定搜索方向 error=sqrt(norm(u-ub)); %error是统计误差 if error<=e, break; end %判断误差是否达到精度 end t=toc %统计计算时间 4. 计算结果： >> format short >> n=10; >> [k,u]=J_CG(n) t = 0.1100 k = 37 u = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0011 1.0022 1.0031 1.0039 1.0045 1.0037 1.0023 1.0000 1.0000 1.0022 1.0042 1.0061 1.0076 1.0088 1.0073 1.0045 1.0000 1.0000 1.0031 1.0061 1.0088 1.0110 1.0128 1.0107 1.0066 1.0000 1.0000 1.0039 1.0076 1.0110 1.0138 1.0162 1.0136 1.0084 1.0000 1.0000 1.0044 1.0087 1.0126 1.0159 1.0189 1.0160 1.0100 1.0000 1.0000 1.0047 1.0091 1.0133 1.0168 1.0203 1.0174 1.0110 1.0000 1.0000 1.0045 1.0088 1.0128 1.0162 1.0201 1.0175 1.0113 1.0000 1.0000 1.0037 1.0073 1.0107 1.0136 1.0175 1.0155 1.0103 1.0000 1.0000 1.0023 1.0045 1.0066 1.0084 1.0113 1.0103 1.0072 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 >> contourf (u, 'DisplayName', 'u'); figure(gcf) 1.0000 1.0000 1.0044 1.0047 1.0087 1.0091 1.0126 1.0133 1.0159 1.0168 1.0183 1.0194 1.0194 1.0208 1.0189 1.0203 1.0160 1.0174 1.0100 1.0110 1.0000 1.0000