数值分析上机报告
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
function [H,y]=holder1(x,k) n=length(x); %取向量长度 M=max(abs(x)); if k>=n, disp('k>=n'); return; %当k大于等于n转出停机 end; if M==0, disp('M=0'); return; %当M等于零转出停机 else z=x/M; %将x规范化 end; s=norm(z); %取z的二范数 if z(k)<0 s=-s; end; p=s*(s+z(k)); u=z; for i=1:k-1 u(i)=0; end; u(k)= s+z(k); H=eye(n,n)-p\u*u'; %求出Householder矩阵 y=H*x;
function[Q,A]=juzhenfenjie(A) [n,n]=size(A); %求矩（方）阵A的阶数 Q=eye(n); %构造正交矩阵Q(1)=I for k=1:n-1 [p,u]=holder2(A(k:n,k)); for jj=k:n %计算上三角阵R,仍内存与A t1=dot(u,A(k:n,jj))/p; %利用向量内积求和 A(k:n,jj)=A(k:n,jj)-t1*u; end for ii=1:n %计算正交矩阵Q t2=dot(u,Q(ii,k:n))/p; Q(ii,k:n)=Q(ii,k:n)-t2*u'; end end
5. 实例示范
> A=[2,-1,7;0,3,10;0,4,5] A = 2 -1 7 0 3 10 0 4 5 >> [Q,A]=juzhenfenjie(A) Q = -1.0000 0 0 0 -0.6000 -0.8000 0 -0.8000 0.6000 A = -2 1 -7 0 -5 -10 0 0 -5
的n阶上三角阵。 3. 变量说明： A- 输入的n阶矩阵，同时用于存储上三角阵R； n-矩（方）阵A的阶数； Q-是具有法正交列向量的n阶矩阵； p，u-向量A（k:n,k），对应初等反射阵的p，u k，jj，ii-循环变量； t1-计算上三角阵R的系数tj； t2-计算正交矩阵Q的系数ti； 4. 程序代码： （1）函数1的程序
（2） 函数2的程序
function [p,u]=holder2(x) n=length(x); %取向量的维数 p=1;u=0; M=max(abs(x)); %求向量x的无穷范数 if M==0, disp('M=0'); return; %如果x=0，终止程序 else x=x/M; %规范化 end; s=norm(x); %求x的二范数 if x(1)<0 s=-s; end; u=x; u(1)= s+x(1); %计算Householder向量u p=s*u(1); if n==1 u=0;end %若x是1*1维向量，则u=0
for i=k+1:n y(i)=0; end;
5. 实例示范： 当x=，k=2时，有 >> x=[4;5;8;6;5] x= 4 5 8 6 5 >> k=2 k= 2 >> [H,y]=holder1(x,k) H= 1.0000 0 0 0 0 0 -0.3881 -0.6209 -0.4657 0 -0.6209 0.7222 -0.2083 0 -0.4657 -0.2083 0.8438 0 -0.3881 -0.1736 -0.1302 y= 4.0000 -11.6423 0 0 0
-0.3881 -0.1736 -0.1302 0.8915
2、 设A为n阶矩阵，编写用Householder变换法对矩 阵A作正交分解的程序。
1. 程序功能： 给定n阶矩阵A，通过本程序用Householder变换法对矩阵A 作正交分解，得出A=QR。 2. 基本原理： 任一实列满值的m*n矩阵A,可以分源自文库成两个矩阵的乘积，即 A=QR，其中Q是具有法正交列向量的m*n矩阵，R是非奇异
1、 给定向量x≠0，计算Householder变换阵，使x后nk个分量化为零。
1. 程序功能： 给定向量x≠0，计算Householder变换阵，使x后n-k个分量化为 零。 2. 方法基本原理： 若x=的分量不全为零，则由Householder变换的基本公式进行编 写程序，当输入K值时要判断与n的大小比较，当在1和n之间时 才可运行程序，否则转出停机。 3. 变量说明： x-输入的n维向量 n-向量x的维数 M-向量x中绝对值最大一项的绝对值 P-Householder初等变换阵的系数 u-Householder初等变换阵的向量 s-向量x的二范数 k-要求输入的k，使得y的第k+1项到最后的项全为零 4. 程序代码
function[Q,A]=juzhenfenjie(A) [n,n]=size(A); %求矩（方）阵A的阶数 Q=eye(n); %构造正交矩阵Q(1)=I for k=1:n-1 [p,u]=holder2(A(k:n,k)); for jj=k:n %计算上三角阵R,仍内存与A t1=dot(u,A(k:n,jj))/p; %利用向量内积求和 A(k:n,jj)=A(k:n,jj)-t1*u; end for ii=1:n %计算正交矩阵Q t2=dot(u,Q(ii,k:n))/p; Q(ii,k:n)=Q(ii,k:n)-t2*u'; end end
5. 实例示范
> A=[2,-1,7;0,3,10;0,4,5] A = 2 -1 7 0 3 10 0 4 5 >> [Q,A]=juzhenfenjie(A) Q = -1.0000 0 0 0 -0.6000 -0.8000 0 -0.8000 0.6000 A = -2 1 -7 0 -5 -10 0 0 -5
的n阶上三角阵。 3. 变量说明： A- 输入的n阶矩阵，同时用于存储上三角阵R； n-矩（方）阵A的阶数； Q-是具有法正交列向量的n阶矩阵； p，u-向量A（k:n,k），对应初等反射阵的p，u k，jj，ii-循环变量； t1-计算上三角阵R的系数tj； t2-计算正交矩阵Q的系数ti； 4. 程序代码： （1）函数1的程序
（2） 函数2的程序
function [p,u]=holder2(x) n=length(x); %取向量的维数 p=1;u=0; M=max(abs(x)); %求向量x的无穷范数 if M==0, disp('M=0'); return; %如果x=0，终止程序 else x=x/M; %规范化 end; s=norm(x); %求x的二范数 if x(1)<0 s=-s; end; u=x; u(1)= s+x(1); %计算Householder向量u p=s*u(1); if n==1 u=0;end %若x是1*1维向量，则u=0
for i=k+1:n y(i)=0; end;
5. 实例示范： 当x=，k=2时，有 >> x=[4;5;8;6;5] x= 4 5 8 6 5 >> k=2 k= 2 >> [H,y]=holder1(x,k) H= 1.0000 0 0 0 0 0 -0.3881 -0.6209 -0.4657 0 -0.6209 0.7222 -0.2083 0 -0.4657 -0.2083 0.8438 0 -0.3881 -0.1736 -0.1302 y= 4.0000 -11.6423 0 0 0
-0.3881 -0.1736 -0.1302 0.8915
2、 设A为n阶矩阵，编写用Householder变换法对矩 阵A作正交分解的程序。
1. 程序功能： 给定n阶矩阵A，通过本程序用Householder变换法对矩阵A 作正交分解，得出A=QR。 2. 基本原理： 任一实列满值的m*n矩阵A,可以分源自文库成两个矩阵的乘积，即 A=QR，其中Q是具有法正交列向量的m*n矩阵，R是非奇异
1、 给定向量x≠0，计算Householder变换阵，使x后nk个分量化为零。
1. 程序功能： 给定向量x≠0，计算Householder变换阵，使x后n-k个分量化为 零。 2. 方法基本原理： 若x=的分量不全为零，则由Householder变换的基本公式进行编 写程序，当输入K值时要判断与n的大小比较，当在1和n之间时 才可运行程序，否则转出停机。 3. 变量说明： x-输入的n维向量 n-向量x的维数 M-向量x中绝对值最大一项的绝对值 P-Householder初等变换阵的系数 u-Householder初等变换阵的向量 s-向量x的二范数 k-要求输入的k，使得y的第k+1项到最后的项全为零 4. 程序代码