基于Matlab的Cramer法则求解线性方程组

利用Matlab进行线性代数问题求解的方法与案例引言线性代数是数学的一个重要分支，广泛应用于工程、物理、计算机科学等领域。

而Matlab作为一种功能强大的数值计算软件，提供了各种实用的工具和函数，可以方便地解决线性代数问题。

本文将介绍一些常用的线性代数问题求解方法，并通过具体的案例来展示Matlab在实际应用中的效果。

一、线性方程组的求解线性方程组是线性代数中最基础的问题之一。

Matlab提供了多种求解线性方程组的函数，如“backslash”操作符（\）和“linsolve”函数等。

下面通过一个实例来说明Matlab的线性方程组求解功能。

案例：假设有以下线性方程组需要求解：2x + 3y - 4z = 53x - 2y + z = 8x + 5y - 3z = 7在Matlab中输入以下代码：A = [2 3 -4; 3 -2 1; 1 5 -3];b = [5; 8; 7];x = A\b;通过以上代码，我们可以得到线性方程组的解x=[1; -2; 3]。

这表明在满足以上方程组的条件下，x=1，y=-2，z=3。

可以看出，Matlab在求解线性方程组时，使用简单且高效。

二、矩阵的特征值和特征向量求解矩阵的特征值和特征向量也是线性代数中的重要概念。

利用特征值和特征向量可以得到矩阵的许多性质和信息。

在Matlab中，我们可以通过“eig”函数来求解矩阵的特征值和特征向量。

案例：假设有一个2x2矩阵A，需要求解其特征值和特征向量。

在Matlab中输入以下代码：A = [2 3; 1 4];[V, D] = eig(A);通过以上代码，我们可以得到矩阵A的特征向量矩阵V和特征值矩阵D。

具体结果如下：特征向量矩阵V = [0.8507 -0.5257; 0.5257 0.8507]特征值矩阵D = [1.5858 0; 0 4.4142]由结果可知，矩阵A的特征向量矩阵V和特征值矩阵D可以提供有关该矩阵的很多信息，如相关线性变换、对称性等。

一、介绍Crout分解法Crout分解法是一种用于解决线性代数问题的数值计算方法。

它可以将一个矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积，从而简化线性方程组的求解过程。

Crout分解法的基本思想是利用矩阵的三角形式，将原线性方程组的求解问题转化为两个三角方程的求解问题。

二、Crout分解法的原理1. 初等行变换Crout分解法的第一步是对系数矩阵进行初等行变换，将其变换为一个上三角矩阵和一个下三角矩阵的乘积。

在进行初等行变换时，需要注意保持方程组的等价性，以便最终求解得到的解与原方程组的解是一致的。

2. LU分解经过初等行变换后，原方程组的系数矩阵可以表示为两个三角矩阵的乘积，即LU分解，其中L是下三角矩阵，U是上三角矩阵。

对于一个n阶矩阵A，LU分解的形式为A=LU。

Crout分解法通过逐步迭代，可以求解L和U矩阵的具体数值。

3. 解方程组经过LU分解后，原线性方程组可以转化为两个三角方程的求解问题。

首先求解下三角方程LUx=b，得到中间向量y，然后再求解上三角方程Ux=y，得到最终的解x。

通过这种分解和迭代的方法，可以有效地求解复杂的线性方程组。

三、Crout分解法的实现1. Matlab实现在Matlab中，可以使用内置的lu函数来实现Crout分解法。

lu函数可以对一个矩阵进行LU分解，返回分解后的下三角矩阵L和上三角矩阵U。

通过这两个矩阵，可以进一步求解原线性方程组。

2. 算法步骤（1）输入系数矩阵A和常数向量b；（2）调用lu函数，对矩阵A进行LU分解，得到下三角矩阵L和上三角矩阵U；（3）将原线性方程组转化为两个三角方程的求解问题，分别求解Ly=b和Ux=y，得到最终的解x。

3. 代码示例```Matlab定义系数矩阵A和常数向量bA = [1 2 3; 4 5 6; 7 8 9];b = [1; 2; 3];调用lu函数进行LU分解[L, U] = lu(A);求解Ly=b和Ux=y，得到最终解xy = L\b;x = U\y;打印解xdisp(x);```四、Crout分解法的应用Crout分解法在实际问题中有着广泛的应用，特别是在工程和科学计算领域。

克拉默法则解方程组一、什么是克拉默法则？克拉默法则（Cramer's Rule）是一种用于解决线性方程组的数学算法。

它通过将方程组转换为一个特殊的行列式，然后求解行列式来解决方程组。

二、克拉默法则的步骤1.确定方程组的变量：首先，我们需要确定方程组的变量。

例如，我们考虑下面的方程组：2x+3y=44x+5y=6这里有两个变量x和y。

2.编写方程组的行列式：接下来，我们需要编写方程组的行列式。

我们可以用下面的矩阵表示方程组：[2 3][4 5]然后，我们可以计算这个矩阵的行列式，它是2×5-3×4=2。

3.计算变量的行列式：接下来，我们需要计算变量x和y的行列式。

为了计算变量x的行列式，我们需要把变量x的系数放到矩阵的第一列，把变量y的系数放到矩阵的第二列，然后计算行列式：[2 4][3 5]这里的行列式是2×5-3×4=2，同样的，我们可以计算变量y的行列式：[3 4][2 5]这里的行列式是3×5-2×4=6。

4.计算变量的值：最后，我们可以根据克拉默法则计算变量x和y的值：x=6/2=3y=2/2=1这样，我们就可以得到方程组的解：x=3，y=1。

三、克拉默法则的优点1.计算方便：克拉默法则可以让我们更轻松地解决复杂的线性方程组，而不需要花费大量的时间和精力。

2.结果准确：克拉默法则可以提供准确的结果，因为它是基于行列式的计算。

四、克拉默法则的缺点1.只能用于线性方程组：克拉默法则只能用于解决线性方程组，不能用于解决非线性方程组。

2.计算量大：对于大型方程组，克拉默法则可能需要大量的计算，这可能会耗费大量的时间和精力。

五、例子假设我们有以下方程组：2x+3y=44x+5y=6我们可以使用克拉默法则来解决这个方程组：1.确定方程组的变量：x和y。

