matlab上机作业报告(计算初等反射阵-用Householder变换法对矩阵A作正交分解-连续函数最佳平方逼近等).doc

一、给定向量x≠0，计算初等反射阵Hk。
1．程序功能：
给定向量x≠0，计算初等反射阵Hk。
2．基本原理：
若的分量不全为零，则由

确定的镜面反射阵H使得；当时，由

有
算法：
(1)输入x，若x为零向量，则报错
(2)将x规范化，

如果M＝0，则报错同时转出停机
否则
(3)计算，如果，则
(4)
(5)计算
(6)
(7)
(8)按要求输出,结束

3．变量说明：
x － 输入的n维向量；
n － n维向量x的维数；
M － M是向量x的无穷范数，即x中绝对值最大的一项的绝对值；
p － Householder初等变换阵的系数ρ；
u － Householder初等变换阵的向量U
s － 向量x的二范数；
x － 输入的n维向量；
n － n维向量x的维数；
p － Householder初等变换阵的系数ρ；
u － Householder初等变换阵的向量U
k － 数k，H*x=y，使得y的第k+1项到最后项全为零；
4．程序代码：
(1)
function [p,u]=holder2(x)
%HOLDER2 给定向量x≠0,计算Householder初等变换阵的p,u
%程序功能：函数holder2给定向量x≠0,计算Householder初等变换阵的p,u；
%输入：n维向量x；
%输出：[p,u]。p是Householder初等变换阵的系数ρ，
% u是Householder初等变换阵的向量U。
n=length(x); % 得到n维向量x的维数；
p=1;u=0; % 初始化p,u；
M=max(abs(x)); % 得到向量x的无穷范数，即x中绝对值最大的一项的绝对值；
if M==0 % 如果x=0，提示出错,程序终止;
disp('Error: M=0');
return;
else
x=x/M; % 规范化
end;
s=norm(x); % 求x的二范数
if x(1)<0 % 首项为负，s值要变号
s=-s;
end
u=x; % 除首项外，其余各项x,u相同
u(1)=s+x(1); % 计算u的首项
p=s*u(1); % 计算p
if n==1 u=0; end % 若x是1×1维向量，则u=0

(2)
function H=holderk(x,k)
%HOLDERK 给定向量x≠0,数k，计算初等反射阵Hk，使HkX=Y，其中Y的第k+1项到最后项全为零；
%程序功能：函数holderk给定向量x≠0，数k，计算初等反射阵Hk，使HkX=Y，%程序功能：函数holder2给定向量x≠0,计算Householder初等变换阵的p,u；
%输入：n维向量x,数k；
%输出：H。H是Householder初等变换阵，H*x=y，使得y的第k+1项到最后项全为零；
%引用函数：holder2；
n=length(x); % 得到n维向量x的维数；
if k>n %如果k值溢出，报错；
disp('Error: k>n');
end

H=eye(n); % 初始化H，并使H(1:k,1:k)=I；
[p,u]=holder2(x(k:n)); % 得到计算Householde初等变换阵的系数ρ、向量U；
H(k:n,k:n)=eye(n-k+1)-p\u*u'; % 计算H(k:n,k:n)=I-p\u*u'；
5．使用示例:
情形1:
X为零向量
>> x=[0,0,0,0]';
>> H=holderk(x,1)
Error: M=0

H =

1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
情形2：
K值溢出：
>> x=[1,2,3,4]';
>> H=holderk(x,5)
Error: k>n
情形3：
K值为1：
>> x=[2,3,4,5]';
>> H=holderk(x,1)

H =

-0.2722 -0.4082 -0.5443 -0.6804
-0.4082 0.8690 -0.1747 -0.2184
-0.5443 -0.1747 0.7671 -0.2911
-0.6804 -0.2184 -0.2911 0.6361
检验:
>> det(H)

ans =

-1.0000

>> H*x

ans =

-7.3485
0.0000
0.0000
0.0000
情形4：
(1)K值为3：
>> x=[4,3,2,1]';
>> H=holderk(x,3)

H =

