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

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

《matlab上机作业报告计算初等反射阵用Householder变换法对矩阵A作正交分解连续函数最佳平方逼近等.docx》由会员分享，可在线阅读，更多相关《matlab上机作业报告计算初等反射阵用Householder变换法对矩阵A作正交分解连续函数最佳平方逼近等.docx（36页珍藏版）》请在冰点文库上搜索。

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

一、给定向量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中绝对值最大的一项的绝对值；

ifM==0%如果x=0，提示出错,程序终止;

disp（'Error:

M=0'）;

return;

else

x=x/M;%规范化

end;

s=norm（x）;%求x的二范数

ifx

（1）<0%首项为负，s值要变号

s=-s;

end

u=x;%除首项外，其余各项x,u相同

u

（1）=s+x

（1）;%计算u的首项

p=s*u

（1）;%计算p

ifn==1u=0;end%若x是1×1维向量，则u=0

（2）

functionH=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的维数；

ifk>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=

1000

0100

0010

0001

情形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.40820.8690-0.1747-0.2184

-0.5443-0.17470.7671-0.2911

-0.6804-0.2184-0.29110.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.0000000

01.000000

00-0.8944-0.4472

00-0.44720.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.0000000

0-0.8018-0.5345-0.2673

0-0.53450.8414-0.0793

0-0.2673-0.07930.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；

fork=1:

n-1

[p,u]=holder2（A（k:

n,k））;%向量A（k:

n,k），对应初等反射阵的ρ,u

forjj=k:

n%计算上三角阵R（仍存贮于A）

t1=dot（u,A（k:

n,jj））/p;%利用向量内积求和

A（k:

n,jj）=A（k:

n,jj）-t1*u;

end

forii=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=[123;230;345]

A=

123

230

345

>>[q,r]=qrhh（A）

q=

-0.26730.87290.4082

-0.53450.2182-0.8165

-0.8018-0.43640.4082

r=

-3.7417-5.3452-4.8107

00.65470.4364

-0.00000.00003.2660

检验：

>>q*r

ans=

1.00002.00003.0000

2.00003.00000.0000

3.00004.00005.0000

（2）A为4阶矩阵：

>>A=[1234;2301;3456;1680]

A=

1234

2301

3456

1680

>>[q,r]=qrhh（A）

q=

-0.25820.0597-0.2660-0.9268

-0.5164-0.10450.8434-0.1049

-0.7746-0.2688-0.46620.3323

-0.25820.9556-0.02220.1399

r=

-3.8730-6.7132-6.7132-6.1968

04.46476.4805-1.4783

0-0.0000-3.3070-3.0178

00.00000-1.8187

检验：

>>q*r

ans=

1.00002.00003.00004.0000

2.00003.0000-0.00001.0000

3.00004.00005.00006.0000

1.00006.00008.00000

数值求解正方形域上的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

第三步：

比较不同的w的迭代次数，用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为精度

fork=1:

N%k为记录迭代次数

error=0;

forj=2:

m-1

forjj=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

iferror

end

end

iferror<=e,break;end%判断是否达到精度

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;

fori=2:

n+1

forj=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计算

fori=1:

m%用不同的w计算迭代

[u,k]=SOR（mu,b,w（i））;%调用[u,k]=SOR（u,b,w），返回迭代次数及精确解

ifkk>k,kk=k;ww=w（i）;uu=u;end%把最少迭代次数付给kk,及其w,u赋给ww,uu

end

t=toc%统计程序运算时间

4.计算结果：

>>formatshort

>>n=10;

>>[kk,ww,uu]=SOR_5dianchafen（n）

t=

0.0310

kk=

48

ww=

1.6000

uu=

Columns1through8

1.00001.00001.00001.00001.00001.00001.00001.0000

1.00001.00111.00221.00311.00391.00441.00471.0045

1.00001.00221.00421.00611.00761.00871.00911.0088

1.00001.00311.00611.00881.01101.01261.01331.0128

1.00001.00391.00761.01101.01381.01591.01681.0162

1.00001.00441.00871.01261.01591.01831.01941.0189

1.00001.00471.00911.01331.01681.01941.02081.0203

1.00001.00451.00881.01281.01621.01891.02031.0201

1.00001.00371.00731.01071.01361.01601.01741.0175

1.00001.00231.00451.00661.00841.01001.01101.0113

1.00001.00001.00001.00001.00001.00001.00001.0000

Columns9through11

1.00001.00001.0000

1.00371.00231.0000

1.00731.00451.0000

1.01071.00661.0000

1.01361.00841.0000

1.01601.01001.0000

1.01741.01101.0000

1.01751.01131.0000

1.01551.01031.0000

1.01031.00721.0000

1.00001.00001.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=zeros（n+1,n+1）;%计算场源b

fori=2:

n+1

forj=2:

n+1

b（i,j）=（i-1）*（j-1）*h^2;

end

end

b=h^2*b;

fori=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

fori=1:

n-1

ifi>1,A（i,i-1）=-1;end

ifi

A（i,i）=4;

end

e=0.000000001;%e是精度

tic;

fork=1:

1000%k是迭代次数

error=0;%误差

forj=2:

n%进行迭代循环

ub=u;%ub是第k-1次迭代结果，用于和第k次迭代结果比较

ifj==2,u（2:

n,2）=pinv（A）*（u（2:

n,3）+b（2:

n,2））;end

ifj==n,u（2:

n,n）=pinv（A）*（u（2:

n,n-1）+b（2:

n,n））;end

ifj>2&j

n,j）=pinv（A）*（u（2:

n,j-1）+u（2:

n,j+1）+b（2:

n,j））;end

end

error=max（max（abs（u-ub）））;%error是前后两次迭代结果的对应元素的最大误差

iferror<=e,break;end%判断误差是否达到精度

end

t=toc%统计程序运算时间

4.计算结果：

>>formatshort

>>n=10;

>>[k,u]=kuai_GaussSeidel（n）

t=

1.3280

k=

93

u=

Columns1through8

1.00001.00001.00001.00001.00001.00001.00001.0000

1.00001.00111.00221.00311.00391.00441.00471.0045

1.00001.00221.00421.00611.00761.00871.00911.0088

1.00001.00311.00611.00881.01101.01261.01331.0128

1.00001.00391.00761.01101.01381.01591.01681.0162

1.00001.00441.00871.01261.01591.01831.01941.0189

1.00001.00471.00911.01331.01681.01941.02081.0203

1.00001.00451.00881.01281.01621.01891.02031.0201

1.00001.00371.00731.01071.01361.01601.01741.0175

1.00001.00231.00451.00661.00841.01001.01101.0113

1.00001.00001.00001.00001.00001.00001.00001.0000

Columns9through11

1.00001.00001.0000

1.00371.00231.0000

1.00731.00451.0000

1.01071.00661.0000

1.01361.00841.0000

1.01601.01001.0000

1.01741.01101.0000

1.01751.01131.0000

1.01551.01031.0000

1.01031.00721.0000

1.00001.00001.0000

>>contourf（u,'DisplayName','u'）;figure（gcf）

图二块Gauss-Sediel迭代法

三（预条件）共轭斜量法求解线性方程组Au=f。

1.基本原理

（1）预条件共轭斜量法原理

（2）预优矩阵的选取

2