数值分析第一次上机实验练习报告.docx

数值分析第一次上机实验练习报告 ——线性方程组求解相关问题 材博141 李根 2014310937 问题的描述 设Hn=[hij]∈Rn×n 是Hilbert矩阵，即 hij=1i+j-1. 对n=2,3,4…（建议算到15左右） 取x=1?1∈Rn，及bn=Hnx. 再用Gauss消去法和Cholesky分解方法来求解Hny=bn，看看误差有多大. 计算条件数：cond(Hn)2. 使用某种正则化方法改善(a)中的结果. 方法描述与方案设计 我们使用MATLAB来完成此次上机实验。MATLAB是一款强大的科学计算软件，内置有多种函数可以方便地调用，同时软件界面非常友好，便于学习和使用。我们的实验方案设计如下： Hilbert矩阵的生成与条件数的计算 MATLAB内置有Hilbert矩阵的生成函数hilb()和矩阵条件数的计算函数cond()，直接调用即可。 Gauss消去法的MATLAB实现 对线性方程组Hy=b，将H和b写为增广矩阵H|b，令之为H(1)|b(1)。顺序消元是对增广矩阵进行初等行变换，经n-1步消元将其化为H(n)|b(n)，其中H(n)位一个上三角矩阵。过程如下： 经第k-1步消元（1≤k≤n），对应增广矩阵化为 H(k)|b(k)=h11(1)h12(1)… h22(2)… ?h1k(1)…h1n(1)h2k(2)…h2n(2)? ? hkk(k)…hkn(k)? ?hnk(k)…hnn(k)b1(1)b2(2)?bk(k)?bn(k) 假设akk(1)≠0，令 lik=hik(1)hkk(1)，i=k+1, k+2, …, n. 前k-1行不动，将第k行乘以-lik加到第i行，进行第k步消元，即 hij(k+1)=hij(k)-likhkjk，i,j=k+1, k+2, …, n.bi(k+1)=bi(k)-likbk(k)，i=k+1, k+2, …, n. 这样经n-1步后，得到了矩阵如下H(n)|b(n)，其中H(n)为一个上三角矩阵，当方程有唯一解时，H可逆，而经初等变换后的系数矩阵行列式不变，即 detH(n)=detH≠0 H(n)|b(n)=h11(1)h12(1) h22(2)…h1n(1)…h2n(2) ?? hnn(n)b1(n)b2(n)?bn(n) 所以hnn(1)≠0。解方程组H(n)y=b(n)，从最后一个方程开始，依次向前一方程代入，即可求解，此过程为回代，计算公式为 yn=bn(n)/hnn(n)，yi=1hiinbin-j=i+1nhijnxj，i=n-1, …, 2, 1. 通过查询MATLAB应用的相关书籍，我们发现，MATLAB解线性方程组常用的方法有两种，一是在方程Ax=b两边同时乘以A的逆矩阵A-1，二是利用预先定义好的矩阵“左除”：x=A\b，而后者的算法即为Gauss消元法，因此对于某个给定的Hilbert矩阵Hn，我们进行矩阵左除运算所得的结果即为Gauss消去法的结果。 Cholesky分??法的MATLAB实现 Cholesky分解法又称平方根法，相较于Gauss消去法，Cholesky分解法求解线性方程组时具有较小的计算量，但要求系数矩阵对称且正定。对系数矩阵H作Cholesky分解即寻找下三角矩阵L，使得 H=LLT L=l11 l21l22 ??ln1ln2? …lnn hij=k=1j-1likljk+lijljj，i=j, j+1, …, n. 这样的L对角元全为正数且唯一。分解时，由上式首先可计算出第一列的元素，当我们计算出前j-1列的元素时，第j列的对角元素由下式算出 ljj=hjj-k=1j-1ljk212 从而计算出第j列元素 lij=1ljjhjj-k=1j-1likljk ，i=j+1, …, n. 然后我们可依次求解方程组Lx=b、LTy=x，最终求解，即 x1=b1/l11，xi=1liibi-k=1i-1likxk，i=2, 3, …,n. yn=xn/lnn，xi=1liixi-k=i+1nlkixk，i=n-1, n-2, …,1. MATLAB内置有矩阵Cholesky分解的函数chol()，但需要注意，函数给出的矩阵Rn为上三角阵，和课本上的定义略有不同。得到Rn后，再利用Gauss消去法求解两个新的线性方程组即可。 Тихонов正则化方法的MATLAB实现 所谓Тихонов正则化方法其实是一种近似方法，其核心是通过引入了一个小的纯量矩阵来改善原病态矩阵的性质，使其条件数显著降低从而得到更精确的解。 对本实验而言，我们令 An=αIn+HnTH