1.0000 0 0 0
0 1.0000 0 0
0 0 -0.8944 -0.4472
0 0 -0.4472 0.8944
检验:

>> det(H)

ans =

-1

>> H*x

ans =

4.0000
3.0000
-2.2361
0
(2)K值为2:
>> x=[4,3,2,1]';
>> H=holderk(x,2)

H =

1.0000 0 0 0
0 -0.8018 -0.5345 -0.2673
0 -0.5345 0.8414 -0.0793
0 -0.2673 -0.0793 0.9604

>> det(H)

ans =

-1.0000

>> H*x

ans =

4.0000
-3.7417
0.0000
0

二、设A为n阶矩阵，编写用Householder变换法对矩阵A作正交分解的程序。
1．程序功能：
给定n阶矩阵A，通过本程序用Householder变换法对矩阵A作正交分解，得出A＝QR
2．基本原理：
任一实列满秩的m×n矩阵A，可以分解成两个矩阵的乘积，即A＝QR，其中Q是具有法正交列向量的m×n矩阵，R是非奇异的n阶上三角阵。
算法：
(1)输入n阶矩阵A
(2)对，求Househoulder初等反射阵的。
(3)计算上三角阵R，仍然存储在A

(4)计算正交阵Q

(5)按要求输出,结束

3．变量说明：
A － 输入的n阶矩阵，同时用于存储上三角阵R；
n － 矩（方）阵A的阶数；
Q － Q是具有法正交列向量的n阶矩阵；
p,u － 向量A（k:n,k），对应初等反射阵的ρ,u
k，jj，ii － 循环变量；
t1 － 计算上三角阵R的系数tj；
t2 － 计算正交矩阵Q的系数ti；
4．程序代码：
function [Q,A]=qrhh(A)
%QRHH 用Householder变换法对n阶矩阵A作正交分解A=QR；
%程序功能：函数qrhh用Householder变换法对矩阵A作正交分解A=QR；
%输入：n阶矩阵A；
%输出：[Q,A]。Q是具有法正交列向量的n阶矩阵，
% A（即R）是非奇异的n阶上三角阵，仍用输入的矩阵A存储。
%引用函数：
% holder2；示例 [p,u]=holder2(x);
[n,n]=size(A); %求矩（方）阵A的阶数；
Q=eye(n); %构造正交矩阵Q(1)=I；
for k=1:n-1
[p,u]=holder2(A(k:n,k)); %向量A（k:n,k），对应初等反射阵的ρ,u

for jj=k:n %计算上三角阵R（仍存贮于A）
t1=dot(u,A(k:n,jj))/p; %利用向量内积求和
A(k:n,jj)=A(k:n,jj)-t1*u;
end
for ii=1:n %计算正交矩阵Q
t2=dot(u,Q(ii,k:n))/p; %利用向量内积求和
Q(ii,k:n)=Q(ii,k:n)-t2*u';
end
end

5．使用示例:
(1)A为3阶矩阵：
>> A=[1 2 3; 2 3 0; 3 4 5]

A =

1 2 3
2 3 0
3 4 5

>> [q,r]=qrhh(A)

q =

-0.2673 0.8729 0.4082
-0.5345 0.2182 -0.8165

-0.8018 -0.4364 0.4082

r =

-3.7417 -5.3452 -4.8107
0 0.6547 0.4364
-0.0000 0.0000 3.2660
检验：
>> q*r

ans =

1.0000 2.0000 3.0000
2.0000 3.0000 0.0000
3.0000 4.0000 5.0000

(2)A为4阶矩阵：
>> A=[1 2 3 4; 2 3 0 1; 3 4 5 6;1 6 8 0]

A =

1 2 3 4
2 3 0 1
3 4 5 6
1 6 8 0

>> [q,r]=qrhh(A)

q =

-0.2582 0.0597 -0.2660 -0.9268
-0.5164 -0.1045 0.8434 -0.1049
-0.7746 -0.2688 -0.4662 0.3323
-0.2582 0.9556 -0.0222 0.1399

r =

