计算方法A上机报告学院（系）：电气工程学院学生姓名：陶然学号：授课老师：完成日期：2019年12月03日西安交通大学Xi'an Jiaotong University目录1 QR分解法求解线性方程组 (2)1.1 算法原理 (2)1.1.1 基于吉文斯变换的QR分解 (2)1.1.2 基于豪斯霍尔德变换的QR分解 (3)1.2 程序流程图 (4)1.2.1 基于吉文斯变换的QR分解流程图 (4)1.2.2 基于豪斯霍尔德变换的QR分解流程图 (5)1.3 程序使用说明 (5)1.3.1 基于吉文斯变换的QR分解程序说明 (5)1.3.2 基于豪斯霍尔德变换的QR分解程序说明 (7)1.4 算例计算结果 (8)2 共轭梯度法求解线性方程组 (10)2.1 算法原理 (10)2.2 程序流程图 (10)2.3 程序使用说明 (11)2.4 算例计算结果 (12)3 三次样条插值 (14)3.1 算法原理 (14)3.2 程序流程图 (16)3.3 程序使用说明 (17)3.4 算例计算结果 (19)4 龙贝格积分 (21)4.1 算法原理 (21)4.2 程序流程图 (22)4.3 程序使用说明 (23)4.4 算例计算结果 (24)结论 (26)1 QR 分解法求解线性方程组1.1 算法原理矩阵的QR 分解是指，可以将矩阵A 分解为一个正交矩阵Q 和一个上三角矩阵R 的乘积，实际中，QR 分解经常被用来解决线性最小二乘问题，分解情况如图1.1所示。

=⨯图1.1 QR 分解示意图本次上机学习主要进行了两个最基本的正交变换—吉文斯（Givens ）变换和豪斯霍尔德（Householder ）变换，并由此导出矩阵的QR 分解以及求解线性方程组的的方法。

1.1.1 基于吉文斯变换的QR 分解吉文斯矩阵也称初等旋转阵，如式（1.1）所示，它把n 阶单位矩阵I 的第,i j 行的对角元改为c ，将第i 行第j 列的元素改为s ，第j 行第i 列的元素改为s −后形成的矩阵。

吉文斯变换矩阵具有以下性质：（1） (),,P i j θ是正交矩阵，()()1T,,,,P i j P i j θθ−=；（2） 仅改变向量的第,i j 个元素； （3） 仅改变矩阵的第,i j 行元素；（4） 可以使用有限个吉文斯矩阵乘积P ，使得12Px x e =。

()111,,111n ncsP i j sc θ⨯⎛⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪=⎪ ⎪ ⎪− ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭i j i j（1.1）矩阵A 按列分块表示为：()12,,,n a a a ，利用性质1,2,4，可以对每一列的各个元素进行计算消去，最终可以得到图1.1的QR 矩阵。

1.1.2 基于豪斯霍尔德变换的QR 分解给定一个n 阶的向量u ，满足21u=，矩阵T 2H I uu =−就是所说的豪斯霍尔德变换矩阵，也称为初等反射阵。

该矩阵存在一定的几何意义，当3n =时，如图1.2所示，对于空间给定的向量x 可以分解为12x x x =+，其中，u 为垂直于π平面的单位法向量，于是有：120x u x ku ==， （1.2）x 1x 2x uHx2x −图1.2 空间向量分解图将空间给定的向量x 和豪斯霍尔德矩阵相乘，有：()()()T 11T T22222Hx I uu x Hx I uu x Hx I uu x ⎧=−⎪=−⇒⎨=−⎪⎩ （1.3）结合式（1.2）可以进行化简，得到：1122Hx x Hx x =⎧⎨=−⎩ （1.4） 得到的向量关于平面对称，因此也称为初等反射矩阵。

豪斯霍尔德矩阵同样存在着优良的性质： （1） 矩阵H 为对称正交矩阵T 1H H H −==；（2） 对于任意给定的单位向量u ，存在初等反射矩阵H ，使得Hx u σ=； 矩阵A 按列分块表示为：()12,,,n a a a ，利用上述性质，可以对每一列进行变换，最终可以得到图1.1的QR 矩阵。