2.编写方程组的行列式：[2 3] [4 5]，行列式的值为2。

3.计算变量的行列式：x的行列式为[2 4] [3 5]，值为2；y的行列式为[3 4] [2 5]，值为6。

1
9—2 Matlab 中线性代数方程组的求解
2
z 定解方程组
定解方程组指未知数的个数和方程个数相同的方程
组，其求解可采用矩阵除法和初等变换法，当方程组的系数矩阵非奇异时，还可通过求系数矩阵的逆来进行
z 不定方程组
不定方程组指未知数的个数大于方程个数的方程组，用矩阵除法将给出一个具有最多零元素的特解，其通解则可通过null 命令求一个齐次组的基础解系来得到，或者由初等变换得到，
3z 超定方程组
超定方程组指未知数的个数小于方程个数的方程组，用矩阵除法可求得一最小二乘近似解。

z 奇异方程组
奇异方程组指方程组系数矩阵是奇异的，用矩阵除法不能直接求解，但经同解变形仍可求得一特解。

z 符号方程组
符号方程组的求解可通过函数solve 、linsolve 来实现。

4
z null()z rref()z A\b
z sym(A)\sym(b)z inv()
5
z rank( )z det( )z inv( )z eig( )z norm( )z jordan( )z poly( ) %求特征单项式
%P=poly(A), poly2str(P,’x’)。

科学计算—理论、方法及其基于ＭＡＴＬＡＢ的程序实现与分析 三、 解线性方程组（线性矩阵方程）解线性方程组是科学计算中最常见的问题。

所说的“最常见”有两方面的含义：１） 问题的本身是求解线性方程组；２） 许多问题的求解需要或归结为线性方程组的求解。

关于线性方程组B A x B Ax 1-=⇒=（１）其求解方法有两类：１） 直接法：高斯消去法（Gaussian Elimination ）； ２） 间接法：各种迭代法（Iteration ）。

１、高斯消去法１） 引例考虑如下（梯形）线性方程组：()⎪⎩⎪⎨⎧==+==+-=⇒⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--⇔⎪⎩⎪⎨⎧==-=+-5.0141315.3221122004301211214322332321321332321x x x x x x x x x x x x x x x 高斯消去法的求解思路：把一般的线性方程组（１）化成（上或下）梯形的形式。

２）高斯消去法——示例考虑如下线性方程组：⎪⎪⎪⎭⎫ ⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---⇔⎪⎩⎪⎨⎧=++-=-+-=+-306015129101.2001.221113060129501.2001.221321321321321x x x x x x x x x x x x １） 第一个方程的两端乘12加到第二个方程的两端，第一个方程的两端乘－１加到第三个方程的两端，得⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--3060031110001.0001.00111321x x x２） 第二个方程的两端乘001.010-加到第三个方程的两端，得 ⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--60600311010001.0001.00111321x x x３） 从上述方程组的第三个方程依此求解，得()⎪⎩⎪⎨⎧==+-==+-=600300001.03100024011332321x x x x x x ３）高斯消去法的不足及其改进——高斯（全、列）主元素消去法在上例中，由于建模、计算等原因，系数2.001而产生0.0005的误差，实际求解的方程组为⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---306015129101.20005.22111321x x x ⎪⎩⎪⎨⎧===⇒70.4509.30142.2565321x x x注：数值稳定的算法高斯列主元素消去法就是在消元的每一步选取（列）主元素—一列中绝对值最大的元取做主元素，高斯列主元素消去法是数值稳定的方法。

matlab解方程组方法在MATLAB中，有多种方法可以解方程组。

以下是其中几种常用的方法：1.solve函数：这是最直接的方法，适用于解线性方程组。

假设你有以下线性方程组：(Ax = b)你可以使用solve函数来求解。

例如：2.matlab复制代码A = [1, 2; 3,4];b = [5; 6];x = solve(A,b);3.\和/运算符：这两个运算符也可以用于解线性方程组。

例如：4.matlab复制代码A = [1, 2; 3, 4];b = [5; 6];x = A\b; % 使用左除运算符或者matlab复制代码x = b/A; % 使用右除运算符5.gaussj函数：这个函数使用高斯-约当消元法来解方程组。

使用方法如下：6.matlab复制代码A = [1, 2; 3,4];b = [5; 6];x = gaussj(A,b);7.mldivide函数：这个函数与\运算符相同，也是用于解线性方程组。

例如：8.matlab复制代码A = [1, 2; 3, 4];b = [5; 6];x = mldivide(A, b); % 等价于A\b9.lyap函数：对于非线性方程组，可以使用lyap函数来求解。

这个函数用于解决Lyapunov方程，通常用于控制系统和稳定性分析。

使用方法如下：10.matlab复制代码A = [1, 2; 3, 4];lyap(A); % 对于给定的A矩阵，求解Lyapunov方程。

11.fzero和root函数：这两个函数用于求解非线性方程的根。

例如，如果你有一个非线性方程(f(x) = 0)，你可以使用fzero或root来找到这个方程的根。

使用方法如下：12.matlab复制代码f = @(x) x^2 - 4; % 非线性方程 f(x) = x^2 - 4x = fzero(f, [1, 2]); % 在区间[1,2]内寻找方程的根或者：matlab复制代码root(f) % 使用root函数求解非线性方程的根。

cramer法则的证明
Cramer法则是一种用于解决线性方程组的方法，它可以根据系数矩阵的行列式的值来求解未知变量的值。

证明：
假设给定的线性方程组为：
Ax = b
其中A是n阶系数矩阵，x是n维未知向量，b是n维常数向量。