-3.8730 -6.7132 -6.7132 -6.1968
0 4.4647 6.4805 -1.4783
0 -0.0000 -3.3070 -3.0178
0 0.0000 0 -1.8187
检验：
>> q*r

ans =

1.0000 2.0000 3.0000 4.0000
2.0000 3.0000 -0.0000 1.0000
3.0000 4.0000 5.0000 6.0000
1.0000 6.0000 8.0000 0

数值求解正方形域上的Poisson方程边值问题

用MATLAB语言编写求解此辺值问题的算法程序，采用下列三种方法，并比较三种方法的计算速度。1、用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子；2、用块 Gauss-Sediel迭代法求解线性方程组Au=f；3、（预条件）共轭斜量法。
用差分代替微分，对Poisson方程进行离散化，得到五点格式的线性方程组

写成矩阵形式Au=f。
其中

一 用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。

1. 基本原理：
Gauss-Seidel迭代法计算简单，但是在实际计算中，其迭代矩阵的谱半径常接近1，因此收敛很慢。为了克服这个缺点，引进一个加速因子（又称松弛因子）对Gauss-Seidel方法进行修正加速。
假设已经计算出第k步迭代的解（i=1，2，···，n），要求下一步迭代的解（i=1，2，···，n）。首先，用Gauss-Seidel迭代格式计算

然后引入松弛因子，用松弛因子对和作一个线性组合。
，i=1,2,…,n
将二者合并成为一个统一的计算公式：

2. 算法
(1)Gauss-Seidel迭代法引入松弛因子w：

五点格式即为：

(2)计算步骤：
第一步：给松弛因子赋初值w=1.1～1.8,给场值u和场源b赋初值
第二步：用不同的w进行迭代计算。置error=0;
计算
在计算机上采用动态计算形式

如果|du|>error则error=|du|,如果error用kk存放最小迭代次数，用ww和uu存放相应的w及u。

3． 程序
① [u,k]=SOR(u,b,w) %%%%%%%（被下面程序调用）
%输入场初值u0、场源b及松弛因子w，通过五点差分格式进行迭代运算，
%如果第k+1次的迭代结果与第k次的差小于精度，则可以近似认为第k+1次的迭代
%结果是精确解，然后返回迭代次数k和迭代解
function [u,k]=SOR(u,b,w) %输出迭代结果u，及迭代次数k
m=length(u); %m为u的维数
h=1/(m-1); %h为步长
N=10000;e=0.0000001; %e为精度
for k=1:N %k为记录迭代次数
error=0;
for j=2:m-1
for jj=2:m-1
sum=4*u(jj,j)-u(jj-1,j)-u(jj+1,j)-u(jj,j-1)-u(jj,j+1);
du=w*(h^2*b(jj,j)-sum)/4; %计算u的修正量
u(jj,j)= u(jj,j)+du; %修正u
if error
end

② [kk,ww,uu]=SOR_5dianchafen(n)
%用超松弛迭代法求解正方形域上的Poisson方程边值问题
%用5点差分格式求取泊松问题
%输入n,对x、y轴进行n等分；先确定场u的边界及场源b,在调用[u,k]=SOR(u,b,w)；
%用不同w计算的迭代次数不同，用kk存放最小的迭代次数，
%用ww和uu分别存放最佳松弛因子w和精确解
function [kk,ww,uu]=SOR_5dianchafen(n)
w=[1.1:0.1:1.8];m=length(w); %w为松弛因子
kk=1000; ww=0; %kk是最少迭代次数，ww是最松弛因子
h=1/n; %h步长
b=zeros(n+1,n+1); %计算场源b
tic;
for i=2:n+1
for j=2:n+1
b(i,j)=(i-1)*(j-1)*h^2;
end
end
uu=zeros(n+1,n+1); 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;
mu=u; %初值mu以便不同的w计算
for i=1:m %用不同的w计算迭代
[u,k]=SOR(mu,b,w(i)); %调用[u,k]=SOR(u,b,w)，返回迭代次数及精确解
if kk>k, kk=k;ww=w(i);uu=u;end %把最少迭代次数付给kk,及其w,u赋给ww,uu
end
t=toc %统计程序运算时间