1.2 程序流程图1.2.1 基于吉文斯变换的QR 分解流程图1.2.2 基于豪斯霍尔德变换的QR 分解流程图计算中间参数计算中间参数()()()()2TTsgn (,):,=,,+1:,R i i R i n i R i i R i n i σωσ=−⎡⎤+⎣⎦T 2,2100u H I uu H H ωω==−⎡⎤=⎢⎥⎣⎦图1.3 基于豪斯霍尔德变换的QR 分解流程图由上述流程图也可以明显的看出来，吉文斯变换的QR 分解需要的循环步骤大于豪斯霍尔德的步骤，与理论分析的一致。

1.3 程序使用说明1.3.1 基于吉文斯变换的QR 分解程序说明基于吉文斯变换的QR 分解程序主要有两个函数构成，主函数如下所示，包含相应的注释，主函数主要就是进行QR 求解，以及方程组的求解。

function [Q,R,x]=Givens(A,b)%基于吉文斯变换的QR分解与方程组求解%输入参数说明：A:参数矩阵% b：列向量m=size(A,1); %矩阵的行数n=size(A,2); %矩阵的列数R=A; %给定QR分解矩阵的初始矩阵RQ=eye(m); %给定QR分解矩阵的初始正交矩阵Qfor i=1:n-1 %分别循环对每一列[a1,a2...an] 求解[P1,P2..Pn] for j=i+1:m %对aij元素消0，求解PijAi=R(:,i); %读取第i列的所有元素Pij=Pij_cal(Ai,i,j); %调用Pij矩阵计算函数Q=Q*Pij'; %求解正交矩阵QR=Pij*R; %求解Rendend%根据QR分解结果求解方程组b=Q'*b; %方程组两边同时左乘正交矩阵Q的逆x(n)=b(n)/R(n,n); %计算初值；for i=n:-1:1 %循环求解x的值sum=0; %每次对求和变量进行初始化for j=i+1:nsum=sum+R(i,j)*x(j);endx(i)=(b(i)-sum)/R(i,i);endend子函数的主要是进行每一个元素的消去的矩阵参数计算，具体如下所示：function [Pij]=Pij_cal(Ai,i,j)%计算每个元素的正交矩阵%输入参数说明：Ai:第i行的列向量% i，j：需要消零元素的编号ai=Ai(i); %需要消元行的首元素赋值；aj=Ai(j); %消元的元素赋值；c=sqrt(ai^2+aj^2); %用于计算角度值；cos_theta=ai/c; %计算cos值sin_theta=aj/c; %计算sin值Pij=eye(length(Ai)); %给定m*m的初始单位矩阵；Pij(i,i)=cos_theta;Pij(i,j)=sin_theta;Pij(j,i)=-sin_theta;Pij(j,j)=cos_theta;end1.3.2 基于豪斯霍尔德变换的QR分解程序说明在进行豪斯霍尔德变换的编程时，我们所需奥解决的主要问题是进行QR矩阵的求解，而对于矩阵的相乘以及矩阵基本的运算和求解向量的范数而言，在进行编程的时候并没有使用自己编写的函数，而是调用了MATLAB中相应的算法。

