matlab上机作业报告(计算初等反射阵,用Householder变换法对矩阵A(4)

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

百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说综合文库matlab上机作业报告(计算初等反射阵,用Householder变换法对矩阵A(4)在线全文阅读。