matlab上机作业报告计算初等反射阵用Householder变换法对矩阵A作正交分解连续函数最佳平方逼近等.docx
- 文档编号:23154327
- 上传时间:2023-05-08
- 格式:DOCX
- 页数:36
- 大小:444.58KB
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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 上机 作业 报告 计算 初等 反射 Householder 变换 矩阵 正交 分解 连续函数 最佳 平方 逼近
![提示](https://static.bdocx.com/images/bang_tan.gif)
链接地址:https://www.bdocx.com/doc/23154327.html