具体的程序及相应的注释如下所示：function [Q,R,x]=householder(A,b)%基于豪斯霍尔德的QR求解与方程组求解m=size(A,1); %矩阵的行数n=size(A,2); %矩阵的列数R=A; %给定QR分解矩阵的初始矩阵RQ=eye(m); %给定QR分解矩阵的初始正交矩阵Qfor i=1:n-1if R(i,i)>=0 %对zgama符号进行判断sigma=-norm(R(i:n,i));elsesigma=norm(R(i:n,i));endw=[R(i,i)-sigma,R(i+1:n,i)']'; %求解中间计算参数wu=w/norm(w,2);H=eye(m-i+1)-2*u*u'; %求解矩阵HiH=blkdiag(eye(i-1),H); %将矩阵扩展为m*m 的矩阵 Q=Q*H'; %求解Q 矩阵 R=H*R; %求解R 矩阵 end%根据QR 分解结果求解方程组b=Q'*b; %方程组两边同时左乘正交矩阵Q 的逆 x(m)=b(m)/R(m,m); %计算初值；for i=n:-1:1 %循环求解x 的值sum=0; %每次对求和变量进行初始化 for j=i+1:nsum=sum+R(i,j)*x(j); endx(i)=(b(i)-sum)/R(i,i); end end1.4 算例计算结果为了验证所编写的程序是否正确，采用了作业题中第74页的计算实习2.4进行验证，该题具体的数据如式（1.5）所示。

54756753941287886537810987556 , 5791197553688910189587877810105756759101052A b ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦（1.5）将相应的数据输入到MATLAB 中，分别运行基于吉文斯变换的QR 分解程序“[Q,R,x]=Givens(A,b)”和豪斯霍尔德程序“[Q,R,x]=householder(A,b)”可以得到相应的方程组的解如式（1.6）所示：[]T1111111x = （1.6）经过比较发现，采用两种算法得到的QR 分解矩阵，如式（1.7）所示，经过比较发现两者得到的结果在正负号方面存在差异，主要是豪斯霍尔德算法中对于σ取值确定上存在正负的情况，因此略有不同。

0.33330.33720.28770.43670.09820.22220.66530.26670.91200.01360.24620.11810.07060.13190.46670.14980.19610.14070.57040.44460.42250.33330.06570.56920.60350.05840.42360.11660.40000.0251givensQ −−−−−−−−−−−=−0.13600.37780.38150.62130.38290.46670.14980.72940.15600.15470.42000.05500.33330.06860.06400.44450.69150.07780.44520.33330.33720.28770.43670.09820.22220.6householderQ ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥−−⎢⎥−−−−⎢⎥⎢⎥−−−⎣⎦−−−=6530.26670.91200.01360.24620.11810.07060.13190.46670.14980.19610.14070.57040.44460.42250.33330.06570.56920.60350.05840.42360.11660.40000.02510.13600.37780.38150.62130.38290.46670.149−−−−−−−−−−−−−−−−−−80.72940.15600.15470.42000.05500.33330.06860.06400.44450.69150.07780.4452⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥−⎢⎥⎢⎥−−−⎣⎦得到的上三角R 矩阵的结果如下所示15.000019.533320.933319.933321.600021.266719.80000.00007.4464 2.6996 2.90553.0995 2.3624 1.10660.00000.0000 3.2416 3.3580 1.68830.4811 2.30360.00000.00000.0000 3.73420.7405 1.6506 1.11390.00000.givensR −−−−−=−−−−householder00000.00000.0000 3.2303 3.2048 3.90190.00000.00000.00000.00000.0000 1.98010.07380.00000.00000.00000.00000.00000.00000.978615.000019.533320.933319.93321.60R ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥−−⎢⎥−−−−⎢⎥⎢⎥−−−−⎣⎦−−−−−=0021.266719.80000.00007.4464 2.6996 2.9055 3.0995 2.3624 1.10660.00000.0000 3.2416 3.3580 1.68830.4811 2.30360.00000.00000.0000 3.73420.7405 1.6506 1.11390.00000.00000.00000.0000 3.2303 3.2048 3.−−−−−−−−−−−−−−−−−−−−90190.00000.00000.00000.00000.0000 1.98010.07380.00000.00000.00000.00000.00000.00000.9786⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥−−−⎢⎥⎢⎥⎣⎦为了验证所计算结果是否正确，采用MATLAB 中关于QR 分解的库算法，在命令栏中输入“[Q,R]=qr(A)”得到的结果与采用豪斯霍尔德的计算结果相同，在命令栏中输入“[x]=inv(A)*b ”得到的方程组的解与计算结果相同，验证了算法的正确性。