上机矩阵实验
大连理工线代上机实验
基本操作
四则运算、转置、求逆、求秩、求行列式、组合、 化为行最简形、求特征值
常见任务
① 矩阵的赋值和其加、减、乘、除（求逆）命令； ② 矩阵化为最简行阶梯型的计算命令；[U0,ip]=rref(A) ③ 多元线性方程组MATLAB求解的几种方法；x=inv(A)*b, U=rref(A) ④ 行列式的几种计算机求解方法； D=det(A),[L,U]=lu(A);D=prod(diag(L)) ⑤ n个m维向量组的相关性及其秩的计算方法和命令; r=rank(A),U=rref(A) ⑥ 求欠定线性方程组的基础解系及超定方程解的MATLAB 命令；xb=null(A) ⑦ 矩阵的特征方程、特征根和特征向量的计算命令； f=poly(A);[P,D]=eig(A) ⑧ 化二次型为标准型的MATLAB命令；yTDy=xTAx; 其中 y=P-1x,
• • • •
-66.5556 25.6667 -18.7778 26.5556
• >>
例三、求秩
• • • • • • • • • • • • >> A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; >> r=rank(A); >> r= % = 计算机不显示r的值 ??? r= | Error: Expression or statement is incomplete or incorrect. >> rank(A) ans = 4 >> r r= 4 %不打；则计算机将显示rank(A)的值
例1 用直接解法求解下列线性方程组. 命令如下： A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; x=A\b
数学的上机实验报告
实验题目：线性代数求解方程组一、实验目的1. 理解线性代数中方程组的求解方法。
2. 掌握利用计算机求解线性方程组的算法。
3. 熟悉数学软件（如MATLAB、Python等）在数学问题中的应用。
二、实验内容本次实验主要利用数学软件求解线性方程组。
线性方程组是线性代数中的一个基本问题，其求解方法有很多种，如高斯消元法、矩阵求逆法等。
本实验以高斯消元法为例，利用MATLAB软件求解线性方程组。
三、实验步骤1. 编写高斯消元法算法程序。
2. 输入方程组的系数矩阵和常数项。
3. 调用程序求解方程组。
4. 输出解向量。
四、实验代码及分析1. 高斯消元法算法程序```matlabfunction x = gaussElimination(A, b)[n, m] = size(A);assert(n == m, 'The matrix A must be square.');assert(n == length(b), 'The length of b must be equal to the number of rows in A.');% 初始化解向量x = zeros(n, 1);% 高斯消元for i = 1:n-1% 寻找最大元素[~, maxIdx] = max(abs(A(i:n, i)));maxIdx = maxIdx + i - 1;% 交换行A([i, maxIdx], :) = A([maxIdx, i], :);b([i, maxIdx]) = b([maxIdx, i]);% 消元for j = i+1:nfactor = A(j, i) / A(i, i);A(j, i:n) = A(j, i:n) - factor A(i, i:n); b(j) = b(j) - factor b(i);endend% 回代求解for i = n:-1:1x(i) = (b(i) - A(i, i+1:n) x(i+1:n)) / A(i, i); endend```2. 输入方程组的系数矩阵和常数项```matlabA = [2, 1, -1; 1, 2, 1; -1, 1, 2];b = [8; 5; 2];```3. 调用程序求解方程组```matlabx = gaussElimination(A, b);```4. 输出解向量```matlabdisp('解向量为：');disp(x);```五、实验结果与分析实验结果：```解向量为：2-13```实验分析：通过高斯消元法，我们成功求解了给定的线性方程组。
天津工业大学matlab上机实验题
“MATLAB及其在通信中的应用”上机实验二1——矩阵操作进阶及图形绘制1.利用基本矩阵产生5*4的单位阵、全1阵、全0阵、[2 5]之间均匀分布随机矩阵。
，正态分布随机矩阵（均值1，方差为9）。
将所编程序截图，填入到下面空白处。
（10分）答案:2.矩阵运算：已知，求矩阵c，MATLAB命令窗口中编写语句，实现下述功能，将所编程序截图，填入到下面空白处。
（1）对c中所有元素求和赋值给SUM。
（2）先令矩阵d=c，仅保留d中主对角线上方第1条对角线元素，d中其他元素赋值0。
（3+7=10分）1完成后以附件形式发送到邮箱commu_matlab2014@。
邮件主题为“班级”加“下划线”加“姓名”3. 编写.m（1） 求出方程组的解，并保存在向量 B（2） 将B 的转置作为行向量，对该行向量进行复制拼接，生成6*4阶矩阵的H（3） 对H 矩阵进行重排，生成3*8阶矩阵K ，删除矩阵K 的第3列，生成3*7阶矩阵J 。
将所编程序和矩阵B 、H 、J 的结果截图，填入到下面空白处。
（6+6+8=20分）答案：4、假设有函数22a b ea --=+，其中1,1.5,2,2.5...,6a =1*2排布的图形子窗口，按要求完成绘图。
子窗口1绘制图1： x 轴用a 为坐标刻度，y 轴用b 为坐标刻度，曲线颜色为绿色，类型为-*，线宽为1，加入栅格线。
加入标题“图1”。
x 轴标注 “a ” ，y 轴标注“b ” 。
x 轴坐标范围[08]，y 轴坐标范围[0 58 ]子窗口2绘制图2： x轴用a 为坐标刻度，y 轴用log 10(b)为坐标刻度，曲线颜色为红色，类型为-+，线宽为1.5，加入栅格线。
加入标题“图2”。
x 轴标注 “a ” ，y 轴标注“log10(b)” 。
x 轴坐标范围[0 8]，y 轴坐标范围[0 100 ]将所编程序和输出图形截图，填入到下面空白处。
（5+5=10分）答案：5、打开 一个图形窗口， 222x 0*2*32180180180ππππ=，，，，...，，绘制如下图所示的三条曲线：（1）曲线1：y=sin(x); 正弦曲线颜色 红色，线型为 * ， 线宽为1 （2）曲线2：y=cos(x); 余弦曲线颜色 蓝色，线型为 -- ，线宽为2（3）曲线3：y=sin(x)+3; 正弦曲线纵轴搬移，曲线颜色 绿色，线型为钻石，线宽为1.5 将所编程序截图，填入到下面空白处。
上机实验1-2-3 - 副本
上机操作 一（1）设A 和B 是两个同维同大小的矩阵，问：1）A*B 和A.*B 的值是否相等？ 2）A./B 和B.\A 的值是否相等？ 3）A/B 和B\A 的值是否相等？4）A/B 和B\A 所代表的数学含义是什么？ （2）1）生成一个10*10的二维随机矩阵A2）将矩阵A 第2—5行中第1，3，5列元素赋给矩阵B 。
3）将矩阵A 的每个元素值加30。
4）求矩阵A 的大小和维数。
6）将矩阵A 转换成5*5*4三维矩阵。
（3）下列命令执行后，L1、L2、L3、L4的值分别是多少/ A=1:9;B=10-A; L1=A==B; L2=A<=5; L3=A>3&A<7; L4=find(A>3&A<7); （4）已知⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=7613870451A ，⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=023352138B 求下列表达式的值：1） A+6B 和A 2-B+I （I 为单位矩阵） 2）A*B ，A .*B 和B*A 3）A/B 和B\A4）[A,B]和 [A([1,3]，：);B^2] （5）已知⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡---=14.35454.9632053256545410778.01023A ，取出其前三行构成矩阵B ，其前两列构成矩阵C ，其右下角3×2子矩阵构成矩阵D ，B 与C 的乘积构成矩阵E ，分别求E<D,E&D,E|D 、～E|～D 和find(A>=10&A<25)（6）使用函数，实现方阵左旋90°或右旋90°的功能。
例如，原矩阵为A,A 左旋后得到B,右旋后得到C 。
⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=129631185210741A ，⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=321654987121110B ，⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=101112789456123C （7）建立一个方阵A ，求A 的逆矩阵和A 的行列式的值，并验证A 与A -1是互逆的。
大连理工大学矩阵与数值分析上机作业
end
case2%2-范数
fori=1:n
s=s+x(i)^2;
end
s=sqrt(s);
caseinf%无穷-范数
s=max(abs(x));
end
计算向量x，y的范数
Test1.m
clearall;
clc;
n1=10;n2=100;n3=1000;
x1=1./[1:n1]';x2=1./[1:n2]';x3=1./[1:n3]';
xlabel('x');ylabel('p(x)');
运行结果：
x=2的邻域：
x =
1.6000 1.8000 2.0000 2.2000 2.4000
相应多项式p值：
p =
1.0e-003 *
-0.2621 -0.0005 0 0.0005 0.2621
p(x)在 [1.95,20.5]上的图像
程序：
[L,U]=LUDe.(A);%LU分解
xLU=U\(L\b)
disp('利用PLU分解方程组的解:');
[P,L,U] =PLUDe.(A);%PLU分解
xPLU=U\(L\(P\b))
%求解A的逆矩阵
disp('A的准确逆矩阵:');
InvA=inv(A)
InvAL=zeros(n);%利用LU分解求A的逆矩阵
0 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0625
0 0 0 0.5000 -0.2500 -0.1250 -0.1250
0 0 0 0 0.5000 -0.2500 -0.2500
大连理工大学矩阵分析matlab上机作业
x(i)=1/i; %按要求给向量 x 赋值，其值递减 end normx1=norm(x,1); %求解向量 x 的 1 范数 normx1 normx2=norm(x,2); %求解向量 x 的 2 范数 normx2 normxinf=norm(x,inf); %求解向量 x 的无穷范数 normxinf normy1=norm(y,1); %求解向量 y 的 1 范数 normy1 normy2=norm(y,2); %求解向量 y 的 2 范数 normy2 normyinf=norm(y,inf); %求解向量 y 的无穷范数 normyinf z1=[normx1,normx2,normxinf]; z2=[normy1,normy2,normyinf]; end
for i=2:n
for j=i:n U(i,j)=A(i,j)-L(i,1:i-1)*U(1:i-1,j);
式
%Doolittle 分解计算上三角矩阵的公
L(j,i)=(A(j,i)-L(j,1:i-1)*U(1:i-1,i))/U(i,i); %Doolittle 分解计算下三角矩 阵的公式
end
1 1 1 ������ x = (1, 2 , 3 , … , ������) ,
������ = (1,2, … , ������)������.
对n = 10，100，1000甚至更大的n计算其范数，你会发现什么结果?你能否修改
你的程序使得计算结果相对精确呢？
1.1 源代码
function [z1,z2]=norm_vector(n) %向量 z1 的值为向量 x 的是三种范数，向量 z2 的值为向量 y 的三 种范数，n 为输入参数
MATLAB上机实验1答案
实验1 Matlab 初步一、问题已知矩阵A 、B 、b 如下：⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡-------------=031948118763812654286174116470561091143A ⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡------=503642237253619129113281510551201187851697236421B []1187531=b应用Matlab 软件进行矩阵输入及各种基本运算。
二、实验目的学会使用Matlab 软件构作已知矩阵对应的行（列）向量组、子矩阵及扩展矩阵，实施矩阵的初等变换及线性无关向量组的正交规范化，确定线性相关相关向量组的一个极大线性无关向量组，且将其余向量用极大线性无关向量组线性表示，并能编辑M 文件来完成所有的实验目的。
三、预备知识1、 线性代数中的矩阵及其初等变换、向量组的线性相关性等知识。
2、 Matlab 软件的相关命令提示如下；（1） 选择A 的第i 行做一个行向量：ai=A(i,:)；（2） 选择A 的第j 行做一个列向量：ai=A(j,:)；（3） 选择A 的某几行、某几列上的交叉元素做A 的子矩阵：A([行号],[列号])；（4） n 阶单位阵：eye(n)；n 阶零矩阵：zeros(n)；（5） 做一个n 维以0或1为元素的索引向量L ，然后取A(:,L)，L 中值为1的对应的列将被取到。
（6） 将非奇异矩阵A 正交规范化，orth(A) ；验证矩阵A 是否为正交阵，只需做A*A'看是否得到单位阵E 。
（7） 两个行向量a1和a2的内积：a1*a2'。
（8） 让A 的第i 行与第j 列互换可用赋值语句：A([i,j],:)=A([j,i],:)；（9）让K乘以A的第i行可用赋值语句：A(i,:)=K*A(i,:)；（10）让A的第i行加上第j行的K倍可用赋值语句：A(i,:)=A(i,:)+K*A(j,:)；（11）求列向量组的A的一个极大线性无关向量组可用命令：rref(A)将A化成阶梯形行的最简形式，其中单位向量对应的列向量即为极大线性无关向量组所含的向量，其它列向量的坐标即为其对应向量用极大线性无关组线性表示的系数。
MATLAB上机实验(4)
第一次上机一、目的通过亲自上机，使同学们巩固近期课程所学到的矩阵初等运算、流程控制以及二/三维绘图等知识。
二、步骤1）给同学们一小段时间，让同学们首先熟悉一下MATLAB 运行环境，包括其桌面环境、菜单、工具栏等。
2）给同学们在黑板上列出每道习题，然后根据同学们在实际解算过程中遇到的问题给予解答，最后将同学们遇到的共性问题在黑板上作出详细的解释。
三、习题1、矩阵的初等运算1）3 2.5347.8 4.32A ⎛⎫= ⎪--⎝⎭，试用两种方法在MATLAB 环境下构造这个矩阵。
2）[0.2 0.45 0.7 0.95 1.2 1.45]I =，试用“:”关键符构造这个向量。
3） 1.54,37.5A B -⎛⎫⎛⎫== ⎪ ⎪-⎝⎭⎝⎭，试通过行列聚合的方式构造两个新矩阵。
4） 2.33.1281.90.35 5.669.018.27A ⎛⎫ ⎪=- ⎪ ⎪--⎝⎭，试将其1~2行、2~3列的元素取出，构造成一个新矩阵。
5）随机生成两个5×5阶的矩阵，求两个矩阵的+、-、×、÷、&、≥操作。
6）0.20.90.730.560.980.010.910.730.27A --⎛⎫ ⎪=- ⎪ ⎪-⎝⎭，试求取该矩阵的正弦。
2、流程控制1）计算1-100所有整数的和。
2）计算前n 个自然数的阶乘，直到超过1030。
3）随机生成一个1×10000阶的向量，判断该向量中元素值等于0.3的元素个数。
3、基本绘图1）绘制函数221y x x =+在定义域0.5~12.3上的图形，并将线型颜色调整为红色虚线，在标题上添加一个合适的标题说明这个曲线。
2）绘制函数2cos()*sin()2x t t y t z t =⎧⎪=+⎨⎪=⎩在定义域t=0~2π上的图形。
3）绘制函数2211z x y =+在定义域x=1~10,y=1~10上的图形。
MATLAB上机实验报告02实验名称：MATLAB绘图和矩阵的初等计算
MATLAB上机实验报告02实验名称：MATLAB绘图和矩阵的初等计算MATLAB上机实验报告02实验名称：MATLAB绘图和矩阵的初等计算平顶山学院计算机语言类课程实验报告课程名称院系学号实验日期MATLAB语言及应用实验机房专业姓名任课教师实验学时2班级机器号实验成绩一．实验名称：MATLAB绘图和矩阵的初等计算二．实验目的和要求1、掌握MATLAB构建组合组合矩阵的方法；2、掌握MATLAB流程控制语句的运用；3、掌握MATLAB绘图方法三．实验内容教材（《MATLAB及其在理工课程中的应用指南，陈怀琛，西安电子科技大学出版社》）P93-3,4四．实验设计方案（实验步骤或开发过程）1、a:使用eye(),magic(),ones(),zeros()这几个函数，直接求出单位矩阵I，魔方矩阵M，全幺矩阵A，全零矩阵Bb:根据题目要求，先计算出C的各个元素，最后再组合为矩阵c:直接从得到的C矩阵中，利用C1=C([2468];:)，再用C2=C1(:;[2468]),分别求出C1和C2d:要求C1和C2的乘积，则首先判断C1的列数和C2的行数是否相同，若相同，则求出D，否则，不能进行计算；C2和C1相乘也按上述所述2、以X在0到2*pi的范围里，插入101个点为横坐标，以y 所得结果为纵坐标，利用plot()函数画出该曲线五．实验中存在问题及解决办法在第三题的b中，[A不能只通过简单的[AB]’就可以得到，必须和I的列数相同才行B]在第四题中，多注意x平方时，因x是行向量，因此采用元素乘方在x后边加"."六．实验结果1、I=[1000M=[162313010051110800109761201*1]414151]B=[0000 A=[11110000]1111]C=[1000110001001100001011000001110011 11162313111151110800009761201*0414151]C1=[010********* 110011115111080000414151]C2=[1010D：C1和C2不能相乘01101111800141]D1=[1211612108111261210811121112892352309614141414 74168155113]2、七．附录（源程序清单）1、第三题clccleardisp第三题aI=eye(4)M=magic(4)A=ones(2,4)B=zeros(2,4)disp第三题bD=[A"B"]E=[A"B"]"C=[ID;EM]disp第三题cC1=C([2468],:)C2=C1(:,[2468])disp第三题d[mC1nC1]=size(C1);[mC2nC2]=size(C2);ifnC1==mC2D=C1*C2 elsedispC1和C2不能相乘endifnC2==mC1D1=C2*C1elsedispC2和C1不能相乘end2、第四题clcclearx=linspace(0,2*pi,101);y=cos(x.*(0.5+3*sin(x)/ (1+x.^2)));%多注意x平方时，因x是行向量，因此采用元素乘方在x后边加"."plot(x,y,"*g")title("y随x的变化曲线")xlabel("X")ylabel("Y")grid扩展阅读：MATLAB上机实验报告03实验名称：MATLAB简单数据处理和构建特殊矩阵平顶山学院计算机语言类课程实验报告课程名称院系学号实验日期MATLAB语言及应用实验机房专业姓名任课教师实验学时2班级机器号实验成绩一．实验名称：MATLAB简单数据处理和构建特殊矩阵二．实验目的和要求1、掌握MATLAB求解统计学问题；2、掌握MATLAB函数创建特殊矩阵三．实验内容教材（《MATLAB及其在理工课程中的应用指南，陈怀琛，西安电子科技大学出版社》）P93-9,10,25四．实验设计方案（实验步骤或开发过程）1、利用randn()函数产生正态分布的随机数矩阵R1，再用mean()和std()函数求出R1的各列的平均值和均方差；求出整体的平均值和均方差，要先对各列的平均数和均方差转置后再使用mean()和std()2、使用rand()函数产生均匀分布的随机数矩阵R,然后抽出该矩阵的前四列，并用inv()函数求出其逆矩阵3、a:根据所给矩阵得到：使用一个简单的行向量，再使用repmat()函数，就可得到该矩阵b:采用两个增量式a和b，a=[1:1:4],b=[0:1:3]，设c=a’,C=[cccc],B=[bbbb]再利用C.^B来获得该矩阵五．实验中存在问题及解决办法在第十题中使用rand函数获得-16到16之间的整数值，不会表示，经过上网查找，可以使用六．实验结果1、R1=-1.06670.18250.09830.23232.02371.00010.9337-1.56510.04140.4264-2.2584-1.66420.3503-0.0845-0.7342-0.37282.2294-0.5900-0.02901.6039-0.0308-0.23650.3376-0.2781a=0.04710.0342-0.15630.01240.5831-0.3831b=0.84151.29840.38880.37862.07521.0967a1=0.0229b 1=0.63852、R=-5-9-10-3-7-14-12-9-5-12-2-1-14-5-9-6-12-6-6-4-3-7-11-8a=-5-9-10-3-12-9-5-12-14-5-9-6-6-4-3-7inva=0.0709-0.1679-0.09990.3430-0.0622-0.31190.07740.4949-0.08250.3075-0.0369-0.46020.01010.19030.0572-0.52243、a:a=[-3-2-10123]A=[-3-2-10123-3-2-10123-3-2-10123-3-2-10123]b:a=[1234]b=[0123]c=[1C=[1111B=[0123s=[111112222201*3124816333330123 13927814]4444]0123]141664256]七．附录（源程序清单）201*.10.11%第9题%clcclearallR1=randn(4,6)a=mean(R1)%求各列的平均数b=std(R1)%求各列的标准差（均方差）a1=mean(mean(R1)")%求整体的平均数（mean2（R1））b1=std(std(R1)")%求整体的标准差（std2（R1））%第10题%clcclearallR=round(-16+16*rand(4,6))a=R(:,[1234])inva=inv(a)%第25题a，用行向量表示矩阵%clcclearalla=[-3:1:3]A=repmat(a,4,1)%a 代表要复制的矩阵，4是行数，1是列数%第25题b，用行向量表示该矩阵%clcclearalla=[1:1:4]c=a"b=[0:1:4]C=[ccccc]B=[b;b;b; b]s=C.^B。
上机矩阵实验讲解
矩阵的LU 分解1、原理 ： 设A n nC⨯∈若A 可以表示成一个下三角矩阵和一个上三角矩阵的乘积A=LU ，则称其为矩阵A 的LU 分解（三角分解）。
矩阵的LU 分解在求解线性方程组时将十分简便。
如对线性方程组Ax=b ，设A=LU 是其LU 分解。
我们先求解方程组Ly=b 。
由于L 是下三角矩阵，则解向量y 可以通过依次求出其分量y1，y2，……，yn 而求出，再求解方程组Ux=y 。
解向量x 可以通过该方程组依次求出分量xn ，xn-1，……，x2，x1而快速得出。
于是由两个方程组Ux=y ，Ly=b 的求解而给出LUx=Ly=b=Ax 的解。
若矩阵A 非奇异，则A 能分解为LU 的充分必要条件是A 的顺序主子行列式不为0。
则存在惟一的主对角线上元素全为1的下三角阵L 与惟一的上三角阵U ，使得A=LU 。
2、算法当n 阶矩阵A 满足LU 分解的条件时，我们可以用初等变换的方法求出L 和U.因为当A=LU 时，由于L 可逆，故必存在可逆阵P 组PL=I即 PA=PLU=U.也就是说，我们可以先对A 施行行的初等变换得出上三角阵U,而矩阵P 可以通过对单位阵I 进行相同的行初等变换得出。
即有P(A,I)=(PA,PI)=(U,P)于是A=P^-1U,为保持P 为上三角阵，在进行行初等变换时，不能进行行的对换，上行的倍数应加到下行的对应元。
我们知道对于第i 行j 列元素a[ij]的第k+1次变换，有设A=L*U ，那么有如下数学表达：3、程序如下：A=[;;;];b=[1 1 1 1];n=length(b);%方程个数nx=zeros(n,1);%未知向量A(2:n,1)=A(2:n,1)./A(1,1);for i=2:n-1A(i,i)=A(i,i)-sum(A(i,1:i-1)'.*A(1:i-1,i));for j=i+1:nA(i,j)=A(i,j)-sum(A(i,1:i-1)'.*A(1:i-1,j));A(j,i)=(A(j,i)-sum(A(j,1:i-1)'.*A(1:i-1,i)))/A(i,i);endendA(n,n)=A(n,n)-sum(A(n,1:n-1)'.*A(1:n-1,n));AU=A;L=A;for i=1:nL(i,i)=1;endfor i=1:n-1for j=i+1:nL(i,j)=0;endendL %下三角阵for i=2:nfor j=1:i-1U(i,j)=0;endendU %上三角阵4、例子以矩阵A=[2 1 1；4 1 0；-2 2 1]为例说明上述程序的正确性，在matlab中输入以下指令：5、结果输入上述指令以后，得到如下的结果：此时可以说，上述程序是正确的。
离散数学上机实验报告 离散数学实验报告：建立关系矩阵实验
离散数学上机实验报告离散数学实验报告：建
立关系矩阵实验
__大学离散数学实验报告建立关系矩阵实验姓名：____
专业：
软件工程班级：
3 学号：
1325116025 日期：
20__.10月7日 1、摘要：建立关系矩阵实验的目的是理解并掌握关系的矩阵表示方法、为用序偶集合表示的关系建立相应的关系矩阵。
学会用所学过的程序设计语言编程，解决关系矩阵的自动建立问题。
实验的内容是用二维数组或向量存储关系矩阵，根据输入的用序偶集合表示的关系，建立相应的关系矩阵。
用建立二维数组的方法构造关系矩阵。
分别输入两个用序偶集合表示的关系作为实验数据，然后建立两个数组之间的关系，得到一个关系矩阵。
关系矩阵一开始初始化为0，建立成功的关系体现为1。
最后得到一个完整的矩阵。
一．导言 2、 1) 问题的描述。
实验的目标是如何为用序偶集合表示的关系建立相应的关系矩阵，解决关系矩阵的自动建立问题。
2) 拟采用的方法用建立二维数组的方法来解决建立关系矩阵。
首先建立两个数组分别代表行和列，然后建立一个新的二维数组，将其初始化为零，集合之间的关系对应真值表，所以在这个二维数组中两个集合的关系就被表示为1，然后就相应地建立了两个集合的关系矩阵。
二．实验过程 1) 算法思想流程 1.申请两个字符型数组用来储存集合元素。
2.建立二维数组然后初始化为0.
3.判断关系是否存在，存在则赋值为1。
4输入数值然后输出关系矩阵。
2）程序流程图开始 P=0 switch P=0
P=’a’&&j=’a’&&j=’a’||j。
大连理工大学矩阵与数值分析MATLAB上机实验
二、解线性方程组 1．分别 Jacobi 迭代法和 Gauss-Seidel 迭代法求解线性方程组
3 1 0 0 x1 1 1 3 1 0 x2 0 , 0 1 2 1 x3 0 0 0 1 3 x4 0
Gauss 列主元消去法程序:
clc; clear; format long a=[2,4,3,1;8,2,0,0;5,0,4,0;9,0,0,5]; %系数矩阵 b=[12;6;23;16]; [n,m]=size(a); nb=length(b); det=1; for k=1:n-1 amax=0; for i=k:n if abs(a(i,k))>amax amax=abs(a(i,k)); r=i; end end if amax<1e-10 return; end if r>k for j=k:n z=a(k,j); a(k,j)=a(r,j); a(r,j)=z; end z=b(k);
从小到大求和程序计算结果：
N 100 10000 1000000 从小到大求和程序得 到的 SN 0.497512437810945 0.499975001249937 0.499999750000134 真实值������������ = ������ 2������ + 1 0.497512437810945 0.499975001249937 0.499999750000125 计算值有效位数 15 15 13
8
2
1 dx x
复化梯形公式程序
clc; clear; format long syms t m=int(1/t,2,8); %真实值 a=2; b=8; n=300; h=(b-a)/n; sum=0; f=inline('1/x'); for i=1:n-1 sum=sum+f(a+i.*h); end T=h/2*(f(a)+2*sum+f(b))
C_上机实验(含作业)总的目的、要求和评分标准
上机实验（含作业）总的目的、要求和评分标准一、实验目的实验作为教学的一个重要环节，其目的在于更深入地理解和掌握课程教学中的有关基本概念，应用基本技术解决实际问题，从而进一步提高分析问题和解决问题的能力。
C程序设计课程实践性很强，即要求独立编写程序，学会独立上机调试程序。
学会独立上机调试程序。
也就是要善于发现程序中的错误，并且能很快地排除这些错误，使程序能正确运行。
计算机技术是实践性很强的技术，要求从事这一领域的人不仅能了解和熟悉有关理论和方法，还要求自己动手实践。
对程序设计来说，要求会编写程序并上机调试通过。
因此调试程序本身是程序设计课程的一个重要的内容和基本要求，应给予充分的重视。
调试程序的经验固然可以借鉴他人的现成经验，但更重要的是通过自己的直接实践来积累，而且有些经验是只能“会意”难以“言传”。
因此，在实验时不但要达到通过程序完成每一次的实验任务，而且应当在已通过的程序基础上作进一步的修改、提高和完善。
甚至于“自设障碍”，即把正确的程序改为有错的（如用scanf函数为输入变量输入数据时，漏写“&”符号，double变量使用格式符“%f”；使数组下标出界；使整数溢出等等），观察和分析所出现的情况。
这样的学习才会有真正的收获。
实验目的可归纳如下：⒈验证自己已建立起来的概念或所编写的程序是否正确；⒉加深对课堂所学内容的理解和语法规则的记忆；⒊理解和掌握运用计算机高级语言进行编程的思想方法；⒋掌握常用算法的设计与应用实现；⒌熟悉Turbo C 2.0程序开发环境，掌握C程序常用的调试手段；⒍学会上机调试程序的方法，不断积累调试经验，提高排错能力；⒎使自己具有独立的应用编程和熟练的程序调试能力。
二、要求：⒈做好每一次上机前的准备以提高上机效率：①预先认真阅读相关实验内容，做到心中有明确的目的要求和任务，要有备而来；②按照实验内容规定的习题题目，事先在实验预习报告上编写好源程序及运行程序所需的典型数据，并经人工静态检查认为无误；手编程序应书写整齐，应在每个题目之间留出一定的空间，以备记录上机调试情况和运行结果等；对程序中自己有疑问的地方，应作出记号，以便上机时给以注意。
大连理工大学矩阵与数值分析上机作业13388
共享知识分享快乐大连理工大学矩阵与数值分析上机作业课程名称：矩阵与数值分析研究生姓名：12 交作业日时间：日20 月年2016卑微如蝼蚁、坚强似大象．共享知识分享快乐第1题1.1程序：Clear ;all n=input('请输入向量的长度n：') for i=1:n;v(i)=1/i;endY1=norm(v,1)Y2=norm(v,2)Y3=norm(v,inf)1.2结果n=10 Y1 =2.9290Y2 =1.2449Y3 =1n=100 Y1 =5.1874Y2 =1.2787Y3 =1n=1000 Y1 =7.4855Y2 =1.2822Y3 =1N=10000 Y1 =9.7876Y2 =1.2825Y3 =11.3 分析一范数逐渐递增，随着n的增加，范数的增加速度减小；二范数随着n的增加，逐渐趋于定值，无群范数都是1.卑微如蝼蚁、坚强似大象．共享知识分享快乐第2题2.1程序;clear all x(1)=-10^-15;dx=10^-18;L=2*10^3; i=1:L fory1(i)=log(1+x(i))/x(i); d=1+x(i); d == 1ify2(i)=1;elsey2(i)=log(d)/(d-1);endx(i+1)=x(i)+dx;end x=x(1:length(x)-1););'r'plot(x,y1,on holdplot(x,y2);卑微如蝼蚁、坚强似大象．共享知识分享快乐2.2 结果2.3 分析红色的曲线代表未考虑题中算法时的情况，如果考虑题中的算法则数值大小始终为1，这主要是由于大数加小数的原因。
第3题3.1 程序;clear all A=[1 -18 144 -672 2016 -4032 5376 -4608 2304 -512];x=1.95:0.005:2.05; i=1:length(x);for y1(i)=f(A,x(i)); y2(i)=(x(i)-2)^9;end figure(3);plot(x,y1);;on hold卑微如蝼蚁、坚强似大象．共享知识分享快乐);'r'plot(x,y2,F.m文件y=f(A,x)function y=A(1); i=2:length(A);for y=x*y+A(i);;end3.2 结果第4题卑微如蝼蚁、坚强似大象．共享知识分享快乐4.1 程序;clear all n=input('请输入向量的长度n：')A=2*eye(n)-tril(ones(n,n),0); i=1:n for A(i,n)=1;end n=length(A);U=A; e=eye(n);for i=1:n-1[max_data,max_index]=max(abs(U(i:n,i))); e0=eye(n);max_index=max_index+i-1; U=e0*U; e1=eye(n); j=i+1:n fore1(j,i)=-U(j,i)/U(i,i);endU=e1*U;中把变换矩阵存到P P{i}=e0;% L{i}=e1; e=e1*e0*e;endk=1:n-2for Ldot{k}=L{k}; i=k+1:n-1forLdot{k}=P{i}*Ldot{k}*P{i};endend Ldot{n-1}=L{n-1};LL=eye(n);PP=eye(n); i=1:n-1for PP=P{i}*PP;LL=Ldot{i}*LL;endb=ones(n,2);解方程 %b=e*b;x=zeros(n,1);x(n)=b(n)/U(n,n); i=n-1:-1:1for卑微如蝼蚁、坚强似大象．共享知识分享快乐x(i)=(b(i)-U(i,:)*x)/U(i,i);end计算逆矩阵%X=U^-1*e^-1*eye(n);AN=X'; result2{n-4,1}=AN;result1{n-4,1}=x;,n)'%d:\n'fprintf(fprintf('%d ',AN);4.2 结果n=51.0625 -0.875 -0.75 -0.5 -0.0625-0.0625 0.0625 -0.75 1.125 -0.5-0.0625 0.125 0.0625 1.25 -0.5-0.0625 0.1250.25 0.06251.50.0625-0.5-0.25-0.0625 -0.125n=101.0625 -0.875 -0.75 -0.5 -0.0625 1.0625 -0.875 -0.75 -0.5 -0.0625 -0.0625 1.125 0.0625 -0.75 -0.5 -0.5 0.0625 1.125 -0.75 -0.0625 -0.0625 0.0625 0.125 1.25 1.25 -0.0625 -0.5 0.0625 0.125 -0.5-0.0625 0.250.250.0625 0.1251.5 1.5 -0.0625 0.1250.06250.0625 -0.0625 -0.125 -0.25 0.0625 -0.5 -0.0625 -0.125 -0.25 -0.5 -0.0625 -0.75 1.0625 -0.5 -0.0625 -0.875 -0.5 -0.75 1.0625 -0.875 -0.0625 -0.5 0.0625 1.125 -0.5 0.0625 1.125 -0.75 -0.0625 -0.75 1.25 0.125 0.0625 -0.0625 -0.0625 -0.5 -0.5 0.0625 0.125 1.250.25-0.0625 -0.0625 1.50.1250.0625 0.0625 0.250.1251.5-0.0625 -0.125 -0.25 0.0625-0.5 0.0625 -0.0625 -0.125 -0.25-0.5同样的方法可以算出n=20,n=30时的结果，这里就不罗列了。
上机实验第二讲
• 均接受原假设 所以 b3=0 成立
• 自己联系:将方程2改为 :price=b0+b1*weight+b2*length+b3*mpg+ b4*foreign+ε
异方差得检验与FGLS
• 异方差就是违背了球型扰动项假设得一种 情形。在存在异方差得情况下:
• 打开程序:doedit myprog、ado • 执行MLE:
例一:
sysuse auto,clear
ml model lf myprog (price= weight length foreign)(sigma:)
ml max 和OLS比较:reg price weight length foreign 回归系数完全相同
• 3。BP 检验:做完回归后,使用命令: • estat hettest ,normal(使用拟合值yˆ ) • estat hettest,rhs (使用方程右边得解释变量,而不
就是yˆ ) • estat hettest [varlist] (指定使用某些解释变量) • 最初得BP 检验假设扰动项服从正态分布,有一
• 2、逐个分层加入 • Stepwise,pe(0、05) hier:reg price mpg rep78 headroom
trunk weight length turn displacement gear_ratio foreign
残差点得图形表示
• rvfplot:残差拟合值图
•
可以加参数yline(0)
• FGLS估计 • predict e,resid • gen e2=e^2 • gen lne2=log(e2) • reg lne2 weight,noc • predict lne2f • gen e2f=exp(lne2f) • reg price mpg weight foreign [aw=1/e2f]
上机习题2教案
实验二基本矩阵操作实验目的：①掌握matlab变量和数据操作；②掌握matlab矩阵的创立、拆分及特殊矩阵；③掌握matlab运算，掌握matlab在矩阵分析中的应用，掌握稀疏矩阵的存储方式和创建方法；④了解字符串处理函数。
实验要求:给出程序和实验结果。
实验内容：1、利用列向量()1,2,3,,10T建立一个范得蒙矩阵A，并利用位于矩阵A的奇数行偶数列的元素建立一个新的矩阵B，须保持这些元素的相对位置不变。
2、矩阵的基本运算与点运算的区别。
3、给出矩阵的两种存储方式的联系和区别，这两种存储方式在实际应用中主要应用于具有什么特点的矩阵？4．将字符串'very good'转换为等值的整数。
5．按水平和竖直方向分别合并下述两个矩阵：6．分别删除第5题两个结果的第2行。
7．分别将第5题两个结果的第2行最后3列的数值改为[11 12 13]。
8．分别查看第5题两个结果的各方向长度。
9．分别判断pi是否为字符串和浮点数。
10．分别将第5题两个结果均转换为2⨯9的矩阵。
11．计算第5题矩阵A的转置。
12．分别计算第5题矩阵A和B的A+B、A.*B和A\B。
13．判断第5题矩阵A和B中哪些元素值不小于4。
14．分别用函数strcat()和矩阵合并符合并如下字符串：' The picture is '和' very good '。
15．创建字符串数组，其中元素分别为‘Picture ’和'Pitch '。
16．在第14题结果中查找字符串'e'。
17．在第15题结果中匹配字符串'Pi'。
18．将十进制的50转换为二进制的字符串；将十六进制的字符串‘50’转换为三进制的整数。
实验结果：1. x=1:10x =1 2 3 4 5 6 7 8 9 10A=Vander(x)B=A(1:2:9,2:2:10)2.矩阵的基本运算中A*B，A列长度必须和矩阵B的行长度一致。
琼斯矩阵上机报告
其中， A0 Ay / Ax ， y x 。其归一化的琼斯矢量为：
E
Ax Ax Ay
2 2
1 A exp( i ) 0
我们假定偏振器件对于入射偏振光（用琼斯矢量表示）的作用是一个线性变 换，也就是说，出射光琼斯矢量的两个分量 A2 和 B2 是入射光琼斯矢量的两个分 量 A1 和 B1 的线性组合，即
1.2 琼斯矢量与琼斯矩阵
在涉及偏振光的问题时，采用矩阵运算方法比较简洁、方便。琼斯矩阵最重 要的作用在于计算偏振光通过偏振器件后偏振态的变化。 一个偏振器件的作用在 于对入射偏振光的光矢量进行一个线性变换，因此这种变换可以用矩阵表示。偏
振器件的特性可以用一个 2x 2 矩阵来描述，该矩阵称为偏振器件的琼斯矩阵。 沿 z 轴方向传播的任一理想单色偏振光，不管是线偏振光、圆偏振光还是椭 圆偏振光， 其光矢量都可以分解为光矢量分别沿 x 轴和 y 轴的两束线偏振光，可 表示为
g11 g12 G g 21 g 22
其中，4 个矩阵元 g11 、 g12 、 g 21 、 g 22 一般是复常数，具体形式与坐标系的选择 有关。
2 程序流程图
给定参数，将液晶盒分层
声明得出每一层的琼斯矩 阵，得出近似琼斯矩阵
设定角度两片偏振片之间 的角度， 计算偏振片琼斯矩 阵
A2 g11 A1 g12 B1 B2 g 21 A1 g 22 B1
写成矩阵形式为
A2 g11 g12 A1 B g g ห้องสมุดไป่ตู้ 2 21 22 1
或写成：
Et GEi
其中，G 为偏振器件的琼斯矩阵：
2011059120016 宋正格
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
矩阵的LU 分解1、原理 ： 设A n nC⨯∈若A 可以表示成一个下三角矩阵和一个上三角矩阵的乘积A=LU ，则称其为矩阵A 的LU 分解（三角分解）。
矩阵的LU 分解在求解线性方程组时将十分简便。
如对线性方程组Ax=b ，设A=LU 是其LU 分解。
我们先求解方程组Ly=b 。
由于L 是下三角矩阵，则解向量y 可以通过依次求出其分量y1，y2，……，yn 而求出，再求解方程组Ux=y 。
解向量x 可以通过该方程组依次求出分量xn ，xn-1，……，x2，x1而快速得出。
于是由两个方程组Ux=y ，Ly=b 的求解而给出LUx=Ly=b=Ax 的解。
若矩阵A 非奇异，则A 能分解为LU 的充分必要条件是A 的顺序主子行列式不为0。
则存在惟一的主对角线上元素全为1的下三角阵L 与惟一的上三角阵U ，使得A=LU 。
2、算法当n 阶矩阵A 满足LU 分解的条件时，我们可以用初等变换的方法求出L 和U.因为当A=LU 时，由于L 可逆，故必存在可逆阵P 组PL=I即 PA=PLU=U.也就是说，我们可以先对A 施行行的初等变换得出上三角阵U,而矩阵P 可以通过对单位阵I 进行相同的行初等变换得出。
即有P(A,I)=(PA,PI)=(U,P)于是A=P^-1U,为保持P 为上三角阵，在进行行初等变换时，不能进行行的对换，上行的倍数应加到下行的对应元。
我们知道对于第i 行j 列元素a[ij]的第k+1次变换，有设A=L*U ，那么有如下数学表达：3、程序如下：A=[;;;];b=[1 1 1 1];n=length(b);%方程个数nx=zeros(n,1);%未知向量A(2:n,1)=A(2:n,1)./A(1,1);for i=2:n-1A(i,i)=A(i,i)-sum(A(i,1:i-1)'.*A(1:i-1,i));for j=i+1:nA(i,j)=A(i,j)-sum(A(i,1:i-1)'.*A(1:i-1,j));A(j,i)=(A(j,i)-sum(A(j,1:i-1)'.*A(1:i-1,i)))/A(i,i);endendA(n,n)=A(n,n)-sum(A(n,1:n-1)'.*A(1:n-1,n));AU=A;L=A;for i=1:nL(i,i)=1;endfor i=1:n-1for j=i+1:nL(i,j)=0;endendL %下三角阵for i=2:nfor j=1:i-1U(i,j)=0;endendU %上三角阵4、例子以矩阵A=[2 1 1；4 1 0；-2 2 1]为例说明上述程序的正确性，在matlab中输入以下指令：5、结果输入上述指令以后，得到如下的结果：此时可以说，上述程序是正确的。
矩阵分解成LU形式是有条件的，首先矩阵必须是非奇异的矩阵，其次矩阵的全部顺序主子式非零的时候才能完全保证矩阵可分解成LU且分解唯一。
矩阵的奇异值分解1、原理设A ∈C m ×n ，s 1，s 2，…，s r 是A 的非零奇异值，则存在m 阶酉矩阵U ∈C m ×n 及n 阶酉矩阵V ，m ×n 矩阵D ，D= 10000000000000r s s ⎛⎫⎪⋅⋅⋅⎪ ⎪⎪⎝⎭= 000∑⎛⎫ ⎪⎝⎭使得A=UDV H 这就是矩阵A 的奇异值分解。
2、算法第一步：求出A H A 的特征值1λ≥2λ≥…≥r λ＞0=1r λ+=…=n λ，确定非零奇异值i s =i λ，i=1,2，…，r 。
第二步：分别求出矩阵A H A 的对应于特征值i λ的特征向量并将其单位正交化，得到标准正交向量组α1，α2，…，αn 令V=（α1，α2，…，αn ）=（V 1，V 2），V 1=（α1，α2，…，αr ），V 2=（αr+1，αr+2，…，αn ）第三步：若U=（γ1，γ2，…，γr ，γr+1，γr+2，…，γm ）=（U 1，U 2），其中U 1=（γ1，γ2，…，γr ），U 2=（γr+1，γr+2，…，γm ）， 则因（A α1，A α2，…，A αr ）=（s 1γ1，s 2γ2，…，s r γr ）即有U 1=AV 1 1-∑。
其中1-∑=11121r s s s ---⎛⎫⎪ ⎪ ⎪⋅⋅⋅⎪ ⎪⎝⎭第四步：解方程组AA H y=0，对基础解系单位正交化可以求得γr+1，γr+2，…，γm ，令U=（γ1，γ2，…，γr ，γr+1，γr+2，…，γm ）。
3、程序在matlab 中有求解矩阵奇异值分解的svd 函数，调用格式为[U ，S ，V] = svd(A)，其中U 就是所求的U 矩阵，S 是所求的对角阵，V 就是所求的酉矩阵V 。
另一程序：function [U1,D,V1]=QYF(A)% 函数[L,D,V]=QYF(S)实现对矩阵A 的奇异值分解 [m,n]=size(A); U1=zeros(m); V1=zeros(n); D=zeros(m,n); [V1,DV]=eig(A'*A); [U1,DU]=eig(A*A');DU=fliplr(DU);DU=flipud(DU);U1=fliplr(U1);V1=fliplr(V1);DV=fliplr(DV);DV=flipud(DV);s=size(DU,1);t=size(DV,1);f=min(s,t);if s<tD(1:f,1:f)=DU.^(0.5);elseD(1:f,1:f)=DV.^(0.5);End4、例子以矩阵A=[1 0；0 1；1 0]为例来说明上述求矩阵的奇异值分解的svd函数的正确性。
或利用第二程序，输入矩阵A，在matlab中输入下述指令：5、结果输入上述指令以后，得到如下的结果：这与教材上给出的结果不是完全一样，只是相差几个符号，在matlab中把上述3个矩阵再次相乘（把V求共轭转置），得到的矩阵和A矩阵相同，此时可以说这个函数确实能实现矩阵的奇异值分解。
指令和结果如下：其中V’是求矩阵V的共轭转置。
另一结果：[U,D,V]=QYF(A)U =0.7071 0 0.70710 -1.0000 00.7071 0 -0.7071D =1.4142 00 1.00000 0V =1 00 1矩阵的QR 分解1、原理矩阵QR 分解是一种特殊的三角分解，在解决矩阵特征值的计算、最小二乘法等问题中起到重要作用。
设A ∈C m ×n ，m ≥n 且rankA=n ，则必存在非奇异的上三角n ×n 矩阵R 及m ×n 矩阵Q ，Q H Q=I n ，使得A=QR 。
用Gram-Schmidt 方法对矩阵进行QR 分解时，所论矩阵必须是列满秩矩阵。
不过，不是列满秩的矩阵只要是方阵，也可以作QR 分解，这就是所谓的Householder 方法。
设w ∈C n 是一个单位向量，令 则称H 是一个Householder 矩阵或Householder 变换。
2、算法设u ∈C n 是一个单位向量，则对于任意的C n 存在Householder 矩阵H ，使得Hx=au 。
其中 为实数。
设A 为任一n 阶矩阵，则必存在n 阶酉矩阵Q 和n 阶上三角阵R ，使得A=QR 。
证，第一步，将矩阵A 按列分块写成A=（α1，α2，…，αn ）。
如果α1≠0，则可得，存在n 阶householder 矩阵H 1使得H 1α1=—a 1e 1，| a 1 |=||α1||，e 1∈C n于是有H1A=（H 1α1，H 1α2，…，H 1αn ）= 11*0n a A --⎛⎫⎪⎝⎭如果α1=0，则直接进行下一步，此时相当于取H 1=I n ，而a 1=0.第二步，将矩阵A n-1按列分块写成A n-1=（α1，α2，…，αn-1）。
如果α1≠0，则可得，存在n-1阶householder 矩阵H ’2使得H 1α1=—a 1e 1，| a 1 |=||α1||，e 1∈C n于是有H ’2 A n-1=（H ’2α1，H ’2α2，…，H ’2αn-1）=22*0n a A --⎛⎫⎪⎝⎭此时，令H 2=2100'T H ⎛⎫⎪⎝⎭则H2是n 阶Householder 矩阵，且使H 2H 1A=122**0*0n a a A --⎛⎫⎪- ⎪ ⎪⎝⎭如果α1=0，则直接进行下一步。
第三步，对n-2阶矩阵继续进行类似的变换，如此下去，之多在第n-1步，我们可以找到Householder 矩阵H 1，H 2，…，H n-1使得H n-1…H 2H 1A=12***0**00*00'nn a a α-⎛⎫ ⎪- ⎪⎪⋅⋅⋅ ⎪⎝⎭令Q= H n-1…H 2H 1，则Q 是酉矩阵之积，从而必有酉矩阵并且A=QRH I H ωωω2)(-=u ax x a H ,2=3、程序此实验中用到的有householder子程序和qrfenjie主程序，程序如下：function [R,y]=householder(x,i,j)xi=x(i);xj=x(j);r=sqrt(xi^2+xj^2);cost=xi/r;sint=xj/r;R=eye(length(x));R(i,i)=cost;R(i,j)=sint;R(j,i)=-sint;R(j,j)=cost;y=x(:);y([i,j])=[r,0];function qrfenjie(A)n=size(A,1);R=A;Q=eye(n);for i=1:n-1for j=2:n-i+1x=R(i:n,i);rt=householder(x,1,j);r=blkdiag(eye(i-1),rt);Q=Q*r',R=r*Rendend4、例子以教材上的矩阵A=[0 -3 1；2 1 -6；0 4 2]为例说明上述程序的正确性，在matlab中输入以下指令：5、结果输入上述指令以后，得到如下的结果：此时可以说，上述程序是正确的。