matlab上机作业报告计算初等反射阵用Householder变换法对矩阵A作正交分解连续函数最佳平方逼近等.docx
《matlab上机作业报告计算初等反射阵用Householder变换法对矩阵A作正交分解连续函数最佳平方逼近等.docx》由会员分享,可在线阅读,更多相关《matlab上机作业报告计算初等反射阵用Householder变换法对矩阵A作正交分解连续函数最佳平方逼近等.docx(36页珍藏版)》请在冰点文库上搜索。
![matlab上机作业报告计算初等反射阵用Householder变换法对矩阵A作正交分解连续函数最佳平方逼近等.docx](https://file1.bingdoc.com/fileroot1/2023-6/14/6c01ada9-5deb-4d59-8eb4-cd3d045897e3/6c01ada9-5deb-4d59-8eb4-cd3d045897e31.gif)
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
iferrorend
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
ifiA(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&jn,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
ifiA(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&jj-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&ii)-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&ii)=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.程序输入:
(