matlab上机作业报告计算初等反射阵用Householder变换法对矩阵A作正交分解连续函数最佳平方逼近等.docx
- 文档编号:30778033
- 上传时间:2023-08-23
- 格式:DOCX
- 页数:36
- 大小:442.64KB
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=
检验:
>>det(H)
ans=
>>H*x
ans=
情形4:
(1)K值为3:
>>x=[4,3,2,1]';
>>H=holderk(x,3)
H=
000
000
00
00
检验:
>>det(H)
ans=
-1
>>H*x
ans=
0
(2)K值为2:
>>x=[4,3,2,1]';
>>H=holderk(x,2)
H=
000
0
0
0
>>det(H)
ans=
>>H*x
ans=
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=
r=
0
检验:
>>q*r
ans=
(2)A为4阶矩阵:
>>A=[1234;2301;3456;1680]
A=
1234
2301
3456
1680
>>[q,r]=qrhh(A)
q=
r=
0
0
00
检验:
>>q*r
ans=
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=~,给场值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=;%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=[: : ];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= kk= 48 ww= uu= Columns1through8 Columns9through11 >>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=;%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= k= 93 u= Columns1through8 Columns9through11 >>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 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 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 fori=1: n-1 ifi>1,A(i,i-1)=-1;end ifi A(i,i)=4; end forj=2: n ifj==2,bb(: 1)=A*u(2: n,2)-u(2: n,3);end%bb=A*u ifj==n,bb(: n-1)=A*u(2: n,n)-u(2: n,n-1);end ifj>2&j j-1)=A*u(2: n,j)-u(2: n,j-1)-u(2: n,j+1);end end r=b(2: n,2: n)-bb;%计算r0 fori=1: n-1%计算z0 z(: i)=pinv(A)*r(: i); end p=z;%计算p0 tic; e=;%e是精度 fork=1: 1000 error=0;a=0;aa=0;cc=0;ub=u; fori=1: n-1 aa=aa+dot(z(: i),r(: i)); ifi==1,cc=cc+dot((A*p(: i)-p(: i+1)),p(: i));end ifi>1&i i)-p(: i+1)-p(: i-1)),p(: i));end ifi==n-1,cc=cc+dot((A*p(: i)-p(: i-1)),p(: i));end end a=aa/cc;%a是 ,确定搜索步长 u(2: n,2: n)=u(2: n,2: n)+a*p;%对u进行更新计算 rr=r;zz=z;%rr和zz存储第k-1次的r和z fori=1: n-1 ifi==1,r(: i)=r(: i)-a*(A*p(: i)-p(: i+1));end ifi>1&i i)=r(: i)-a*(A*p(: i)-p(: i+1)-p(: i-1));end ifi==n-1,r(: i)=r(: i)-a*(A*p(: i)-p(: i-1));end end z(: 1: n-1)=pinv(A)*r(: 1: n-1);% kk=0;mm=0; fori=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是统计误差 iferror<=e,break;end%判断误差是否达到精度 end t=toc%统计计算时间 4.计算结果: >>formatshort >>n=10; >>[k,u]=J_CG(n) t= k= 37 u= >>contourf(u,'DisplayName','u');figure(gcf) 图三J-PCG 四三种迭代法效率分析 有场的等值图可以看出,三种迭代方法的结果(达到相同的精度)比较一致,但是J-PCG只需要37次(耗时秒)迭代即可达到比较好的结果;超松弛迭代法需要48次(耗时秒)迭代即可达到满意的结果;块Gauss-Sediel迭代法需要93次(耗时)迭代才到达满意的结果。 所以对此问题来说,超松弛迭代法和J-PCG法的效率要好于块Gauss-Sediel迭代法。 一.任务: 用MATLAB语言编写连续函数最佳平方逼近的算法程序(函数式M文件)。 并用此程序进行数值试验,写出计算实习报告。 二.程序功能要求: 在书中Page355或Page345的程序(见附一)的基础上进行修改,使其更加完善。 要求算法程序可以适应不同的具体函数,具有一定的通用性。 所编程序具有以下功能: 1.用Lengendre多项式做基,并适合于构造任意次数的最佳平方逼近多项式。 可利用递推关系 2.被逼近函数f(x)不用内联函数构造,而改用M文件建立数学函数。 这样,此程序可通过修改建立数学函数的M文件以适用不同的被逼近函数(要学会用函数句柄)。 3.要考虑一般的情况 。 因此,程序中要有变量代换的功能。 4.计算组合系数时,计算函数的积分采用变步长复化梯形求积法(见附三)。 5.程序中应包括帮助文本和必要的注释语句。 另外,程序中也要有必要的反馈信息。 6.程序输入: (
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 上机 作业 报告 计算 初等 反射 Householder 变换 矩阵 正交 分解 连续函数 最佳 平方 逼近
链接地址:https://www.bdocx.com/doc/30778033.html