由于A是n阶系数矩阵，所以它的行列式为Det(A)，它的值可以用下面的公式来表示：
Det(A) = a11a22...ann - a12a21...ann + a13a21...ann - ... + (-
1)^(n+1)ann-1a2a3...a1
假设A的每一行都乘以一个常数c，那么新矩阵的行列式为：
Det(A') = c^n·Det(A)
由于每一行都乘以常数c，所以A'的第i行可以表示为：
a'i = c·ai
这样，更新后的方程组可以表示为：
A'x = b'
其中b' = c·b。

由于A'x = b'，所以未知向量x的每一个分量可以表示为：
xi = (Det(Ai))/Det(A)
其中Ai是由A的第i行替换为b'后得到的矩阵。

由此可以得出Cramer法则的结论：
xi = (Det(Ai))/Det(A)
即每一个未知向量x的分量可以由A的行列式和替换第i行后得到的行列式的比值来求得。

matlab计算方程组Matlab作为一款试用范围广泛的科学计算软件，其计算方程组的能力也是非常强大的。

在Matlab中，可以通过多种方式计算方程组，比如使用直接法、迭代法、线性方程组求解器等等。

下面将分步骤阐述使用Matlab计算方程组的方法。

一、使用直接法求解直接法是一种将系数矩阵直接求逆再与常数向量相乘的方法，通常在方程组的规模较小时使用。

下面是使用Matlab求解线性方程组的示例代码：```matlab% 定义系数矩阵和常数向量A = [1 2 3; 4 5 6; 7 8 9];b = [3; 6; 9];% 求解方程组x = A\b;disp(x);```这段代码首先定义了一个3x3的系数矩阵A和一个3x1的常数向量b，然后使用反斜线符号来求解方程组。

该符号将A的逆矩阵乘上b，得到解向量x。

二、使用迭代法求解当方程组的规模较大时，直接法的计算量可能会非常大，在这种情况下可以使用迭代法来求解方程组。

迭代法的主要思想是通过反复迭代求解来逼近方程组的解。

常见的迭代法有Jacobi迭代法、Gauss-Seidel迭代法等。

以Jacobi迭代法为例，下面是使用Matlab求解线性方程组的示例代码：```matlab% 定义系数矩阵和常数向量A = [1 2 3; 4 5 6; 7 8 9];b = [3; 6; 9];% 定义Jacobi迭代法函数function [x, k] = jacobi(A, b, x0, tol, max_iter)D = diag(diag(A));L = -tril(A, -1);U = -triu(A, 1);x = x0;for k = 1:max_iterx = inv(D)*(b + L*x + U*x);if norm(A*x - b) < tolreturnendendend% 求解方程组x0 = [0; 0; 0];tol = 1e-6;max_iter = 1000;[x, k] = jacobi(A, b, x0, tol, max_iter);disp(x);```这段代码首先定义了一个3x3的系数矩阵A和一个3x1的常数向量b，然后定义了一个Jacobi迭代法的函数来求解方程组。

实验1.1 用matlab 求解线性方程组第一节 线性方程组的求解 一、齐次方程组的求解rref （A ） %将矩阵A 化为阶梯形的最简式null （A ） %求满足AX ＝0的解空间的一组基，即齐次线性方程组的基础解系【例１】 求下列齐次线性方程组的一个基础解系，并写出通解：我们可以通过两种方法来解： 解法1：>> A=[1 -1 1 -1;1 -1 -1 1;1 -1 -2 2]; >> rref(A) 执行后可得结果： ans=1 -1 0 0 0 0 -1 1 0 0 0 0 由最简行阶梯型矩阵，得化简后的方程⎪⎩⎪⎨⎧=+--=+--=-+-02200432143214321x x x x x x x x x x x x取x2，x4为自由未知量，扩充方程组为即提取自由未知量系数形成的列向量为基础解系，记所以齐次方程组的通解为解法2： clearA=[1 -1 1 -1;1 -1 -1 1;1 -1 -2 2];B=null(A, 'r') % help null 看看加个‘r ’是什么作用,若去掉r ，是什么结果？执行后可得结果： B=1 0 1 0 0 1 0 1⎩⎨⎧=-=-004321x x x x ⎪⎪⎩⎪⎪⎨⎧====44432221x x x x x x x x ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡11000011424321x x x x x x ,00111⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=ε,11002⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=ε2211εεk k x +=易见，可直接得基础解系所以齐次方程组的通解为二、非齐次线性方程组的求解 Matlab 命令的基本格式：X ＝A\b %系数阵A 满秩时，用左除法求线性方程组AX ＝b 的解注意：A/B 即为AB -1, 而A\B 即为A -1B.C =[A,b];D ＝rref(C) % 求线性方程组AX ＝b 的特解，即D 的最后一列元素【例2】 求下列非齐次线性方程组的解：,00111⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=ε,11002⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=ε⎪⎪⎪⎩⎪⎪⎪⎨⎧=+=++=++=++=+150650650651655454354332121x x x x x x x x x x x x x 2211εεk k x +=解： clearA=[5 6 0 0 0;1 5 6 0 0;0 1 5 6 0;0 0 1 5 6;0 0 0 1 5]; b=[1;0;0;0;1];format rational %采用有理数近似输出格式,比较format short 看看x=A\b执行后可得所求方程组的解. 作业：【第一题】 求下列非齐次线性方程组的通解.A=[1 2 3 1;1 4 6 2;2 9 8 3;3 7 7 2] B=[3;2;7;12] format rational x=A\B x =⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++=+++1227737389222643324321432143214321x x x x x x x x x x x x x x x x42/31/2684838239393950-7/3【第二题】计算工资问题一个木工，一个电工，一个油漆工，三个人相互同意彼此装修他们自己的房子。

第六章：函数，极限与连续的MATLAB1 映射与函数。

（1）集合（更多的是用于数组间的运算）：ismember（一个个元素判断是否是子集，返回一个数组）；intersect（求交集，返回结果数组）；setdiff(a,b)（求差集，属于a不属于b的数组）；union （求并集）。

（2）函数：定义方法：y=@(x)f(x)；syms x y=f(x)；y=sym(‘f(x)’);求反函数：finverse(f,t)；求复合函数f(g(x))：y=compose(f,g)；2 求极限。

（1）求数列极限：limit(xn, n, inf)；limit(xn, inf)。

（2）求函数极限：limit(fx, x, x0(, ‘left’) )；limit(fx, x, inf)。

3 函数的连续性与间断点。

（1）判断连续性的函数代码：P144。

（2）判断x0是否是函数f(x)的间断点的函数代码：（P146，文件夹MATLAB学习中的程序储存里）。

实际应用中，可以根据绘图来判定是否是间断点。

（3）求函数区间的方法：P215。

第七章：导数与微分的MATLAB求解1 导数求解：diff(fx,x,n)后面2个可以省略，则是求导函数；隐函数的导数求解见P156的2个例子；稍微总结就是把y定义为y=sym(‘y(x)’)，然后定义隐函数的表达式为F=…，把表达式等号右侧置为0，左侧为F函数表达式，之后：diff(F,x)。

参数方程确定的函数的导数P157。

2 洛必达法则：P168.3 泰勒公式：P172.另外，MATLAB有taylor(fx,x,n,a)。

MATLAB提供了泰勒级数逼近分析界面：taylortool，4 函数的凹凸性与曲线的单调性：求函数单调区间及各个区间单调性的判定：P175。

求凹凸性与拐点的程序：P179。

求方程实根从而可以进行一些特殊数值表达式的求解（比如（-8）^(1/3)的求解）的函数代码：P176。

在这章中我们要学习线性方程组的直接法，特别是适合用数学软件在计算机上求解的方法.3.1 方程组的逆矩阵解法及其MATLAB 程序3.1.3 线性方程组有解的判定条件及其MATLAB 程序 判定线性方程组A n m ⨯b X =是否有解的MATLAB 程序function [RA,RB,n]=jiepb(A,b)B=[A b];n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0,disp('请注意：因为RA~=RB ，所以此方程组无解.') return endif RA==RB if RA==ndisp('请注意：因为RA=RB=n ，所以此方程组有唯一解.') elsedisp('请注意：因为RA=RB<n ，所以此方程组有无穷多解.') end end例3.1.4 判断下列线性方程组解的情况.如果有唯一解，则用表 3-2方法求解.(1) ⎪⎪⎩⎪⎪⎨⎧=-+-=+-+=-++=+-+;0742,0634,0723,05324321432143214321x x x x x x x x x x x x x x x x (2) ⎪⎪⎩⎪⎪⎨⎧=++-=+-+=-+-=+-+;0327,01613114,02332,075434321432143214321x x x x x x x x x x x x x x x x (3) ⎪⎩⎪⎨⎧=+=+-=-+;8311,1023,22421321321x x x x x x x x (4) ⎪⎩⎪⎨⎧=--+=+-+=+-+.12,2224,12w z y x w z y x w z y x解 在MATLAB 工作窗口输入程序>> A=[2 3 -1 5;3 1 2 -7;4 1 -3 6;1 -2 4 -7]; b=[ 0; 0; 0; 0]; [RA,RB,n]=jiepb(A,b)运行后输出结果为请注意：因为RA=RB=n ，所以此方程组有唯一解. RA = 4,RB =4,n =4 在MATLAB 工作窗口输入>>X=A\b,运行后输出结果为 X =(0 0 0 0)’.(2) 在MATLAB 工作窗口输入程序>> A=[3 4 -5 7;2 -3 3 -2;4 11 -13 16;7 -2 1 3];b=[ 0; 0; 0; 0];[RA,RB,n]=jiepb(A,b)运行后输出结果请注意：因为RA=RB<n ，所以此方程组有无穷多解. RA =2,RB =2,n =4(3) 在MATLAB 工作窗口输入程序>> A=[4 2 -1;3 -1 2;11 3 0]; b=[2;10;8]; [RA,RB,n]=jiepb(A,B)运行后输出结果请注意：因为RA~=RB ，所以此方程组无解. RA =2,RB =3,n =3(4)在MATLAB 工作窗口输入程序>> A=[2 1 -1 1;4 2 -2 1;2 1 -1 -1]; b=[1; 2; 1]; [RA,RB,n]=jiepb(A,b)运行后输出结果请注意：因为RA=RB<n ，所以此方程组有无穷多解. RA =2,RB =2,n =33.2 三角形方程组的解法及其MATLAB 程序3.2.2 解三角形方程组的MATLAB 程序 解上三角形线性方程组b AX =的MATLAB 程序function [RA,RB,n,X]=shangsan(A,b)B=[A b]; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0,disp('请注意：因为RA~=RB ，所以此方程组无解.') return endif RA==RB if RA==ndisp('请注意：因为RA=RB=n ，所以此方程组有唯一解.') X=zeros(n,1); X(n)=b(n)/A(n,n); for k=n-1:-1:1X(k)=(b(k)-sum(A(k,k+1:n)*X(k+1:n)))/A(k,k);end elsedisp('请注意：因为RA=RB<n ，所以此方程组有无穷多解.')end end例3.2.2 用解上三角形线性方程组的MATLAB 程序解方程组⎪⎪⎩⎪⎪⎨⎧==+-=-+-=++-.63,456,7472,203254434324321x x x x x x x x x x . 解 在MATLAB 工作窗口输入程序>>A=[5 -1 2 3;0 -2 7 -4;0 0 6 5;0 0 0 3]; b=[20; -7; 4;6];[RA,RB,n,X]=shangsan(A,b)运行后输出结果请注意：因为RA=RB=n ，所以此方程组有唯一解. RA = RB =4, 4, n =4,X =[2.4 -4.0 -1.0 2.0]’3.3 高斯（Gauss ）消元法和列主元消元法及其MATLAB 程序3.3.1 高斯消元法及其MATLAB 程序用高斯消元法解线性方程组b AX =的MATLAB 程序f unction [RA,RB,n,X]=gaus(A,b)B=[A b]; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0,disp('请注意：因为RA~=RB ，所以此方程组无解.') return endif RA==RB if RA==ndisp('请注意：因为RA=RB=n ，所以此方程组有唯一解.') X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1for k=p+1:nm= B(k,p)/ B(p,p);B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1); end endb=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q); endelsedisp('请注意：因为RA=RB<n ，所以此方程组有无穷多解.') end end例3.3.2 用高斯消元法和MATLAB 程序求解下面的非齐次线性方程组，并且用逆矩阵解方程组的方法验证.⎪⎪⎩⎪⎪⎨⎧-=+---=+--=+--=-+-.142,16422,0,13432143214324321x x x x x x x x x x x x x x x 解 在MATLAB 工作窗口输入程序>> A=[1 -1 1 -3; 0 -1 -1 1;2 -2 -4 6;1 -2 -4 1]; b=[1;0; -1;-1]; [RA,RB,n,X] =gaus (A,b)运行后输出结果请注意：因为RA=RB=n ，所以此方程组有唯一解. RA =4RB =4n =43.3.2 列主元消元法及其MATLAB 程序用列主元消元法解线性方程组b AX =的MATLAB 程序function [RA,RB,n,X]=liezhu(A,b)B=[A b]; n=length(b); RA=rank(A);X = 0 -0.5000 0.5000 0RB=rank(B);zhica=RB-RA; if zhica>0,disp('请注意：因为RA~=RB ，所以此方程组无解.') return endif RA==RB if RA==ndisp('请注意：因为RA=RB=n ，所以此方程组有唯一解.') X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1[Y,j]=max(abs(B(p:n,p))); C=B(p,:); B(p,:)= B(j+p-1,:); B(j+p-1,:)=C; for k=p+1:nm= B(k,p)/ B(p,p);B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1); end endb=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q); endelsedisp('请注意：因为RA=RB<n ，所以此方程组有无穷多解.') end end例3.3.3 用列主元消元法解线性方程组的MATLAB 程序解方程组⎪⎪⎩⎪⎪⎨⎧-=+---=+--=-+-=+--.142,16422,13,0432143214321432x x x x x x x x x x x x x x x . 解 在MATLAB 工作窗口输入程序>> A=[0 -1 -1 1;1 -1 1 -3;2 -2 -4 6;1 -2 -4 1]; b=[0;1;-1;-1]; [RA,RB,n,X]=liezhu(A,b)运行后输出结果请注意：因为RA=RB=n ，所以此方程组有唯一解. RA = 4，RB = 4，n = 4，X =[0 -0.5 0.5 0]’3.4 LU 分解法及其MATLAB 程序3.4.1判断矩阵LU 分解的充要条件及其MATLAB 程序 判断矩阵A 能否进行LU 分解的MATLAB 程序function hl=pdLUfj(A)[n n] =size(A); RA=rank(A); if RA~=ndisp('请注意：因为A 的n 阶行列式hl 等于零，所以A 不能进行LU 分解.A的秩RA 如下：'), RA,hl=det(A); returnendif RA==nfor p=1:n,h(p)=det(A(1:p, 1:p));, endhl=h(1:n); for i=1:nif h(1,i)==0disp('请注意：因为A 的r 阶主子式等于零，所以A 不能进行LU 分解.A 的秩RA 和各阶顺序主子式值hl 依次如下：'),hl;RA,returnend endif h(1,i)~=0disp('请注意：因为A 的各阶主子式都不等于零，所以A 能进行LU 分解.A 的秩RA 和各阶顺序主子式值hl 依次如下：')hl;RA end end例3.4.1 判断下列矩阵能否进行LU 分解，并求矩阵的秩.（1）⎪⎪⎪⎭⎫ ⎝⎛6547121321；（2）⎪⎪⎪⎭⎫ ⎝⎛654721321；（3）⎪⎪⎪⎭⎫ ⎝⎛654321321.解 （1）在MATLAB 工作窗口输入程序>> A=[1 2 3;1 12 7;4 5 6];hl=pdLUfj(A)运行后输出结果为请注意：因为A 的各阶主子式都不等于零，所以A 能进行LU 分解.A 的秩RA 和各阶顺序主子式值hl 依次如下：RA = 3， hl = 1 10 -48（2）在MATLAB 工作窗口输入程序>> A=[1 2 3;1 2 7;4 5 6];hl=pdLUfj(A)运行后输出结果为请注意：因为A 的r 阶主子式等于零，所以A 不能进行LU 分解.A 的秩RA 和各阶顺序主子式值hl 依次如下：RA = 3， hl =1 0 12（3）在MATLAB 工作窗口输入程序>> A=[1 2 3;1 2 3;4 5 6];hl=pdLUfj(A)运行后输出结果为请注意：因为A 的n 阶行列式hl 等于零，所以A 不能进行LU 分解.A 的秩RA 如下RA = 2， hl = 03.4.2 直接LU 分解法及其MATLAB 程序 将矩阵A 进行直接LU 分解的MATLAB 程序function hl=zhjLU(A)[n n] =size(A); RA=rank(A); if RA~=ndisp('请注意：因为A 的n 阶行列式hl 等于零，所以A 不能进行LU 分解.A的秩RA 如下:'), RA,hl=det(A);return endif RA==n for p=1:nh(p)=det(A(1:p, 1:p)); endhl=h(1:n); for i=1:nif h(1,i)==0disp('请注意：因为A 的r 阶主子式等于零，所以A 不能进行LU 分解.A的秩RA 和各阶顺序主子式值hl 依次如下：'), hl;RAreturn end endif h(1,i)~=0disp('请注意：因为A 的各阶主子式都不等于零，所以A 能进行LU 分解.A的秩RA 和各阶顺序主子式值hl 依次如下：')for j=1:nU(1,j)=A(1,j); endfor k=2:n for i=2:n for j=2:nL(1,1)=1;L(i,i)=1; if i>jL(1,1)=1;L(2,1)=A(2,1)/U(1,1);L(i,1)=A(i,1)/U(1,1);L(i,k)=(A(i,k)- L(i,1:k-1)*U(1:k-1,k))/U(k,k); elseU(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j); end end end endhl;RA,U,L end end例3.4.3 用矩阵进行直接LU 分解的MA TLAB 程序分解矩阵⎪⎪⎪⎪⎪⎭⎫⎝⎛=3010342110100201A . 解 在MATLAB 工作窗口输入程序>> A=[1 0 2 0;0 1 0 1;1 2 4 3;0 1 0 3]; hl=zhjLU(A)运行后输出结果请注意：因为A 的各阶主子式都不等于零，所以A 能进行LU 分解.A 的秩RA和各阶顺序主子式值hl 依次如下：RA = 4U = 1 0 2 00 1 0 10 0 2 10 0 0 2 3.4.4 判断正定对称矩阵的方法及其MATLAB 程序 判断矩阵A 是否是正定对称矩阵的MATLAB 程序function hl=zddc(A) [n n] =size(A); for p=1:nh(p)=det(A(1:p, 1:p)); endhl=h(1:n);zA=A'; for i=1:nif h(1,i)<=0disp('请注意：因为A 的各阶顺序主子式hl 不全大于零，所以A 不是正定的.A 的转置矩阵zA 和各阶顺序主子式值hl 依次如下：'), hl;zA,returnend endif h(1,i)>0disp('请注意：因为A 的各阶顺序主子式hl 都大于零，所以A 是正定的.A 的转置矩阵zA 和各阶顺序主子式值hl 依次如下：')hl;zA endL = 1 0 0 0 0 1 0 0 1 2 1 0 0 1 0 1 hl = 1 1 2 4例3.4.5 判断下列矩阵是否是正定对称矩阵：（1）⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--98754113211143214321.0；(2) ⎪⎪⎪⎪⎪⎭⎫ ⎝⎛------19631690230311211； (3) ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛----212100212100002121002121；（4）⎪⎪⎪⎭⎫⎝⎛---401061112. 解 （1）在MATLAB 工作窗口输入程序>> A=[0.1 2 3 4;-1 2 -3 4;11 21 13 41;5 7 8 9];hl=zddc (A)运行后输出结果请注意： A 不是对称矩阵请注意：因为A 的各阶顺序主子式hl 不全大于零，所以A 不是正定的.A 的转置矩阵zA 和各阶顺序主子式值hl 依次如下：zA = 1/10 -1 11 5 2 2 21 7 3 -3 13 8 4 4 41 9 hl = 1/10 11/5 -1601/10 3696/5因此,A 即不是正定矩阵，也不是对称矩阵.（2）在MATLAB 工作窗口输入程序>> A=[1 -1 2 1;-1 3 0 -3;2 0 9 -6;1 -3 -6 19],hl=zddc(A)运行后输出结果A = 1 -1 2 1 -1 3 0 -3 2 0 9 -6 1 -3 -6 19 请注意： A 是对称矩阵请注意：因为A 的各阶顺序主子式hl 都大于零，所以A 是正定的.A 的转置矩阵zA和各阶顺序主子式值hl 依次如下：zA = 1 -1 2 1 -1 3 0 -3 2 0 9 -6 1 -3 -6 19 hl = 1 2 6 24 （3）在MATLAB 工作窗口输入程序>> A=[1/sqrt(2) -1/sqrt(2) 0 0; -1/sqrt(2) 1/sqrt(2) 0 0; 00 1/sqrt(2) -1/sqrt(2); 0 0 -1/sqrt(2) 1/sqrt(2)], hl=zddc (A) 运行后输出结果A= 985/1393 -985/1393 0 0 -985/1393 985/1393 0 0 0 0 985/1393 -985/1393 0 0 -985/1393 985/1393 请注意： A 是对称矩阵请注意：因为A 的各阶顺序主子式hl 不全大于零，所以A 不是正定的.A 的转置矩阵zA 和各阶顺序主子式值hl 依次如下：zA = 985/1393 -985/1393 0 0 -985/1393 985/1393 0 0 0 0 985/1393 -985/1393 0 0 -985/1393 985/1393 hl = 985/1393 0 0 0可见，A 不是正定矩阵，是半正定矩阵；因为A = A T因此，A 是对称矩阵.（4）在MATLAB 工作窗口输入程序>> A=[-2 1 1;1 -6 0;1 0 -4];hl=zddc (A)运行后输出结果A = -2 1 11 -6 0 1 0 -4 请注意： A 是对称矩阵请注意：因为A 的各阶顺序主子式hl 不全大于零，所以A 不是正定的.A 的转置矩阵zA 和各阶顺序主子式值hl 依次如下：zA = -2 1 1 hl = -2 11 -38 1 -6 0 1 0 -4可见A 不是正定矩阵，是负定矩阵；因为A = A T因此，A 是对称矩阵.3.5 求解线性方程组的LU 方法及其MATLAB 程序3.5.1 解线性方程组的楚列斯基（Cholesky ）分解法及其MATLAB 程序例3.5.1 先将矩阵A 进行楚列斯基分解,然后解矩阵方程b AX =，并用其他方法验证.⎪⎪⎪⎪⎪⎭⎫⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛------=7531,19631690230311211b A . 解 在工作窗口输入>>A=[1 -1 2 1;-1 3 0 -3; 2 0 9 -6;1 -3 -6 19];b1=1:2:7; b=b1'; R=chol(A);C=A-R'*R,R1=inv(R);R2=R1'; x=R1*R2*b,Rx=A\b-x运行后输出方程组的解和验证结果x = Rx = 1.0e-014 * C = 1.0e-015 * -8.0000 -0.7105 0 0 0 0 0.3333 -0.0833 0 -0.4441 0 0 3.6667 0.2220 0 0 0 0 2.0000 0.1332 0 0 0 03.5.2 解线性方程组的直接LU 分解法及其MATLAB 程序例3.5.2 首先将矩阵A 直接进行LU 分解，然后解矩阵方程b AX =⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=3010342110100201A ，⎪⎪⎪⎪⎪⎭⎫⎝⎛-=5121b . 解 (1) 首先将矩阵A 直接进行LU 分解.在MATLAB 工作窗口输入程序>> A=[1 0 2 0;0 1 0 1;1 2 4 3;0 1 0 3];b=[1;2;-1;5];hl=zhjLU(A),A-L*U 运行后输出LU 分解请注意：因为A 的各阶主子式都不等于零，所以A 能进行LU 分解.A 的秩RA和各阶顺序主子式值hl 依次如下：RA = 4U = 1 0 2 00 1 0 10 0 2 10 0 0 2A 分解为一个单位下三角形矩阵L 和一个上三角形矩阵U 的积 LU A =.(2)在工作窗口输入>> U=[1 0 2 0;0 1 0 1;0 0 2 1;0 0 0 2]; L=[1 0 0 0;0 1 0 0;1L = 1 0 0 0 0 1 0 0 1 2 1 0 0 1 0 1 hl = 1 1 2 42 1 0;0 1 0 1];b=[1;2;-1;5];U1=inv(U); L1=inv(L); X=U1*L1*b,x=A\b运行后输出方程组的解X = 8.50000000000000 0.50000000000000 -3.75000000000000 1.500000000000003.5.3 解线性方程组的选主元的LU 方法及其MATLAB 程序例3.5.3 先将矩阵A 进行LU 分解，然后解矩阵方程b AX = 其中⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--=98754113211143214321.0A ，⎪⎪⎪⎪⎪⎭⎫⎝⎛-=5121b . 解 方法1 根据（3.28）式编写MATLAB 程序，然后在工作窗口输入>> A=[0.1 2 3 4;-1 2 -3 4;11 21 13 41;5 7 8 9];b=[1;2;-1;5]; [L U P]=LU(A), U1=inv(U); L1=inv(L); X=U1*L1*P*b运行后输出结果L = 1.0000 0 0 0 -0.0909 1.0000 0 0 0.0091 0.4628 1.0000 0 0.4545 -0.6512 0.2436 1.0000 U =11.0000 21.0000 13.0000 41.0000 0 3.9091 -1.8182 7.72730 0 3.7233 0.05120 0 0 -4.6171方法2 根据（3.29）式编写MATLAB 程序，然后在工作窗口输入>> A=[0.1 2 3 4;-1 2 -3 4;11 21 13 41;5 7 8 9];b=[1;2;-1;5]; [F U]=LU(A), U1=inv(U); F1=inv(F); X=U1*F1*b运行后输出结果F=0.0091 0.4628 1.0000 0 -0.0909 1.0000 0 0 1.0000 0 0 0 0.4545 -0.6512 0.2436 1.0000 X =[-1.2013 3.3677 0.0536 -1.4440]’ 用LU 分解法解线性方程组A n m ⨯b X =的MATLAB 程序function [RA,RB,n,X,Y]=LUjfcz(A,b)[n n] =size(A);B=[A b]; RA=rank(A); RB=rank(B); for p=1:nh(p)=det(A(1:p, 1:p)); endhl=h(1:n); for i=1:nif h(1,i)==0disp('请注意：因为A 的r 阶主子式等于零，所以A 不能进行LU 分解.A的秩RA 和各阶顺序主子式值hl 依次如下：')hl;RA return end endif h(1,i)~=0disp('请注意：因为A 的各阶主子式都不等于零，所以A 能进行LU 分解.A 的秩RA 和各阶顺序主子式值hl 依次如下：')P = 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 1 X =[-1.2013 3.3677 0.0536 -1.4440]’ U=11.0000 21.0000 13.0000 41.0000 0 3.9091 -1.8182 7.7273 0 0 3.7233 0.0512 0 0 0 -4.6171X=zeros(n,1); Y=zeros(n,1); C=zeros(1,n);r=1:1; for p=1:n-1[max1,j]=max(abs(A(p:n,p))); C=A(p,:); A(p,:)= A(j+P1,:); C= A(j+P1,:); g=r(p); r(p)= r(j+P1); r(j+P1)=g; for k=p+1:nH= A(k,p)/A(p,p); A(k,p) = H; A(k,p+1:n)=A(k,p+1:n)-H* A(p,p+1:n);end endY(1)=B(r(1)); for k=2:nY(k)= B(r(k))- A(k,1:k-1)* Y(1:k-1); endX(n)= Y(n)/ A(n,n); for i=n-1:-1:1X(i)= (Y(i)- A(i, i+1:n) * X (i+1:n))/ A(i,i); end end[RA,RB,n,X,Y]’;3.6 误差分析及其两种MATLAB 程序3.6.1 用MATLAB 软件作误差分析例3.6.2 解下列矩阵方程b AX =，并比较方程（1）和（2）有何区别，它们的解有何变化.其中,13/112/111/110/19/18/17/112/111/110/19/18/17/16/111/110/19/18/17/16/15/110/19/18/17/16/15/14/19/18/17/16/15/14/13/18/17/16/15/14/13/12/17/16/15/14/13/12/11)1(⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=A ;2222221⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=b ,13/112/111/110/19/18/17/112/111/110/19/18/17/16/111/110/19/18/17/16/15/110/19/18/17/16/15/14/19/18/17/16/15/14/13/18/17/16/15/14/13/12/17/16/15/14/13/12/1001.1)2(⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=A .2222221⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=b 解 (1) 矩阵方程b AX =的系数矩阵A 为7阶希尔伯特（Hilbert ）矩阵，我们可以用下列命令计算n 阶希尔伯特矩阵>>h=hilb(n) % 输出h 为n 阶Hilbert 矩阵 在MATLAB 工作窗口输入程序>> A=hilb(7);b=[1;2;2;2;2;2;2];X=A\b运行后输出b AX =的解为 X =（-35 504 -1260 -4200 20790 -27720 12012T).（2）在MATLAB 工作窗口输入程序>> B =[0.001,zeros(1,6);zeros(6,1),zeros(6,6)]; A=(B+hilb(7)); b=[1;2;2;2;2;2;2];X=A\b运行后输出方程的解为 X=（-33 465 -966 -5181 22409 -29015 12413T).在MATLAB 工作窗口输入程序>> X =[-33, 465,-966,-5181,22409,-29015,12413]';X1 =[-35,504,-1260,-4200,20790,-27720,12012]'; wu=X1'- X' 运行后输出方程（1）和（2）的解的误差为=δX 401- 1295 1619- 981 294- 39 -2(T ).方程（1）和（2）的系数矩阵的差为⎪⎪⎭⎫⎝⎛=δ⨯⨯⨯661661001.0O O O A ，常数向量相同，则b Ax =的解的差为=δX 40112951619981294392(----T ).A 的微小变化，引起X 的很大变化，即X 对A 的扰动是敏感的.3.6.2 求P 条件数和讨论b AX =解的性态的MATLAB 程序求P 条件数和讨论b AX =解的性态的MATLAB 程序function Acp=zpjxpb(A)Acw = cond (A, inf);Ac1= cond (A,1);Ac2= cond (A,2);Acf = cond (A,'fro');dA=det(A);if (Acw>50)&(dA<0.1)disp('请注意：AX=b 是病态的，A 的∞条件数,1条件数,2条件数, 弗罗贝尼乌斯条件数和A 的行列式的值依次如下：')Acp=[Acw Ac1 Ac2 Acf dA]';elsedisp(' AX=b 是良态的，A 的∞条件数,1条件数,2条件数,弗罗贝尼乌斯条件数和A 的行列式的值依次如下：')Acp=[Acw Ac1 Ac2 Acf dA]';end例3.6.3 根据定理3.10，讨论线性方程组b AX =解的性态，并且求出A 的4种条件数.其中（1）A 为7阶希尔伯特矩阵；（2）⎪⎪⎪⎪⎪⎭⎫ ⎝⎛-----=7421631472135132A . 解 （1）首先将求P 条件数和讨论b AX =解的性态的MATLAB 程序保存名为zpjxpb.m 的M 文件，然后在MATLAB 工作窗口输入程序>> Acp =zpjxpb(hilb(7)); Acp',det(hilb(7))运行后输出结果请注意：AX=b 是病态的，A 的∞条件数,1条件数,2条件数, 弗罗贝尼乌斯条件数和A 的行列式的值依次如下：ans = 1.0e+008 *9.8519 9.8519 4.7537 4.8175 0.0000ans = 4.8358e-025（2）在MATLAB 工作窗口输入程序>> A=[2 3 -1 5;3 1 2 -7;4 1 -3 6;1 -2 4 -7];Acp=zpjxpb(A); Acp' 运行后输出结果 AX=b 是良态的，A 的∞条件数,1条件数,2条件数, 弗罗贝尼乌斯条件数和A 的行列式的值依次如下：ans =14.1713 19.4954 8.2085 11.4203 327.00003.6.3 用P 范数讨论b AX =解和A 的性态的MATLAB 程序用P 范数讨论b AX =解和A 的性态的MATLAB 程序function Acp=zpjwc(A,jA,b,jb,P)Acp = cond (A,P);dA=det(A); X=A\b;dertaA=A-jA;PndA=norm(dertaA, P);dertab=b-jb;Pndb=norm(dertab, P);if Pndb>0jX=A\jb; Pnb= norm(b, P);PnjX = norm(jX,P); dertaX=X-jX;PnjdX= norm(dertaX, P);jxX= PnjdX/PnjX; PnjX =norm(jX,P);PnX = norm(X,P); jxX= PnjdX/PnjX; xX= PnjdX/PnX;Pndb=norm(dertab,P);xAb=Pndb/Pnb;Pnbj=norm(jb,P); xAbj=Pndb/Pnbj;Xgxx= Acp*xAb;endif PndA>0jX=jA\b; dertaX=X-jX;PnX = norm(X,P);PnjdX= norm(dertaX, P);PnjX = norm(jX,P); jxX= PnjdX/PnjX;xX= PnjdX/PnX;PnjA=norm(jA,P); PnA=norm(A,P);PndA=norm(dertaA,P);xAbj= PndA/PnjA;xAb= PndA/PnA;Xgxx= Acp*xAb;endif (Acp >50)&(dA<0.1)disp('请注意：AX=b 是病态的，A 的P 条件数Acp,A 的行列式值dA ，解X ，近似解jX ，解的相对误差jxX ，解的相对误差估计值Xgxx ，b 或A 的相对误差xAb 依次如下：')Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj'elsedisp('请注意： AX=b 是良态的，A 的P 条件数Acp,A 的行列式值dA ，解X ，近似解jX ，解的相对误差jxX ，解的相对误差估计值Xgxx ，b 或A 的相对误差xAb 依次如下：')Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj'end例3.6.4 根据定理3.10，讨论线性方程组b AX =解的性态，并利用（3.32）式讨论当A 的每个元都取4位有效数字时，其解的相对误差.其中A 为7阶希尔伯特矩阵，()22224311=b T .解 （1）取∞范数和∞条件数，线性方程组b AX =的b 不变时，取∞范数和∞条件数，系数矩阵A 为7阶希尔伯特矩阵，A 中的每个元素取4位有效数字.用P 范数讨论b AX =解和A 的性态的MATLAB 程序保存名为zpjwc.m 的文件,然后在工作窗口输入MATLAB 程序>> jA =[1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.14290.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.12500.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.11110.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.10000.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.09090.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.08330.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769]; A=hilb(7);b=[1;1/3;4;2;2;2;2];jb=[1;0.3333;4;2;2;2;2]; Acp=zpjwc(A,jA,b,jb,inf)运行后输出结果请注意：AX=b 是病态的，A 的P 条件数Acp,A 的行列式值dA ，解X ，近似解jX ，解的相对误差jxX ，解的相对误差估计值Xgxx ，b 或A 的相对误差xAb 依次如下：Acp = dA =9.8519e+008 4.8358e-025ans =1.0e+007 *0.0020 -0.0697 0.6243 -2.3054 4.0677 -3.41231.0943ans =1.0e+004 *0.0349 -0.4807 2.1126 -5.1087 7.6557 -6.32392.1112xX = jxX = Xgxx =0.9981 530.3248 4.9291e+004 xAb = xAbj =5.0032e-005 5.0031e-005由此可见，因为∞条件数(Cond ≈∞)A 985 194 889.719 848 31>>，所以此方程组为病态的b AX =的解T 932)94210320,12334-790,676380,402.305-300,2436 256,697- 565,19( =X ， b X jA =)(的解为T )11221239,63-,55767087,51-126,21 079,48- ,493( =jX 解的相对误差120.998≈∞∞X Xδ,291.2749≤+∞∞X X X δδ , 530.32≈+∞∞X X X δδ,,1025.0035-∞∞⨯≈A Aδ即相对误差放大了约985 194 889.72倍.(2) 如果取2范数和2条件数计算，在MATLAB 工作窗口输入程序>> Acp =zpjwc(A,jA,b,jb,2)运行后输出结果请注意：AX=b 是病态的，A 的P 条件数Acp,A 的行列式值dA ，解X ，近似解jX ，解的相对误差jxX ，解的相对误差估计值Xgxx ，b 或A 的相对误差xAb 依次如下：Acp = dA =4.7537e+008 4.8358e-025ans =1.0e+007 *0.0020 -0.0697 0.6243 -2.3054 4.0677 -3.41231.0943ans =1.0e+004 *0.0349 -0.4807 2.1126 -5.1087 7.6557 -6.32392.1112xX = jxX = Xgxx=0.9981 511.0640 2.9951e+004xAb = xAbj =6.3006e-005 6.3005e-005因为2条件数Cond ≈2)(A 475 367 356.591>>，所以此方程组为病态的.解的相对误差10.998 22≈X X δ,,105300.6522-⨯≈A A δ,06.51122≈+X X X δδ.85.9502922≤+A A A δδ 即相对误差放大了约475 367 356.59倍.。