数值线性代数实验报告

迭代法解线性方程组数值分析实验报告一、实验目的本次实验旨在深入研究和掌握迭代法求解线性方程组的基本原理和方法，并通过数值实验分析其性能和特点。

具体目标包括：1、理解迭代法的基本思想和迭代公式的推导过程。

2、掌握雅克比（Jacobi）迭代法、高斯赛德尔（GaussSeidel）迭代法和超松弛（SOR）迭代法的算法实现。

3、通过实验比较不同迭代法在求解不同类型线性方程组时的收敛速度和精度。

4、分析迭代法的收敛性条件和影响收敛速度的因素。

二、实验原理1、线性方程组的一般形式对于线性方程组＄Ax ＝ b$，其中＄A$ 是＄n×n$ 的系数矩阵，＄x$ 是＄n$ 维未知向量，＄b$ 是＄n$ 维常向量。

2、迭代法的基本思想迭代法是从一个初始向量＄x^｛（0)｝＄出发，按照某种迭代公式逐步生成近似解序列＄＼｛x^｛（k)｝＼｝＄，当迭代次数＄k$ 足够大时，＄x^｛（k)｝＄逼近方程组的精确解。

3、雅克比迭代法将系数矩阵＄A$ 分解为＄A ＝ D L U$，其中＄D$ 是对角矩阵，＄L$ 和＄U$ 分别是下三角矩阵和上三角矩阵。

雅克比迭代公式为：＄x^｛（k+1)｝＝ D^｛－1}（b ＋（L ＋ U)x^｛（k)｝）＄。

4、高斯赛德尔迭代法在雅克比迭代法的基础上，每次计算新的分量时立即使用刚得到的最新值，迭代公式为：＄x_i^｛（k+1)｝＝（b_i ＼sum_{j=1}＾｛i-1}a_{ij}x_j^｛（k+1)｝＼sum_{j=i+1}＾｛n}a_{ij}x_j^｛（k)｝）／a_{ii}＄。

5、超松弛迭代法在高斯赛德尔迭代法的基础上引入松弛因子＄＼omega$，迭代公式为：＄x_i^｛（k+1)｝＝ x_i^｛（k)｝＋＼omega(（b_i ＼sum_{j=1}＾｛i-1}a_{ij}x_j^｛（k+1)｝＼sum_{j=i}＾｛n}a_{ij}x_j^｛（k)｝）／ a_{ii} x_i^｛（k)｝）＄。

数值代数实验数值线性代数实验一一、实验名称：矩阵的LU分解.二、实验目的：用不选主元的LU分解和列主元LU分解求解线性方程组Ax=b, 并比较这两种方法.三、实验内容与要求（1）用所熟悉的计算机语言将不选主元和列主元LU分解编成通用的子程序，然后用编写的程序求解下面的84阶方程组将计算结果与方程组的精确解进行比较，并就此谈谈你对Gauss消去法的看法.（2）写出追赶法求解三对角方程组的过程，并编写程序求该实验中的方程组Gauss消去法：用消去法解方程组的基本思想是用逐次消去未知数的方法把原来方程组Ax=b化为与其等价的三角方程组，而求解三角方程组就容易了。

换句话说，上述过程就是用行的初等变换将原方程组系数矩阵化为简单形式，从而将求解原方程组的问题转化为求解简单方程组的问题。

利用Gauss消去法对线性方程组Ax=b进行求解。

用MATLAB建立m文件DelGauss.m,程序如下：function x=DelGauss(a,b)[n,m]=size(a);nb=length(b);det=1;x=zeros(n,1);for k=1:n-1for i=k+1:nif a(k,k)==0returnendm=a(i,k)/a(k,k);for j=k+1:na(i,j)=a(i,j)-m*a(k,j);endb(i)=b(i)-m*b(k);enddet=det*a(k,k);enddet=det*a(n,n);for k=n:-1:1for j=k+1:nb(k)=b(k)-a(k,j)*x(j);endx(k)=b(k)/a(k,k);End在matlab中输入如下：结果如下：方程组的精确解为x 1=x 2=…=x 84=1.0000，与Gauss 消去法求得的解差距很大，所得结果不够准确，计算简单但其消元过程有时不能进行到底而使求解出现解失真的情况。

数值线性代数实验二一、实验名称：实对称正定矩阵的A的Cholesky分解.二、实验目的：用平方根法和改进的平方根方法求解线性方程组Ax=b.三、实验内容与要求用所熟悉的计算机语言将Cholesky分解和改进的Cholesky分解编成通用的子程序，然后用编写的程序求解对称正定方程组Ax=b，其中（1）b随机的选取，系数矩阵为100阶矩阵（2）系数矩阵为40阶Hilbert矩阵，即系数矩阵A的第i行第j列元素为，向量b的第i个分量为（3）用实验一的程序求解这两个方程组，并比较所有的计算结果，然后评价各个方法的优劣。

一、引言很多情况下在工程中抽象出来的数学方程组是超定的，没有精确解，这样就需要找一个在某种意义下最接近精确解的解。

设A 是m*n 的实矩阵，22A (A )T x bQ x b -=-,这样2min A y b -就等价与2min (A )T Q x b -,为方便求解，需要A T Q 是上三角矩阵，这样引入QR 分解就比较求解方便。

二、豪斯霍尔德变换豪斯霍尔德变换（Householder transformation ）又称初等反射（Elementary reflection ），最初由A.C Aitken 在1932年提出，Alston Scott Householder 在1958年指出了这一变换在数值线性代数上的意义。

这一变换将一个向量变换为由一个超平面反射的镜像，是一种线性变换。

其变换矩阵被称作豪斯霍尔德矩阵，在一般内积空间中的类比被称作豪斯霍尔德算子，超平面的法向量被称作豪斯霍尔德向量。

三、理论依据householder 变换：任意的一个向量x ，经过一个正交变换2**T H I W W =- 后总可以变为一个与之范数相等的另一个向量Hx 。

如上图中所示，记v Hx x =-，2/w v v =，2**TH I W W =-上述H 既为所要求的householder 变换。

具体操作时将需要变化的矩阵的每一列当做一个向量，第一列变为除了第一个元素都为0的向量…第i 列变为除了前i 个元素都为0的向量...为了保证在每次变化时不改变已经变好的0元素，第i 次只变化每列第i 个元素到第n 个元素，每个变化矩阵的形式是[I 0;0 H]。