4.计算结果：
>> format short
>> n=10;
>> [kk,ww,uu]=SOR_5dianchafen(n)
t =
0.0310
kk =
48
ww =
1.6000
uu =
Columns 1 through 8
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.0183 1.0194 1.0189
1.0000 1.0047 1.0091 1.0133 1.0168 1.0194 1.0208 1.0203
1.0000 1.0045 1.0088 1.0128 1.0162 1.0189 1.0203 1.0201
1.0000 1.0037 1.0073 1.0107 1.0136 1.0160 1.0174 1.0175
1.0000 1.0023 1.0045 1.0066 1.0084 1.0100 1.0110 1.0113
1.0000 1.0000 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 (uu, 'DisplayName', 'uu'); figure(gcf)

图一 超松弛

二 用块Jacobi迭代法求解线性方程组Au=f。
1. 基本原理：
对A做自然分解A=D-L-U=D-(L+U)
其中D是有A的对角线元素组成的矩阵，L是由A的对角线以下元素组成的矩阵，U是由A得对角线以上元素组成的矩阵。
于是将M=D,N=L+U,代入得到
Dx=(L+U)x+b
任取x的初值进行迭代
2. 算法：
（1）Gauss-Sediel迭代法原理

五点差分格式：

因为A可以写成块状，即：
如果把每一条线上的节点看作一个组，可以把Au=f表示成块状求解：

（2）计算步骤：
第一步：给场值u和场源b赋初值，及定义
第二步：用公式，进行迭代计算
第三步：把第k次的u赋给ub，即ub=u；然后把第k+1次的u和ub进行比较，看是否达到精度，如果达到精度，则输出迭代次数k和精确解u。
3. 程序

[k,u]=kuai_GaussSeidel(n)
%用块Gauss-Sediel迭代法求解正方形域上的Poisson方程边值问题
%输入n,对x、y轴进行n等分；先确定场u的边界及场源b；
%用k和u分别存放迭代次数和精确解
function [k,u]=kuai_GaussSeidel(n) %对x、y轴进行n等分
h=1/n; %步长
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=zerosn+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
A=zeros(n-1,n-1); %定义矩阵的子块A
for i=1:n-1
if i>1, A(i,i-1)=-1; end
if i %误差
for j=2:n %进行迭代循环
ub=u; %ub是第k-1次迭代结果，用于和第k次迭代结果比较
if j==2, u(2:n,2)=pinv(A)*(u(2:n,3)+b(2:n,2)); end
if j==n, u(2:n,n)=pinv(A)*(u(2:n,n-1)+b(2:n,n)); end
if j>2&jA)*(u(2:n,j-1)+u(2:n,j+1)+b(2:n,j));end
end
error=max(max(abs(u-ub))); %error是前后两次迭代结果的对应元素的最大误差
if error<=e, break; end %判断误差是否达到精度
end
t=toc %统计程序运算时间

4. 计算结果：
>> format short
>> n=10;
>> [k,u]=kuai_GaussSeidel(n)
t =
1.3280
k =
93
u =
Columns 1 through 8
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.0183 1.0194 1.0189
1.0000 1.0047 1.0091 1.0133 1.0168 1.0194 1.0208 1.0203
1.0000 1.0045 1.0088 1.0128 1.0162 1.0189 1.0203 1.0201
1.0000 1.0037 1.0073 1.0107 1.0136 1.0160 1.0174 1.0175
1.0000 1.0023 1.0045 1.0066 1.0084 1.0100 1.0110 1.0113
1.0000 1.0000 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)

图二 块Gauss-Sediel迭代法

三 （预条件）共轭斜量法求解线性方程组Au=f。
1.基本原理
（1）预条件共轭斜量法原理

（2）预优矩阵的选取

2. 算法

3.程序
%用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=2n+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 i2&jnd
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&izz=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
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; %beita是
p=z+beita*p; %对p进行更新，确定搜索方向
error=sqrt(norm(u-ub)); %error是统计误差
if error<=e, break; end %