基于MATLAB的数值分析编程上机作业

基于MATLAB的数值分析编程上机作业（1）班级：学号：姓名：日期：页数：基于MATLAB的数值分析编程上机作业（1）1一、给定向量x≠0，计算初等反射阵H。1、基本原理：若xxxRx,,,的分量不全为零，则由12112212122()x(,,,)1()22nTTsignxexxxxuxuuuHIIuuu确定的镜面反射阵H使得yeHx；当(1)kn时，由21/2k()T1()()()k1()()()(())(0,,0,,,,)1()()=()2()nkiikknkkknkTkkTkkkkkkTkksignxxxxxxuRuuuHIuu有T121(,,,,,0,,0)nkkkxxxHxR算法1：1、输入x2、将x规范化，xxxM,,,max如果M＝0，则转出停机否则niMxxii,,2,1,3、计算2x，如果01x，则4、)(1x5、计算1,(1)xuxu6、1THIuu7、(M,0,,0)y8、按要求输出结束基于MATLAB的数值分析编程上机作业（1）22、程序代码：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。%使用举例：%[p,u]=holder2(x)%Definevariables:%x－输入的n维向量；%n－n维向量x的维数；%M－M是向量x的无穷范数，即x中绝对值最大的一项的绝对值；%p－Householder初等变换阵的系数ρ；%u－Householder初等变换阵的向量U%s－向量x的二范数；n=length(x);%得到n维向量x的维数；p=1;u=0;%初始化p,u；M=max(abs(x));%得到向量x的无穷范数，即x中绝对值最大的一项的绝对值；ifM==0%如果x=0，提示出错,程序终止;disp('M=0');return;elsex=x/M;%规范化end;s=norm(x);%求x的二范数ifx(1)0%首项为负，s值要变号s=-s;endu=x;%除首项外，其余各项x,u相同u(1)=s+x(1);%计算u的首项p=s*u(1);%计算pifn==1u=0;end%若x是1×1维向量，则u=0functionH=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；示例[p,u]=holder2(x);基于MATLAB的数值分析编程上机作业（1）3%使用举例：%H=holderk(x,k)%Definevariables:%x－输入的n维向量；%n－n维向量x的维数；%p－Householder初等变换阵的系数ρ；%u－Householder初等变换阵的向量U%k－数k，H*x=y，使得y的第k+1项到最后项全为零；n=length(x);%得到n维向量x的维数；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'；3、计算实例：x=[2,0,2,1]'x=2021H=holderk(x,3)H=1.000000001.00000000-0.8944-0.447200-0.44720.8944H*xans=2.00000-2.23610基于MATLAB的数值分析编程上机作业（1）4二、用Householder变换法求矩阵A的正交分解A＝QR。1、基本原理：任一实列满秩的m×n矩阵A，可以分解成两个矩阵的乘积，即A＝QR，其中Q是具有法正交列向量的m×n矩阵，R是非奇异的n阶上三角阵。算法：1、输入n阶矩阵A2、对:,Aknk，求Househoulder初等反射阵的)(,kku。1,,2,1nk3、计算上三角阵R，仍然存储在A)1()(kkkAAHnkkj,,1,kijkijkijnkikijkikjutaaaut)()1()(14、计算正交阵Q)1()()1(kkkQHQIQni,,2,1kjikijkijnkjkjkijkiutqqnkkjuqt)()1()(,,1,15、按要求输出结束2、程序代码：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);%使用举例：%[Q,R]=qrhh(A)%变量说明：%A－输入的n阶矩阵，同时用于存储上三角阵R；%n－矩（方）阵A的阶数；基于MATLAB的数值分析编程上机作业（1）5%Q－Q是具有法正交列向量的n阶矩阵；%p,u－向量A（k:n,k），对应初等反射阵的ρ,u%k，jj，ii－循环变量；%t1－计算上三角阵R的系数tj；%t2－计算正交矩阵Q的系数ti；[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），对应初等反射阵的ρ,uforjj=k:n%计算上三角阵R（仍存贮于A）t1=dot(u,A(k:n,jj))/p;%利用向量内积求和A(k:n,jj)=A(k:n,jj)-t1*u;endforii=1:n%计算正交矩阵Qt2=dot(u,Q(ii,k:n))/p;%利用向量内积求和Q(ii,k:n)=Q(ii,k:n)-t2*u';endend3、计算实例：A=[210;131;014]A=210131014[q,r]=qrhh(A)q=-0.89440.40820.1826-0.4472-0.8165-0.36510-0.40820.9129r=-2.2361-2.2361-0.44720-2.4495-2.44950-0.00003.2863q*rans=2.00001.0000-0.00001.00003.00001.000001.00004.0000