matlab上机作业报告计算初等反射阵用Householder变换法对矩阵A作正交分解连续函数最佳平方逼近等Word文档格式.docx
- 文档编号:612217
- 上传时间:2023-04-29
- 格式:DOCX
- 页数:36
- 大小:444.58KB
matlab上机作业报告计算初等反射阵用Householder变换法对矩阵A作正交分解连续函数最佳平方逼近等Word文档格式.docx
《matlab上机作业报告计算初等反射阵用Householder变换法对矩阵A作正交分解连续函数最佳平方逼近等Word文档格式.docx》由会员分享,可在线阅读,更多相关《matlab上机作业报告计算初等反射阵用Householder变换法对矩阵A作正交分解连续函数最佳平方逼近等Word文档格式.docx(36页珍藏版)》请在冰点文库上搜索。
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,%程序功能:
n维向量x,数k;
H。
H是Householder初等变换阵,H*x=y,使得y的第k+1项到最后项全为零;
%引用函数:
holder2;
%得到n维向量x的维数;
ifk>
n%如果k值溢出,报错;
k>
n'
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)=I-p\u*u'
5.使用示例:
情形1:
X为零向量
>
x=[0,0,0,0]'
H=holderk(x,1)
M=0
H=
1000
0100
0010
0001
情形2:
K值溢出:
x=[1,2,3,4]'
H=holderk(x,5)
n
情形3:
K值为1:
x=[2,3,4,5]'
-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
-7.3485
0.0000
情形4:
(1)K值为3:
x=[4,3,2,1]'
H=holderk(x,3)
1.0000000
01.000000
00-0.8944-0.4472
00-0.44720.8944
-1
4.0000
3.0000
-2.2361
0
(2)K值为2:
H=holderk(x,2)
0-0.8018-0.5345-0.2673
0-0.53450.8414-0.0793
0-0.2673-0.07930.9604
-3.7417
二、设A为n阶矩阵,编写用Householder变换法对矩阵A作正交分解的程序。
给定n阶矩阵A,通过本程序用Householder变换法对矩阵A作正交分解,得出A=QR
任一实列满秩的m×
n矩阵A,可以分解成两个矩阵的乘积,即A=QR,其中Q是具有法正交列向量的m×
n矩阵,R是非奇异的n阶上三角阵。
(1)输入n阶矩阵A
(2)对
,求Househoulder初等反射阵的
。
(3)计算上三角阵R,仍然存储在A
(4)计算正交阵Q
(5)按要求输出,结束
A-输入的n阶矩阵,同时用于存储上三角阵R;
n-矩(方)阵A的阶数;
Q-Q是具有法正交列向量的n阶矩阵;
p,u-向量A(k:
n,k),对应初等反射阵的ρ,u
k,jj,ii-循环变量;
t1-计算上三角阵R的系数tj;
t2-计算正交矩阵Q的系数ti;
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:
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'
(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
1.00002.00003.0000
2.00003.00000.0000
3.00004.00005.0000
(2)A为4阶矩阵:
A=[1234;
2301;
3456;
1680]
1234
2301
3456
1680
-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
-3.8730-6.7132-6.7132-6.1968
04.46476.4805-1.4783
0-0.0000-3.3070-3.0178
00.00000-1.8187
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<
e,则停机,输出结果u,k.
第三步:
比较不同的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为精度
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<
abs(du),error=abs(du);
end;
%统计误差
=e,break;
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
n+1
b(i,j)=(i-1)*(j-1)*h^2;
end
uu=zeros(n+1,n+1);
u=zeros(n+1,n+1);
%对u赋初值
u(1,1:
n+1)=1;
u(n+1,1:
u(1:
n+1,1)=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
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
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
contourf(uu,'
DisplayName'
'
uu'
figure(gcf)
图一超松弛
二用块Jacobi迭代法求解线性方程组Au=f。
对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表示成块状求解:
给场值u和场源b赋初值,及定义
用公式
,进行迭代计算
第三步:
把第k次的u赋给ub,即ub=u;
然后把第k+1次的u和ub进行比较,看是否达到精度,如果达到精度,则输出迭代次数k和精确解u。
3.程序
[k,u]=kuai_GaussSeidel(n)
%用块Gauss-Sediel迭代法求解正方形域上的Poisson方程边值问题
先确定场u的边界及场源b;
%用k和u分别存放迭代次数和精确解
function[k,u]=kuai_GaussSeidel(n)%对x、y轴进行n等分
%步长
u=zeros(n+1,n+1);
b=h^2*b;
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);
A=zeros(n-1,n-1);
%定义矩阵的子块A
ifi>
1,A(i,i-1)=-1;
ifi<
n-1,A(i,i+1)=-1;
A(i,i)=4;
e=0.000000001;
%e是精度
1000%k是迭代次数
%误差
n%进行迭代循环
ub=u;
%ub是第k-1次迭代结果,用于和第k次迭代结果比较
ifj==2,u(2:
n,2)=pinv(A)*(u(2:
n,3)+b(2:
n,2));
ifj==n,u(2:
n,n)=pinv(A)*(u(2:
n,n-1)+b(2:
n,n));
ifj>
2&
j<
n,u(2:
n,j)=pinv(A)*(u(2:
n,j-1)+u(2:
n,j+1)+b(2:
n,j));
error=max(max(abs(u-ub)));
%error是前后两次迭代结果的对应元素的最大误差
=e,break;
end%判断误差是否达到精度
4.计算结果:
[k,u]=kuai_GaussSeidel(n)
1.3280
k=
93
u=
contourf(u,'
u'
图二块Gauss-Sediel迭代法
三(预条件)共轭斜量法求解线性方程组Au=f。
1.基本原理
(1)预条件共轭斜量法原理
(2)预优矩阵的选取
2
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 上机 作业 报告 计算 初等 反射 Householder 变换 矩阵 正交 分解 连续函数 最佳 平方 逼近
链接地址:https://www.bingdoc.com/p-612217.html