算法如下：qr_Houserhoilder.m function [Q,R]=qr_Houserhoilder(A) [m,n]=size(A); Q=eye(n); for i=1:nu=house(A(i:m,i));P=eye(m-i+1)-2/(norm(u)^2)*u*u'; A(i:m,i:n)=P*A(i:m,i:n);PP=blkdiag(eye(i-1),P);Q=Q*PP;endR=triu(A);程序house.mfunction u=house(Ai)Ai(1)=Ai(1)+sign(Ai(1))*norm(Ai);u=Ai/norm(Ai);四、数值实验结果：为了保证在一定阶数下，所做的数值试验具有代表性，我们随机产生3个n阶矩阵，随着n的不断变化，得到如下一系列图：（1）A=rand(4,4)A =0.3027 0.7894 0.1793 0.39220.2364 0.7134 0.4350 0.39360.4690 0.0742 0.9747 0.11150.5381 0.6768 0.7245 0.8650[Q,R]=qr_Houserhoilder(A)Q =-0.3735 -0.5369 0.4602 0.6003-0.2916 -0.5439 -0.7840 -0.0668-0.5786 0.6445 -0.2679 0.4219-0.6639 -0.0208 0.3190 -0.6761R =-0.8106 -0.9951 -1.2387 -0.90010 -0.7781 0.2803 -0.37070 0 -0.2886 0.11800 0 0 -0.3286test_qr3err =1.0e-10 *0.1950 0.0120 0.0119 0.3029-111 1.52 2.53 3.54（2）A=rand(5,5)A =0.7575 0.7005 0.6420 0.8233 0.88770.6373 0.4698 0.5815 0.5290 0.37700.6643 0.3036 0.9458 0.7210 0.12950.1056 0.0254 0.3903 0.1048 0.30270.2813 0.3787 0.5302 0.0627 0.1156[Q,R]=qr_Houserhoilder(A)Q =-0.6161 0.4378 0.2968 -0.5805 -0.0598-0.5183 -0.0381 0.2932 0.7098 -0.3743-0.5403 -0.6729 -0.2440 -0.1030 0.4303-0.0859 -0.1847 -0.5513 -0.2507 -0.7693-0.2288 0.5657 -0.6801 0.2928 0.2817R =-1.2295 -0.9280 -1.3629 -1.1944 -0.86480 0.2940 -0.1496 -0.1288 0.29660 0 -0.4455 0.1231 0.09690 0 0 -0.1847 -0.30310 0 0 0 -0.3388test_qr3err =1.0e-10 *0.1722 0.0119 0.0120 0.4239-111 1.52 2.53 3.54（3）A=rand(6,6)A =0.5961 0.8422 0.4279 0.3473 0.2731 0.34860.1335 0.2860 0.2157 0.6181 0.6306 0.19300.4479 0.0467 0.3284 0.4588 0.4228 0.95480.9430 0.7296 0.7951 0.0920 0.0322 0.31140.7581 0.7811 0.4338 0.1512 0.6185 0.12710.8647 0.1359 0.1735 0.9036 0.4393 0.8964 [Q,R]=qr_Houserhoilder(A)Q =-0.3571 -0.5571 0.2402 0.2038 0.4576 0.5035-0.0800 -0.2530 -0.1880 0.8253 -0.1645 -0.4314-0.2684 0.3532 -0.4622 0.2419 -0.3487 0.6399-0.5650 -0.0874 -0.6012 -0.3372 0.2898 -0.3376-0.4542 -0.3266 0.3151 -0.2560 -0.7185 -0.0773-0.5181 0.6217 0.4823 0.1991 0.2045 -0.1760R =-1.6690 -1.1737 -0.9944 -0.8854 -0.7882 -1.09430 -0.7594 -0.2803 0.3166 -0.0940 0.58290 0 -0.3472 0.1832 0.1390 -0.10880 0 0 0.8020 0.5966 0.50220 0 0 0 -0.4714 -0.02290 0 0 0 0 0.4305test_qr3err =1.0e-10 *0.1518 0.0119 0.0120 0.2845-111 1.52 2.53 3.54五、总结与分析由上面的数值试验结果和绘图可知：比较之下，我自己写的houserhoilder程序误差是最大的，CGS与MGS的误差相差不大，系统自带的matlab软件速度是最快的，我自己写的houserhoilder程序还是最慢的。

数值代数实验报告Numerical Linear Algebra And ItsApplications学生所在学院：理学院学生所在班级：计算数学10-1学生姓名：戈东潮指导教师：于春肖教务处2012年12月实验一实验名称： Poisson 方程边值问题的五点差分格式 实验时间: 2012年12月13日 星期四 实验成绩： 一、实验目的：通过上机利用Matlab 数学软件实现Poisson 方程的边值问题的五点差分格式的线性代数方程组来认真解读Poisson 方程边值问题的具体思想与方法，使我们掌握得更加深刻，并且做到学习与使用并存，增加学习的实际动手性，不再让学习局限于书本和纸上，而是利用计算机学习来增加我们的学习兴趣。

二、实验内容：利用Poisson 方程来解下列问题：⎪⎩⎪⎨⎧∂∈=∈=∂∂+∂∂-ΩΩ),(,),(),(,sin sin )(y x y x u y x y x y u x u 0222222πππ 其中}{的边界是ΩΩ∂<<∈Ω,1,0),(y x y x 。

边值问题的解释是y x y x u ππsin sin ),(=（1）用例题中模型问题的方法取N=5,10（也可以取N=20）列出五点差分格式的线性代数方程组三、实验过程：系数矩阵A 的Matlab 实现函数如下：%文件名：poisson.m function T=poisson(N) for i=1:N a(i)=4; end b=diag(a); for j=1:N-1b(j,j+1)=-1;endfor j=2:Nb(j,j-1)=-1;endfor i=1:NI(i)=-1;endI=diag(I);for i=1:N:N*Nfor j=1:NT(i+j-1,i:i+N-1)=b(j,:);endendfor i=1:N:N*N-Nfor j=1:NT(i+j-1,i+N:i+N-1+N)=I(j,:);endendfor i=1+N:N:N*Nfor j=1:NT(i+j-1,i-N:i-1)=I(j,:);endend求解常数项b的Matlab实现函数如下：%函数名：constant.mfunction b=constant(N)n=N+1;h=1/n;i=1;for y=1/n:1/n:N/nfor x=1/n:1/n:N/nf(i)=2*pi*pi*sin(pi*x)*sin(pi*y);i=i+1;endendb=h*h*f;b=b';四、实验结果(总结/方案)在Matlab运行窗口输入>>A=poisson(5) Enter输出结果A =Columns 1 through 114 -1 0 0 0 -1 0 0 0 0 0-1 4 -1 0 0 0 -1 0 0 0 00 -1 4 -1 0 0 0 -1 0 0 00 0 -1 4 -1 0 0 0 -1 0 00 0 0 -1 4 0 0 0 0 -1 0-1 0 0 0 0 4 -1 0 0 0 -10 -1 0 0 0 -1 4 -1 0 0 00 0 -1 0 0 0 -1 4 -1 0 00 0 0 -1 0 0 0 -1 4 -1 00 0 0 0 -1 0 0 0 -1 4 00 0 0 0 0 -1 0 0 0 0 40 0 0 0 0 0 -1 0 0 0 -10 0 0 0 0 0 0 -1 0 0 00 0 0 0 0 0 0 0 -1 0 00 0 0 0 0 0 0 0 0 -1 00 0 0 0 0 0 0 0 0 0 -10 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 Columns 12 through 220 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 00 -1 0 0 0 0 0 0 0 0 00 0 -1 0 0 0 0 0 0 0 00 0 0 -1 0 0 0 0 0 0 0 -1 0 0 0 -1 0 0 0 0 0 04 -1 0 0 0 -1 0 0 0 0 0 -1 4 -1 0 0 0 -1 0 0 0 00 -1 4 -1 0 0 0 -1 0 0 00 0 -1 4 0 0 0 0 -1 0 00 0 0 0 4 -1 0 0 0 -1 0 -1 0 0 0 -1 4 -1 0 0 0 -10 -1 0 0 0 -1 4 -1 0 0 00 0 -1 0 0 0 -1 4 -1 0 00 0 0 -1 0 0 0 -1 4 0 00 0 0 0 -1 0 0 0 0 4 -10 0 0 0 0 -1 0 0 0 -1 40 0 0 0 0 0 -1 0 0 0 -10 0 0 0 0 0 0 -1 0 0 00 0 0 0 0 0 0 0 -1 0 0 Columns 23 through 250 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 -1 0 0 0 -1 0 0 0 -1 0 0 4 -1 0 -1 4 -10 -1 4 在运行窗口输入 >>b=constant(5) Enter 常数项的输出结果b =[ 0.1371 0.2374 0.2742 0.2374 0.1371 0.2374 0.4112 0.4749 0.4112 0.2374 0.2742 0.4749 0.5483 0.4749 0.2742 0.2374 0.4112 0.4749 0.4112 0.2374 0.1371 0.2374 0.2742 0.2374 0.1371]T所以A*u=b 的五点差分格式为j i j i j i j i j i j i f h u u u u u ,,,,,,211114=-----+-+ （i,j=1,2……）实验二实验名称： 用Jacobi 迭代法和SOR 迭代法求解方程组 实验时间: 2012年12月13日 星期四 实验成绩： 一、实验目的：通过上机利用Matlab 数学软件采用Jacobi 迭代法和SOR 迭代法来求解方程组，求解过程中了解两种迭代方法的基本思想与迭代过程，并分析两种迭代方法的收敛性和实用用以及他们的异同点。

数学实验报告题目第一次实验题目一、 实验目的1．熟悉MATLAB 的矩阵初等运算；2．掌握求矩阵的秩、逆、化最简阶梯形的命令； 3．会用MABLAB 求解线性方程组二、 问题求解和程序设计流程1． 已知⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=351503224A ，⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---=112302431B ，在MATLAB 命令窗口中建立A 、B矩阵并对其举行以下操作：(1) 计算矩阵A 的行列式的值det()?A = (2) 分别计算下列各式：B A -2 、 B A *和B A *.、 1-AB 、 B A 1-、 2A 、 T A解：（1） 编写程序如下：A=[4 -2 2;-3 0 5;1 5 3];B=[1 3 4;-2 0 -3;2 -1 1]; a=det(A) 运行结果： a = -158（2）编写程序如下： C=2*A-BD=A*B E=A.*B F=A/BG=A\B H=A*A K=A'运行结果：C =7 -7 0 -4 0 13知识归纳整理求知若饥，虚心若愚。

线性代数实验报告0 11 5D =12 10 247 -14 -7-3 0 -8E =4 -6 86 0 -152 -5 3F =0 0 2.0000-2.7143 -8.0000 -8.14292.42863.0000 2.2857G =0.4873 0.4114 1.00000.3671 -0.4304 0-0.1076 0.2468 0H =24 2 4-7 31 9-8 13 36K =4 -3 1-2 0 52 5 32．在MATLAB中分别利用矩阵的初等变换及函数rank、函数inv求下列矩阵的秩：线性代数实验报告(1) ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----=4211104532361A 求 Rank(A)=? (2) 3501120010201202B ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦求?1=-B 解：（1）编写程如下：format ratA=[1 -6 3 2;3 -5 4 0;-1 -11 2 4]; rref(A) 运行结果： ans =1 0 0 -8/5 0 1 0 0 0 0 1 6/5 由A 经初等变换后得到的行最简型可知：A 的秩为3。

线性代数实验报告一、实验目的线性代数是一门重要的数学基础课程，它在工程、科学、计算机等领域都有着广泛的应用。

本次实验的目的是通过实际操作和计算，加深对线性代数基本概念和方法的理解，提高运用线性代数知识解决实际问题的能力。

二、实验环境本次实验使用了软件名称软件进行计算和绘图。

三、实验内容（一）矩阵的运算1、矩阵的加法和减法给定两个矩阵 A 和 B，计算它们的和 A ＋ B 以及差 A B。

观察运算结果，验证矩阵加法和减法的规则。

2、矩阵的乘法给定两个矩阵 C 和 D，其中 C 的列数等于 D 的行数，计算它们的乘积 CD。

分析乘法运算的结果，理解矩阵乘法的意义和性质。

（二）行列式的计算1、二阶和三阶行列式的计算手动计算二阶和三阶行列式的值，熟悉行列式的展开法则。

使用软件验证计算结果的正确性。

2、高阶行列式的计算选取一个四阶或更高阶的行列式，利用软件计算其值。

观察行列式的值与矩阵元素之间的关系。

（三）线性方程组的求解1、用高斯消元法求解线性方程组给定一个线性方程组，将其增广矩阵通过初等行变换化为行阶梯形矩阵。

求解方程组的解，并验证解的正确性。

2、用矩阵的逆求解线性方程组对于系数矩阵可逆的线性方程组，计算系数矩阵的逆矩阵。

通过逆矩阵求解方程组，并与高斯消元法的结果进行比较。

（四）向量组的线性相关性1、判断向量组的线性相关性给定一组向量，计算它们的线性组合是否为零向量。

根据计算结果判断向量组的线性相关性。

2、求向量组的极大线性无关组对于给定的向量组，通过初等行变换找出极大线性无关组。

（五）特征值和特征向量的计算1、计算矩阵的特征值和特征向量给定一个矩阵，计算其特征值和对应的特征向量。

验证特征值和特征向量的定义和性质。

2、利用特征值和特征向量进行矩阵对角化对于可对角化的矩阵，将其化为对角矩阵。

四、实验步骤（一）矩阵的运算1、首先在软件中输入矩阵 A 和 B 的元素值。

2、然后使用软件提供的矩阵加法和减法功能，计算 A ＋ B 和 A B 的结果。

线性代数方程组的直接解法实验目的：线性方程组求解的直接法编程实现，实验内容：线性方程组求解的高斯消去法算法实现线性方程组求解的主元素消去法算法实现线性方程组求解的LU 分解得方法算法实现线性方程组求解追赶法算法实现实验比较：高斯消去、主元素消去、LU 分解都用实例 ⎪⎩⎪⎨⎧-=-+-=+-=++15.3181533126321321321x x x x x x x x x 这个进行比较。

知识理论解线性方程组的方法大致分为直接法和迭代法。

直接法是指假设计算过程中不产生舍入误差，经过有限次的运算可求的方程组精确解的方法。

方程（2-1）⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++nn nn n n n n n n b x a x a x a b x a x a x a b x a x a x a ...... (22112)222212111212111 记为：AX=b ；一、高斯（Gauss ）消元法(1).Gauss 消元法是最基本的一种方法。

先逐次消去变量，将方程组化成同解的上三角方程组。

消元过程：先逐次消去变量，将方程组化成同解的上三角方程组 基本思想回代过程：按方程的相反顺序求解三角形方程组，得到原方程组的解。

(2) Gauss 消去法的求解思路为：若先让第一个方程组保持不变，利用它消去其余方程组中的,使之变成一个关于变元的n-1阶方程组。

按照（1）中的思路继续运算得到更为低阶的方程组。

经过n-1步的消元后，得到一个三角方程。

利用求解公式回代得到线性方程组的解。

①消元过程：第一次消元，,0)1(11≠a设 )1(11)1(11a a l i i =记（i=1,2,3……，n ）.将（2-1）中第i 个方程减去第一个方程乘以1i l （i=1,2,3……，n ），完成第一次消元，⎪⎩⎪⎨⎧=++=++=+++)2()2(2)2(2)2(2)2(22)2(22)1(1)1(12)1(121)1(11...... ... n n nn n n n n n b x a x a b x a x a b x a x a x a 得等价方程组 （2-2）)2()2(:b x A =简记为其中：)1(11)1()2(j i ij ij a l a a -=；)1(11)1()2(b l b b i i i -=次消元后经过1-k ：⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎨⎧=+++=+++=+++=+++++=++++++++++++++++++++... ...... ......... ......)()(1)(1,)()(1)(,11)(1,1)(,1)()(1)(1,)()2(2)2(21)2(1,2)2(22)2(22)1(1)1(11)1(1,1)1(12)1(121)1(11k n n k nn k k k n k k nk k k n k n k k k k k k k k k k kn k kn k k k k k k kk n n k k k k n n k k k k b x a x a x a b x a x a x a b x a x a x a b x a x a x a x a b x a x a x a x a x a 简记为：)1()1(++=k k b x A其中：)()()1(k kjik k ij k ija l a a -=+ )()()1(k k ik k i k ib l b b -=+),,1,(n k j i +=按上述方法完成n-1次消元后得到。

（1）估计5到20阶Hilbert 矩阵的∞范数条件数（2）设n n R A ⨯∈⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡------=111111111011001ΛΛO O MM M O OΛ,先随机地选取n R x ∈,并计算出x A b n =;然后再用列主元Gauss 消去法求解该方程组,假定计算解为∧x 。

试对n 从5到30估计计算解∧x 的精度,并且与真实相对误差作比较。

解(1)分析:利用for 使n 从5循环到20,利用()hilb 函数得到Hilbert 矩阵A ;先将算法2、5、1编制成通用的子程序,利用算法2、5、1编成的子程序)(B opt v =,对TAB -=求解,得到∞-1A的一个估计值v v =~;再利用inf),(A norm 得到∞A ;则条件数inf),(1A norm v A A K *==∞∞-。

另,矩阵A 的∞范数条件数可由inf),(A cond 直接算出,两者可进行比较。

程序为1 算法2、5、1编成的子程序)(B opt v =function v=opt(B)k=1;n=length(B); x=1、/n*ones(n,1);while k==1 w=B*x;v=sign(w); z=B'*v;if norm(z,inf)<=z'*x v=norm(w,1); k=0; elsex=zeros(n,1);[s,t]=max(abs(z)); x(t)=1; k=1; end end end2 问题(1)求解 ex2_1for n=5:20A=hilb(n);B=inv(A、');v=opt(B);K1=v*norm(A,inf);K2=cond(A,inf);disp(['n=',num2str(n)])disp(['估计条件数为',num2str(K1)])disp(['实际条件数为',num2str(K2)])end计算结果为n=5估计条件数为943656实际条件数为943656n=6估计条件数为29070279、0028实际条件数为29070279、0028n=7估计条件数为985194887、5079实际条件数为985194887、5079n=8估计条件数为33872789099、7717实际条件数为33872789099、7717n=9估计条件数为16、422实际条件数为16、422n=10估计条件数为35353368771750、67实际条件数为35353368771750、67n=11估计条件数为1232433965549344实际条件数为1232433965549344Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 2、547634e-17、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 2、547634e-17、> In cond at 47In ex2_1 at 6n=12估计条件数为3、9245e+16实际条件数为3、9245e+16Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 7、847381e-19、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 7、847381e-19、> In cond at 47In ex2_1 at 6n=13估计条件数为1、2727e+18实际条件数为1、2727e+18Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 2、246123e-18、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 2、246123e-18、> In cond at 47In ex2_1 at 6n=14估计条件数为4、8374e+17实际条件数为4、8374e+17Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 8、491876e-19、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 8、491876e-19、> In cond at 47In ex2_1 at 6n=15估计条件数为4、6331e+17实际条件数为5、234289848563619e+17Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 9、137489e-19、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 9、137489e-19、> In cond at 47In ex2_1 at 6n=16估计条件数为8、3166e+17实际条件数为8、3167e+17Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 6、244518e-19、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 6、244518e-19、 > In cond at 47 In ex2_1 at 6 n=17估计条件数为1、43e+18 实际条件数为1、43e+18Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 4、693737e-19、 > In ex2_1 at 3Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 4、693737e-19、 > In cond at 47 In ex2_1 at 6 n=18估计条件数为2、5551e+18 实际条件数为2、8893e+18Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 4、264685e-19、 > In ex2_1 at 3Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 4、264685e-19、 > In cond at 47 In ex2_1 at 6 n=19估计条件数为2、411858563109357e+18 实际条件数为2、411858563109357e+18Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 1、351364e-19、 > In ex2_1 at 3Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 1、351364e-19、 > In cond at 47 In ex2_1 at 6 n=20估计条件数为2、31633670586674e+18 实际条件数为6、37335273308473e+18结果分析随着矩阵阶数增加,估计值误差开始出现,20,17,16,15 n 时估计条件数与实际值存在误差;且条件数很大,Hilbert 矩阵为病态的。

课程实验报告COURSE PAPER课程名称：数值代数与计算方法课程号：08305114授课教师：学号：姓名：所属：计算机科学与工程打印时间：2015评语：题目一：算法：1、对于问题一：2、对于问题二：直接编写递归函数程序，算出三个差分方程的10个近似值程序:1）%main.mclc;clear all;a=1;b=-2;c=-3;[x1,x2]=roots(a,b,c)%roots.mfunction [x1,x2]=roots(a,b,c)d=sqrt(b*b-4*a*c);if d>0x1=(-2*c)/(b+d);x2=(-b-d)/(2*a);elseif d<0x1=(-b+d)/(2*a);x2=(-2*c)/(b-d);end2）%solu.mfunction [X,R,P,Q]=solu(X0,R0,P0,P1,Q0,Q1) X(1)=X0;R(1)=R0;P(1)=P0;P(2)=P1;Q(1)=Q0;Q(2)=Q1;for i=1:9X(i+1)=X(i)/2;X(i)=X(i+1);endfor i=1:9R(i+1)=R(i)/2;R(i)=R(i+1);endfor i=3:10P(i)=3/2*P(i-1)-1/2*P(i-2);P(i-1)=P(i);P(i-2)=P(i-1);endfor i=3:10Q(i)=5/2*Q(i-1)-Q(i-2);Q(i-1)=Q(i);Q(i-2)=Q(i-1);end%Xn.mclc;clear all;X0=1;R0=0.994;P0=1;P1=0.497;Q0=1;Q1=0.497;[X,R,P,Q]=solu(X0,R0,P0,P1,Q0,Q1);for i=1:10x(i)=i;endplot(x,X-R,'r*');hold on结果：1）x1 =3x2 =-12）题目二：1、此题我们以第一个式子为例，求出g(Xn+1)= (2/(2+3*Xn-Xn.^3))^(1/2)2、利用二分法：3、利用牛顿迭代法，并根据题目提示的立方根算法：此题我们以第三个式子为例：，设A=71/3 ,求出x的3次方便是近似值。

1 / 8数值分析实验六：解线性方程组的迭代法2016113 张威震1 病态线性方程组的求解1.1 问题描述理论的分析表明，求解病态的线性方程组是困难的。

实际情况是否如此，会出现怎样的现象呢？实验内容：考虑方程组Hx=b 的求解，其中系数矩阵H 为Hilbert 矩阵，,,1(),,,1,2,,1i j n n i j H h h i j n i j ⨯===+-这是一个著名的病态问题。

通过首先给定解（例如取为各个分量均为1）再计算出右端b 的办法给出确定的问题。

实验要求：（1）选择问题的维数为6，分别用Gauss 消去法、列主元Gauss 消去法、J 迭代法、GS 迭代法和SOR 迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何？（2）逐步增大问题的维数（至少到100），仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？（3）讨论病态问题求解的算法1.2 算法设计首先编写各种求解方法的函数，Gauss 消去法和列主元高斯消去法使用实验5中编写的函数myGauss.m 即可，Jacobi 迭代法函数文件为myJacobi.m ，GS 迭代法函数文件为myGS.m ，SOR 方法的函数文件为mySOR.m 。

1.3 实验结果1.3.1 不同迭代法球求解方程组的结果比较选择H 为6*6方阵，方程组的精确解为x* = (1, 1, 1, 1, 1, 1)T ，然后用矩阵乘法计算得到b ，再使用Gauss 顺序消去法、Gauss 列主元消去法、Jacobi 迭代法、G-S 迭代法和SOR 方法分别计算得到数值解x1、x2、x3、x4，并计算出各数值解与精确解之间的无穷范数。

Matlab 脚本文件为Experiment6_1.m 。

迭代法的初始解x 0 = (0, 0, 0, 0, 0, 0)T ，收敛准则为||x(k+1)-x(k)||∞<eps=1e-6，SOR方法的松弛因子选择为w=1.3，计算结果如表1。

（1）估计5到20阶Hilbert 矩阵的∞数条件数（2）设n n R A ⨯∈⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡------=111111111011001，先随机地选取n R x ∈，并计算出x A b n =；然后再用列主元Gauss 消去法求解该方程组，假定计算解为∧x 。

试对n 从5到30估计计算解∧x 的精度，并且与真实相对误差作比较。

解（1）分析：利用for 使n 从5循环到20，利用()hilb 函数得到Hilbert 矩阵A ；先将算法2.5.1编制成通用的子程序，利用算法2.5.1编成的子程序)(B opt v =，对TA B -=求解，得到∞-1A的一个估计值v v =~；再利用inf),(A norm 得到∞A ；则条件数inf),(1A norm v A A K *==∞∞-。

另，矩阵A 的∞数条件数可由inf),(A cond 直接算出，两者可进行比较。

程序为1 算法2.5.1编成的子程序)(B opt v =function v=opt(B)k=1;n=length(B);x=1./n*ones(n,1); while k==1 w=B*x;v=sign(w); z=B'*v;if norm(z,inf)<=z'*x v=norm(w,1); k=0; elsex=zeros(n,1);[s,t]=max(abs(z)); x(t)=1; k=1; end end end2 问题（1）求解 ex2_1for n=5:20A=hilb(n);B=inv(A.');v=opt(B);K1=v*norm(A,inf);K2=cond(A,inf);disp(['n=',num2str(n)])disp(['估计条件数为',num2str(K1)])disp(['实际条件数为',num2str(K2)])end计算结果为n=5估计条件数为943656实际条件数为943656n=6估计条件数为29070279.0028实际条件数为29070279.0028n=7估计条件数为985194887.5079实际条件数为985194887.5079n=8估计条件数为.7717实际条件数为.7717n=9估计条件数为86.422实际条件数为86.422n=10估计条件数为750.67实际条件数为750.67n=11估计条件数为49344实际条件数为49344Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.547634e-17.> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.547634e-17.> In cond at 47In ex2_1 at 6n=12估计条件数为3.3713e+16实际条件数为3.3713e+16Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.847381e-19.Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.847381e-19.> In cond at 47In ex2_1 at 6n=13估计条件数为1.5327e+18实际条件数为1.5327e+18Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.246123e-18.> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.246123e-18.> In cond at 47In ex2_1 at 6n=14估计条件数为4.8374e+17实际条件数为4.8374e+17Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.491876e-19.> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.491876e-19.> In cond at 47In ex2_1 at 6n=15估计条件数为4.9674e+17实际条件数为5.3619e+17Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 9.137489e-19.> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 9.137489e-19.> In cond at 47In ex2_1 at 6n=16估计条件数为8.3166e+17实际条件数为8.3167e+17Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.244518e-19.> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.244518e-19.> In cond at 47n=17估计条件数为1.093e+18 实际条件数为1.093e+18Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.693737e-19. > In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.693737e-19. > In cond at 47 In ex2_1 at 6 n=18估计条件数为2.0651e+18 实际条件数为2.7893e+18Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.264685e-19. > In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.264685e-19. > In cond at 47 In ex2_1 at 6 n=19估计条件数为2.9357e+18 实际条件数为2.9357e+18Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.351364e-19. > In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.351364e-19. > In cond at 47 In ex2_1 at 6 n=20估计条件数为2.674e+18 实际条件数为6.473e+18结果分析随着矩阵阶数增加，估计值误差开始出现，20,17,16,15=n 时估计条件数与实际值存在误差；且条件数很大，Hilbert 矩阵为病态的。

数值线性代数实验大报告指导老师：赵国忠姓名：1108300001 刘帅1108300004 王敏1108300032 郭蒙一、实验名称：16题P75上机习题二、实验目的：编制通用的子程序，完成习题的计算任务三、实验内容与要求：P75上机习题先用熟悉的计算机语言将算法2.5.1编制成通用的子程序，然后再用所编制的子程序完成下面两个计算任务：（1） 估计5到20阶Hilbert 矩阵的无穷范数条件数。

（2） 设A n = 11...111................1-1 (01)-- 先随机地选取x ∈R n ，并计算出b=An x;然后再用列主元Gauss 消去法求解该方程组，假定计算解为∧x .试对n 从5到30估计计算解∧x 的精度，并且与真实的相对误差作比较。

四、 实验原理：（1）矩阵范数（martix norm ）是数学上向量范数对矩阵的一个自然推广。

利用for循环和cond （a ）Hilbert 求解Hilbert 矩阵的无穷范数，再利用norm(a,inf)求矩阵的无穷范数条件数。

（2）本题分为4步来求解。

先运用rand 随机选取x ∈R n，输入A n 矩阵，编制一个M 文件计算出b 。

第二步用列主元高斯消去法求解出方程的解X2。

第三步建立M 文件: soluerr.m 估计计算解∧x 的精度。

第四步, 建立M 文件: bijiao.m ,与真实相对误差作比较。

五、 实验过程：(1)程序：clearfor n=5:20for i=1:nfor j=1:na(i,j)=1/(i+j-1);endendc=cond(a);f=norm(c,inf);fprintf('n=%3.0f\nnorm(c,inf)%e\n',n,f) end运行结果：n= 5norm(c,inf)4.766073e+005n= 6norm(c,inf)1.495106e+007n= 7norm(c,inf)4.753674e+008n= 8norm(c,inf)1.525758e+010n= 9norm(c,inf)4.931542e+011n= 10norm(c,inf)1.602467e+013n= 11norm(c,inf)5.224376e+014n= 12norm(c,inf)1.698855e+016n= 13norm(c,inf)3.459404e+017n= 14norm(c,inf)4.696757e+017n= 15norm(c,inf)2.569881e+017n= 16norm(c,inf)7.356249e+017n= 17norm(c,inf)4.362844e+017n= 18norm(c,inf)1.229633e+018n= 19norm(c,inf)9.759023e+017n= 20norm(c,inf)1.644051e+018（2）程序：M文件：matrix1.mfunction [a,b,x1]=matrix1(n) format longA1=-1*ones(n,n)A2=tril(A1)for i=1:nA2(i,i)=1endA2(:,n)=1a=A2x1=rand(n,1)b=A2*x1end运行结果：>> A1 =-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1A2 =-1 0 0 0 0 -1 -1 0 0 0 -1 -1 -1 0 0 -1 -1 -1 -1 0 -1 -1 -1 -1 -1A2 =1 0 0 0 0 -1 -1 0 0 0 -1 -1 -1 0 0 -1 -1 -1 -1 0 -1 -1 -1 -1 -1A2 =1 0 0 0 0-1 -1 -1 0 0 -1 -1 -1 -1 0 -1 -1 -1 -1 -1A2 =1 0 0 0 0 -1 1 0 0 0 -1 -1 1 0 0 -1 -1 -1 -1 0 -1 -1 -1 -1 -1A2 =1 0 0 0 0 -1 1 0 0 0 -1 -1 1 0 0 -1 -1 -1 1 0 -1 -1 -1 -1 -1A2 =1 0 0 0 0 -1 1 0 0 0 -1 -1 1 0 0 -1 -1 -1 1 0 -1 -1 -1 -1 1A2 =1 0 0 0 1 -1 1 0 0 1 -1 -1 1 0 1 -1 -1 -1 1 1 -1 -1 -1 -1 1a =-1 1 0 0 1 -1 -1 1 0 1 -1 -1 -1 1 1 -1 -1 -1 -1 1x1 =0.8147236863931790.9057919370756190.1269868162935060.9133758561390190.632359246225410b =1.4470829326185890.723427496907850-0.961169560949882-0.301767337397875-2.128519049675914a =1 0 0 0 1 -1 1 0 0 1 -1 -1 1 0 1 -1 -1 -1 1 1 -1 -1 -1 -1 1b =1.4470829326185890.723427496907850-0.961169560949882-0.301767337397875-2.128519049675914x1 =0.8147236863931790.9057919370756190.1269868162935060.9133758561390190.632359246225410M文件：LZYgauss.mfunction[x2]=LZYgauss(a,b)format longn=length(a);x2=zeros(n,1);a=[a b];for k=1:n-1max=k;for i=k+1:nif a(i,k)>a(max,k)max=i;endendtemp=a(k,k:n+1);a(k,k:n+1)=a(max,k:n+1);a(max,k:n+1)=temp;for i=k+1:na(i,k)=-a(i,k)/a(k,k);a(i,k+1:n+1)=a(i,k+1:n+1)+a(i,k)*a(k,k+1:n+1);endendx2(n,1)=a(n,n+1)/a(n,n);for i=n-1:-1:1sum=0;for j=i+1:nsum=sum+x2(j,1)*a(i,j);endx2(i,1)=(a(i,n+1)-sum)/a(i,i);end运行结果：>> LZYgauss(a,b)ans =0.8147236863931790.9057919370756190.1269868162935060.913375856139020 0.632359246225410估计计算解x的精度：M文件: soluerr.mfunction [x,error]=soluerr(a,b)format long%估计计算解的精度% 算法：列主元Gauss消去法，其中% A --- 系数矩阵% b-右端项% index --- index=0表示计算成败；index=1表示计算成功%输出结果：error--本算法给出的计算解的估计% normA--逆矩阵无穷范数估计% rnorm--计算解的残量[n,m]=size(a); nb=length(b);if n~=merror('The rows and columns of matrix a must be equal!');return;endif m~=nberror('The columns of a must be equal the dimension of b!'); return;endindex=1;%列主元矩阵三角分解[L,U,u,index_col]=Gauss_col(a);%解下三角方程组Ly=Pb[y,index_low]=Gauss_low(L,b(u));%解上三角方程组Ux=y[x,index_upp]=Gauss_upp(U,y);%输出数值解xpause(0.3)%估计矩阵逆的无穷大范数normA=normAinv(L,U,u);%估计计算解的残量rnorm=norm(b-a*x,inf);%计算右端项bnorm=norm(b,inf);%计算矩阵A的范数Anorm=norm(a,inf);%计算解的精度error=normA*Anorm*rnorm/bnorm;运行结果：>> [x,error]=soluerr(a,b)x =0.8147236863931790.9057919370756190.1269868162935060.9133758561390200.632359246225410x =0.8147236863931790.9057919370756190.1269868162935060.9133758561390200.632359246225410error =5.215941218821592e-016X1与真实相对误差做比较M文件：bijiao.mfunction [x,x1,error,error1]=bijiao(a,b,n) [a,x1]=matrix1(n)[x,error]=soluerr(a,b)error1=abs((x-x1)/x1)运行结果：>>x =0.8147236863931790.9057919370756190.1269868162935060.9133758561390200.632359246225410x1 =0.9578935505663481.6132601689718680.629241553693582-0.799716694069468-1.770467991515150error =5.215941218821592e-016error1 =Columns 1 through 40 0 0 00 0 0 00 0 0 00 0 0 00 0 0 0Column 50.0808655478999350.3995939126190040.2836847318376260.9675930648949141.357170674226221六、结果分析：1、矩阵范数（martix norm）是数学上向量范数对矩阵的一个自然推广。

数值线性代数实验题目：数值线性代数专业：信息与计算科学班级：班姓名：山东科技大学2013年 1 月16日实验报告说明学院：信息学院专业：信息班级10-2 姓名：一、主要参考资料：（1）《Matlab数值计算-案例分析》北京航空出版（2）《Matlab数值分析》机械工业出版二、课程设计应解决的主要问题：（1）平方根（2）QR方法（3）最小二乘法三、应用软件：（1）Matlab7.0（2）数学公式编辑器四、发出日期：课程设计完成日期：指导教师签字：系主任签字：指导教师对课程设计的评语指导教师签字：年月日一、问题描述先用你所熟悉的计算机语言将平方根和改进的平方根法编成写通用的子程序，然后用你编写的程序求解对称正定方程组b x =A ，其中 （1）b 随机的选取，系数矩阵位100阶矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡1011101110111011101110（2）系数矩阵为40阶Hilbert 矩阵，即系数矩阵A 的第i 行第j 列元素为11-+=j i a ij ，向量b 的第i 个分量为∑=-+=nj i j i b 111。

二、分析与程序1. 平方根法函数程序如下：function [x,b]=pingfanggenfa(A,b) n=size(A); n=n(1);x=A^-1*b; disp('Matlab 自带解即为x'); for k=1:nA(k,k)=sqrt(A(k,k));A(k+1:n,k)=A(k+1:n,k)/A(k,k); for j=k+1:n;A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k); endend for j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n);A=A';for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1);disp('平方根法的解即为b');endfunction [x]=ave(A,b,n)求解Ax=bL=zeros(n,n);D=diag(n,0);S=L*D;for i=1:n %L的主对角元素均为1L(i,i)=1;endfor i=1:nfor j=1:nif (eig(A)<=0)|(A(i,j)~=A(j,i))disp('wrong');break;endendendD(1,1)=A(1,1);for i=2:nfor j=1:i-1S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)');L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1);endD(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1)');endy=zeros(n,1);x=zeros(n,1);for i=1:ny(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)))/D(i,i); endfor i=n:-1:1x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n));end2.改进平方根法函数程序如下：function b=gaijinpinfanggenfa(A,b)n=size(A);n=n(1);v=zeros(n,1);for j=1:nfor i=1:j-1v(i)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1))/A(j,j);end %LDL'分解B=diag(A);D=zeros(n);for i=1:nD(i,i)=B(i);A(i,i)=1;EndA=tril(A);for j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j); endb(n)=b(n)/A(n,n);A=D*(A');for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1);disp('改进平方根法解得的解即为b'); end3.调用函数解题：clear;clc;n=input('请输入矩阵维数：');b=zeros(n,1); A=zeros(n);for i=1:nfor j=1:nA(i,j)=1/(i+j-1);b(i)=b(i)+1/(i+j-1);endend[x,b]=pingfanggenfa(A,b)b=gaijinpinfanggenfa(A,b)4.运行结果：请输入矩阵维数：40Warning: Matrix is close to singular or badly scaled.Results may be inaccurate. RCOND = 6.570692e-020. > In pingfanggenfa at 4In qiujie at 10Matlab自带解即为x平方根法的解即为bx =1.60358.96850.85621.01950.9375-50.2500-3.0000-16.000024.0000-49.5000-30.000039.000022.0000-64.0000 -12.00002.000010.2500 -10.5000-1.0000 -10.875083.000046.0000 -98.000012.0000 -69.000068.000021.000017.0000 -50.7188-8.7500-8.0000 112.00006.0000 -68.750022.000044.0000 -28.00008.0000 -44.000012.0000b =1.0e+007 *0.0000-0.00000.0001-0.0004-0.00140.0424-0.29801.1419-2.73354.2539-4.30182.7733-1.19890.5406-0.36880.3285-0.44380.4621-0.25130.05650.0000-0.00510.0071-0.0027-0.00310.0036-0.00190.00090.0002-0.0002-0.00060.00040.0001-0.00020.00010.0000-0.00000.0000-0.0000-0.0000改进平方根法解得的解即为bb =1.0e+024 *0.0000-0.00000.0001-0.00120.0139-0.09540.4208-1.21012.0624-1.0394-3.33436.2567-0.2463-7.45942.80303.69900.7277-1.7484-0.4854-3.60100.25325.1862-2.12991.44100.8738-4.56541.04224.0920-2.7764-2.2148-0.89530.36654.89671.04160.1281 -4.3387 -1.1902 -2.8334 8.4610 -3.6008一、问题描述先用你所熟悉的计算机语言将算法2.5.1编成写通用的子程序，然后用你编写的程序完成下面两个计算任务：（1） 估计5到20阶Hilbert 矩阵的∞范数条件数；（2）设n *n n R 11-1-1-111-1-101-1001A ∈⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=先随机的选取n R x ∈，并计算出x A b n =；然后再用列主元Gauss 消去法求解该方程组，假定计算解为xˆ。

线性代数方程组的数值解法实验1. 主元的选取与算法的稳定性问题提出：Gauss 消去法是我们在线性代数中已经熟悉的。

但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss 消去法作为数值算法的稳定性呢？Gauss 消去法从理论算法到数值算法，其关键是主元的选择。

主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。

实验内容：考虑线性方程组 n n n R b R A b Ax ∈∈=⨯,,编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss 消去过程。

实验要求：（1）取矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157,6816816816 b A ，则方程有解T x )1,,1,1(* =。

取n=10计算矩阵的条件数。

让程序自动选取主元，结果如何？（2）现选择程序中手动选取主元的功能。

每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。

若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。

（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。

（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。

重复上述实验，观察记录并分析实验结果。

1.1程序清单n=input('矩阵A 的阶数:n=');A=6*diag(ones(1,n))+diag(ones(1,n-1),1)+8*diag(ones(1,n-1),-1); b=A*ones(n,1);p=input('计算条件数使用p-范数,p='); cond_A=cond(A,p) [m,n]=size(A); Ab=[A b];r=input('选主元方式(0:自动;1:手动)，r=');Abfor i=1:n-1switch rcase(0)[aii,ip]=max(abs(Ab(i:n,i)));ip=ip+i-1;case (1)ip=input(['第',num2str(i),'步消元，请输入第',num2str(i),'列所选元素所处的行数：']);end;Ab([i ip],:)=Ab([ip i],:);aii=Ab(i,i);for k=i+1:nAb(k,i:n+1)=Ab(k,i:n+1)-(Ab(k,i)/aii)*Ab(i,i:n+1);end;if r==1Abendend;x=zeros(n,1);x(n)=Ab(n,n+1)/Ab(n,n);for i=n-1:-1:1x(i)=(Ab(i,n+1)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);endx1.2运行结果(1)n=10，矩阵的条件数及自动选主元Cond(A,1) =2.5575×103Cond(A,2) = 1.7276×103Cond(A,inf) =2.5575×103程序自动选择主元（列主元）a.输入数据矩阵A的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=0b.计算结果x=[1，1，1，1，1，1，1，1，1，1]T (2)n=10，手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元 矩阵A 的阶数:n=10 计算条件数使用p-范数,p=1 选主元方式(0:自动;1:手动)，r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元，请输入第1列所选元素所处的行数：1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元，请输入第2列所选元素所处的行数：2 …（实际选择时，第k 步选择主元处于第k 行） 最终计算得x=[1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000001, 0.999999999999998, 1.000000000000004, 0.999999999999993, 1.000000000000012, 0.999999999999979, 1.000000000000028]Tb. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=10 计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元，请输入第1列所选元素所处的行数：2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元，请输入第2列所选元素所处的行数：3 …（实际选择时，第k 步选择主元处于第k+1行） 最终计算得x=[1，1，1，1，1，1，1，1，1，1]T (3)n=20，手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元 矩阵A 的阶数:n=20 计算条件数使用p-范数,p=1 选主元方式(0:自动;1:手动)，r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元，请输入第1列所选元素所处的行数：1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元，请输入第2列所选元素所处的行数：2 …（实际选择时，第k 步选择主元处于第k 行） 最终计算得x=[1.000000000000000，1.000000000000000，1.000000000000000，1.000000000000001，0.999999999999998，1.000000000000004，0.999999999999993，1.000000000000014，0.999999999999972，1.000000000000057，0.999999999999886，1.000000000000227，0.999999999999547，1.000000000000902，0.999999999998209，1.000000000003524，0.999999999993179，1.000000000012732，0.999999999978173，1.000000000029102]T b. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=20 计算条件数使用p-范数,p=1 选主元方式(0:自动;1:手动)，r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元，请输入第1列所选元素所处的行数：2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元，请输入第2列所选元素所处的行数：3…（实际选择时，第k步选择主元处于第k+1行）最终计算得x=[1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1]T(4)A分别为幻方矩阵，Hilbert矩阵，pascal矩阵和随机矩阵将计算结果列于下表：1.3简要分析计算(1)表明：对于同一矩阵，不同范数定义的条件数是不同的；Gauss消去法在消去过程中选择模最大的主元能够得到比较精确的解。

数值代数课程设计实验报告姓名： 班级： 学号： 实验日期：一、实验名称 代数的数值解法 二、实验环境实验一、平方根法与改良平方根法一、实验要求：用熟悉的运算机语言将不选主元和列主元Gasuss 消元法编写成通用的子程序，然后用编写的程序求解以下方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--⨯14151515157681681681681681612321n n n n n x x x x x x 用所编的程序别离求解40、84、120阶方程组的解。

二、算法描述及实验步骤GAuss 程序如下：（1）求A 的三角分解：LU A =;（2）求解b y =L 得y ； （3）求解y x =U 得x ；列主元Gasuss 消元法程序如下： 1求A 的列主元分解：LU PA =;2求解b y P L =得y ；3求解y x U 得x ；三、调试进程及实验结果：%----------------方程系数---------------->> A1=Sanduijiaozhen(8,6,1,40); >> A2=Sanduijiaozhen(8,6,1,84); >> A3=Sanduijiaozhen(8,6,1,120); >> b1(1)=7;b2(1)=7;b3(1)=7; >> for i=2:39 b1(i)=15; end>> b1(40)=14; >> for i=2:83 b2(i)=15; end>> b2(40)=14; >> for i=2:119 b1(i)=15; end>> b3(120)=14;%----------------方程解---------------->> x11=GAuss(A1,b1') >> x12=GAuss Zhu (A1,b1') >> x21=GAuss(A2,b2') >> x22=GAuss Zhu (A3,b3') >> x31=GAuss(A3,b3') >> x32=GAuss Zhu (A3,b3') 运行结果：（n=40） GAuss 消元法的解即为x11 =列主元GAuss消元法的解即为x12 =1 1 1 1 1 1 1 1 1 1 111111111111111111111111111111六、源程序：function A=Sanduijiaozhen(a,b,c,n)%生成n阶以a,b,c为元素的三对角阵A=diag(b*ones(1,n),0)+diag(c*ones(1,n-1),1)+diag(a*ones(1,n-1),-1);function x=GAuss(A,b)n=length(b);x=b;%-------分解---------------for i=1:n-1for j=i+1:nmi=A(j,i)/A(i,i);b(j)=b(j)-mi*b(i);for k=i:nA(j,k)=A(j,k)-mi*A(i,k);endAB=[A,b];endend%-----------回代------------------x(n)=b(n)/A(n,n);for i=n-1:-1:1s=0;for j=i+1:ns=s+A(i,j)*x(j);endx(i)=(b(i)-s)/A(i,i);endfunction x=GAussZhu(A,b)n=length(b);x=b;%----------------------选主元---------------------for k=1:n-1a_max=0;for i=k:nif abs(A(i,k))>a_maxa_max=abs(A(i,k));r=i;endendif r>kfor j=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;end%--------------消元-----------------for i=k+1:nm=A(i,k)/A(k,k);for j=k:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);endendif abs(A(n,n))==0return;endAbZhu=[A,b];%----------------回代-----------------------x(n)=b(n)/A(n,n);for i=n-1:-1:1for j=i+1:nb(i)=b(i)-A(i,j)*x(j);endx(i)=b(i)/A(i,i);end实验二、平方根法与改良平方根法一、实验要求：用运算机语言将平方根法和改良的平方根法编成通用的子程序，然后用编写的程序求解对称正定方程组100阶方程组AX=b，二、算法描述及实验步骤：平方根法函数程序如下：一、求A 的Cholesky 分解：L L A T=; 二、求解b y =L 得y ；3、求解y x =TL 得x ；改良平方根法函数程序如下：一、求A 的Cholesky 分解：T=LDL A ; 二、求解b y =L 得y ； 3、求解y x =TDL 得x ；三、调试进程及实验结果：clear;clc;%----------------方程系数---------------->> A=Sanduijiaozhen(1,10,1,100); >> b(1)=11; >> for i=2:99 b(i)=12; end>> b(100)=11;>> x1=Cholesky(A,b') >> x2=GJCholesky(A,b')运行结果：平方根法的解即为 x1 =改良平方根法解得的解即为x2 =四、源程序：function x=Cholesky(A,b)n=size(A);n=n(1);% x=A^-1*b;% disp('Matlab自带解即为x');%-----------------Cholesky分解-------------------for k=1:nA(k,k)=sqrt(A(k,k));A(k+1:n,k)=A(k+1:n,k)/A(k,k);for j=k+1:n;A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endend%------------------前代法求解Ly=b----------------------------for j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n);%-----------------回代法求解L'x=y-----------------------------A=A';for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1);disp('平方根法的解即为');function b=GJCholesky(A,b)n=size(A);n=n(1);v=zeros(n,1);%----------------------LDL'分解-----------------------------for j=1:nfor i=1:j-1v(i)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1))/A(j,j);endB=diag(A);D=zeros(n);for i=1:nD(i,i)=B(i);A(i,i)=1;end%-------------------前代法---------------------------A=tril(A); %取得L和Dfor j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n);%-----------------回代法-----------------------------A=D*(A');for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1);disp('改良平方根法解得的解即为');实验三、二次多项式拟合一、实验要求：用运算机语言编制利用QR分解求解线性最小二乘问题的通用子程序，用编写的程序求解一个二次多项式使在残向量的范数最小的意义下拟合下面的数据t-1 0iyi二、算法描述及实验步骤：QR分解求解程序如下：一、求A 的QR 分解; 二、计算b c 11T=Q ;3、求解上三角方程1c x =R 得x ；三、调试进程及实验结果：>> t=[-1 0 ]; >> y=[ ]; >> plot(t,y,'r*');>> legend('实验数据(ti,yi)'); >> xlabel('t'), ylabel('y');>> title('二次多项式拟合的数据点(ti,yi)的散点图');运行后屏幕显示数据的散点图（略）.（3）编写以下MATLAB 程序计算)(x f 在),(i i y x 处的函数值，即输入程序 >> syms a b c >> t=[-1 0 ]; >> fi=a.*t.^2+ b.*t+c%运行后屏幕显示关于 ,,a b c 的线性方程组fi =[a-b+c,9/16*a-3/4*b+c,1/4*a-1/2*b+c,c,1/16*a+1/4*b+c,1/4*a+1/2*b+c,9/16*a+3/4*b +c]编写构造残向量2范数的MATLAB 程序>> y=[ ]; >> y=[ ];>> fy=fi-y; fy2=fy.^2; J=sum(fy.^2); 运行后屏幕显示误差平方和如下 J=(a-b+c-1)^2+(9/16*a-3/4*b+c-13/16)^2+(1/4*a-1/2*b+c-3/4)^2+(c-1)^2+(1/16*a+1/4*b+c-21/16)^2+(1/4*a+1/2*b+c-7/4)^2+(9/16*a+3/4*b+c-37/16)^2为求,,a b c 使J 达到最小，只需利用极值的必要条件0J a ∂=∂，0J b ∂=∂，0Jc∂=∂，取得关于,,a b c 的线性方程组，这能够由下面的MATLAB 程序完成，即输入程序 >> Ja1=diff(J,a); Ja2=diff(J,b); Ja3=diff(J,c);>> Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3) 运行后屏幕显示J 别离对,,a b c 的偏导数如下 Ja11 =451/128*a-63/32*b+43/8*c-887/128 Ja21 =-63/32*a+43/8*b-3/2*c-61/32Ja31 =43/8*a-3/2*b+14*c-143/8解线性方程组112131000Ja Ja Ja ===，，，输入以下程序 >> A=[451/128, -63/32, -3/2 ;-63/32,43/8,-3/2;43/8,-3/2,14]; >> B=[887/128,61/32,143/8]; >> C=B/A, f=poly2sym(C)运行后屏幕显示拟合函数f 及其系数C 如下 C =f =924/2999*x^2+10301/11996*x+4204/2999 故所求的拟合曲线为2()0.30810.8581 1.4018f x x x =++四、源程序：>> t=[-1 0 ]; >> y=[ ]; >> plot(t,y,'r*');>> legend('实验数据(ti,yi)'); >> xlabel('t'), ylabel('y');>> title('二次多项式拟合的数据点(ti,yi)的散点图'); >> syms a b c>> t=[-1 0 ]; >> fi=a.*t.^2+ b.*t+c fi =[ a-b+c, 9/16*a-3/4*b+c, 1/4*a-1/2*b+c, c, 1/16*a+1/4*b+c, 1/4*a+1/2*b+c, 9/16*a+3/4*b+c]>> y=[ ]; >> y=[ ];>> fy=fi-y; fy2=fy.^2; J=sum(fy.^2) J =(a-b+c-1)^2+(9/16*a-3/4*b+c-13/16)^2+(1/4*a-1/2*b+c-3/4)^2+(c-1)^2+(1/16*a+1/4*b+c-21/16)^2+(1/4*a+1/2*b+c-7/4)^2+(9/16*a+3/4*b+c-37/16)^2>> Ja1=diff(J,a); Ja2=diff(J,b); Ja3=diff(J,c);>> Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3)Ja11 =451/128*a-63/32*b+43/8*c-887/128Ja21 =-63/32*a+43/8*b-3/2*c-61/32Ja31 =43/8*a-3/2*b+14*c-143/8>> A=[451/128, -63/32, -3/2 ;-63/32,43/8,-3/2;43/8,-3/2,14]; >> B=[887/128,61/32,143/8];>> C=B/A, f=poly2sym(C)C =f =924/2999*x^2+10301/11996*x+4204/2999>>。