解线性方程组直接方法matlab用法
Matlab求解线性方程组、非线性方程组
求解线性方程组solve，linsolve例：A=[5 0 4 2;1 -1 2 1;4 1 2 0;1 1 1 1];%矩阵的行之间用分号隔开，元素之间用逗号或空格B=[3;1;1;0]X=zeros(4,1);%建立一个4元列向量X=linsolve(A,B)diff（fun，var，n）：对表达式fun中的变量var求n阶导数。
例如：F=sym（'u(x,y)*v(x,y)'）; %sym（）用来定义一个符号表达式diff(F); %matlab区分大小写pretty(ans) %pretty（）：用习惯书写方式显示变量；ans是答案表达式非线性方程求解fsolve(fun,x0,options)其中fun为待解方程或方程组的文件名；x0位求解方程的初始向量或矩阵；option为设置命令参数建立文件fun.m：function y=fun(x)y=[x(1)-0.5*sin(x(1))-0.3*cos(x(2)), ...x(2) - 0.5*cos(x(1))+0.3*sin(x(2))];>>clear;x0=[0.1,0.1];fsolve(fun,x0,optimset('fsolve'))注：...为续行符m文件必须以function为文件头，调用符为；文件名必须与定义的函数名相同；fsolve（）主要求解复杂非线性方程和方程组，求解过程是一个逼近过程。
Matlab求解线性方程组AX=B或XA=B在MATLAB中，求解线性方程组时，主要采用前面章节介绍的除法运算符“/”和“\”。
如：X=A\B表示求矩阵方程AX＝B的解；X＝B/A表示矩阵方程XA=B的解。
对方程组X＝A\B，要求A和B用相同的行数，X和B有相同的列数，它的行数等于矩阵A的列数，方程X＝B/A同理。
如果矩阵A不是方阵，其维数是m×n，则有：m＝n 恰定方程，求解精确解；m>n 超定方程，寻求最小二乘解；m<n 不定方程，寻求基本解，其中至多有m个非零元素。
matlab在科学计算中的应用 线性方程组的直接求解-分析方法-文档资料
>> jacobi(a,b,[0;0;0])
n= 11
ans =
0.9958
0.9579 0.7916
4.3.2Байду номын сангаасGauss-Seidel迭代法
Jacobi迭代：x(k1) Bx(k) f
x(k+1)
1
x(k+1) 2
b21x1(k )
b12 x2(k )
b1n xn(k ) f1 b2n xn(k ) f2 ,
16 4 8
例：进行Cholesky分解。
A
4
5
4
8 4 22
>> A=[16 4 8; 4 5 -4; 8 -4 22];
>> D=chol(A)
D=
412
0 2 -3
003
●利用矩阵的LU、QR和cholesy分解求方程组的解
（1）LU分解：
A*X=b
变成 L*U*X=b
a1n xn b1 a2n xn b2
ann xn bn
其中det( A) det(aij ) 0,不妨设aii 0(i 1, 2,..., n)
x1
x2
b21x1
b12 x2
b1n xn f1 b2n xn f2 ,
迭代公式为： X(k+1)= (D-wL)-1((1-w)D+wU)x(k) + w(D-wL)-1 b 其中：w最佳值在[1, 2）之间，不易计
算得到，因此 w通常有经验给出。
function y=sor(a,b,w,x0) D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1); M=(D-w*L)\((1-w)*D+w*U); f=(D-w*L)\b*w; y=M*x0+f; n=1; while norm(y-x0)>=1.0e-6
MatLab解线性方程组
MatLab解线性方程组MatLab解线性方程组当齐次线性方程AX=0,rank(A)=r<n时,该方程有无穷多个解,怎样用matlab求它的一个基本解呢?< p="">用matlab 中的命令 x=null(A, r )即可.其中:r=rank(A)A=[ 1 1 1 1 -3 -1 11 0 0 0 1 1 0-2 0 0 -1 0 -1 -2]用matlab 求解程序为:A=[1 1 1 1 -3 -1 1;1 0 0 0 1 1 0;-2 0 0 -1 0 -1 -2];r=rank(A);y=null(A, r )得到解为:y=[ 0 -1 -1 0-1 2 1 11 0 0 00 2 1 -20 1 0 00 0 1 00 0 0 1]其列向量为Ay=0的一个基本解一：基本概念1．N级行列式A：|A|等于所有取自不同行不同列的n个元素的积的代数和。
2．矩阵B：矩阵的概念是很直观的，可以说是一张表。
3．线性无关：一向量组(a ,a ,…. a )不线性相关，即没有不全为零的数k ,k ,……kn 使得：k1* a +k2* a +…..+kn*an=04. 秩：向量组的极在线性无关组所含向量的个数称为这个向量组的秩。
5．矩阵B的秩：行秩，指矩阵的行向量组的秩；列秩类似。
记：R(B)6．一般线性方程组是指形式： (1)其中x1,x2,…….xn为n个未知数，s为方程个数。
记：A*X=b7．性方程组的增广矩阵： =8. A*X=0 (2)二：基本理论三种基本变换：1，用一非零的数乘某一方程；2，把一个方程的倍数加到另一个方程；3互换两个方程的位置。
以上称初等变换。
消元法（理论上分析解的情况，一切矩阵计算的基础）首先用初等变换化线性方程组为阶梯形方程组，把最后的一些恒等式”0=0”（如果出现的话）去掉，1：如果剩下的方程当中最后的一个等式是零等于一非零数，那么方程组无解；否则有解，在有解的情况下，2：如果阶梯形方程组中方程的个数r等于未知量的个数，那么方程组有唯一的解，3：如果阶梯形方程组中方程的个数r小于是未知量的个数，那么方程组就有无穷个解。
MATLAB计算方法3解线性方程组计算解法ppt课件
1 c1
y1
2 c2
y2
n1
cn1
yn1
n yn
回代过程为 xi
xn yn / n
x2
d2
给定
an1
bn1
cn1
xn1
d
n1
an bn xn dn
（三对角方程组）
且按行对角占优：b 1 c1 ， bi ai ci (aici 0)， bn an （i 2, , n 1）
36
利用Gauss消元法得到同解的三角方程为
L ( LK L1 )1 ，则
A LU （下三角 上三角）
（三角因子分解） 18
定义3.1 A LU 叫 A 的三角（因子）分解，其中 L是
下三角, U是上三角。
定义3.2 若 L为单位下三角阵（对角元全为1），
U 为上三角阵，则称 A LU 为Doolittle分解;
若L是L 下三角，U 是单位上三角，则称 A LU
设方程组Ax b, 对其增广矩阵作LU分解时,为了避免
k 1
用小ukk作除数,引进量si aik liquqk (i k, k 1, , n) q1
于 是ukk
sk ,比较 |
si
| 的大小, 取 max kin
|
si
| (|
st
|)为ukk ,
并交换矩阵的第t行与第k行, 且元素的足码也相应改变。
l n1
l
n2
l
nn
30
由矩阵乘法
(1)
1)
l2 11
a11
l11
a11
(取正)
2) L第1行 LT第j列 (j 2, ,n)
l11l j1
matlab解线性方程组
设有n个变量，m个方程，方程组的系数矩阵为A，常数项列向量为b，则A为m×n矩阵，b为m×l矩阵，方程组可写为Ax=b其中x为n个变量构成的列向量，若rank（A）=m，且m=n，则方程有唯一解，称为恰定方程组；设B=（A|b）为增广矩阵，且若rank （A）≠rank（B），则方程组无解，称为超定方程组；rank（A）=rank （B）=r<m，则方程有无穷多组解，称为欠定方程组。
若b中元素全为0，称方程组为齐次线性方程组。
齐次线性方程组的求解对于齐次线性方程组Ax=0，由线性代数知识知，它至少有一零解，若rank（A）<m，则它有无穷多组解。
在MATLAB中有一个非常有用的函数B=null（A），它返回了矩阵A的零空间的标准正交基组成的矩阵B，使得BTxB=E，B的列数等于矩阵A的零空间的维数，即B 的列向量构成了线性方程组的一组标准正交基础解系。
例：解方程组：在MATLAB中输入：因而，原方程组的通解为其中k1与k2为任意常数。
例5.2.2 求方阵A的含有最多零元素个数的解。
在例5.2.1求解后，运行：则有故方程组的通解为恰定方程组的求解恰定方程组Ax=b的求解比较简单。
一般可用两种方法：一种是利用逆矩阵求解：x=inv（A）b；另一种是用除法求解x=A\b。
两种方法的异同点是：算法上都采用Guass消去法，但用除法求解时，无需求A的逆，这样可以很好地保证求解时的计算精度，还能节省大量的计算时间。
当然也可以用Cramer法则求解方程组。
比较多种形式求解恰定方程组：结果为：由结果可知，用除法耗时最少，大约只有逆矩阵法的1/3，而Cramer 法则约是它的85倍；而从误差分析来看，仍然是除法最精确。
可见，用除法求解恰定方程组，既省时又精确。
值得注意的是，利用逆矩阵inv（A）求解时，若A不是方阵，则要用广义逆pinv（A）来求。
超定方程组的求解由于超定方程组无解，而在实际应用中，能求得其最小二乘解也是有意义的，这时，方程组的求解仍可用除法和广义逆矩阵法，不过这样求得的解不会满足Ax=b，而是其最小二乘意义下的解，即解x=inv （ATA）ATb。
用matlab解线性方程组
用matlab解线性方程组2008-04-12 17:00一。
高斯消去法1.顺序高斯消去法直接编写命令文件a=[]d=[]'[n,n]=size(a);c=n+1a(:,c)=d;for k=1:n-1a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); %消去endx=[0 0 0 0]' %回带x(n)=a(n,c)/a(n,n);for g=n-1:-1:1x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g)end2.列主高斯消去法*由于“[r,m]=max(abs(a(k:n,k)))”返回的行是“k:n,k”内的第几行，所以要通过修正来把m 改成真正的行的值。
该程序只是演示程序，真正机器计算不需要算主元素所在列以下各行应为零的值。
直接编写命令文件a=[]d=[] '[n,n]=size(a);c=n+1a(:,c)=d; %(增广)for k=1:n-1[r,m]=max(abs(a(k:n,k))); %选主m=m+k-1; %(修正操作行的值)if(a(m,k)~=0)if(m~=k)a([k m],:)=a([m k],:); %换行enda(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); %消去endendx=[0 0 0 0]' %回带x(n)=a(n,c)/a(n,n);for g=n-1:-1:1x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g)end3．分别用顺序高斯消去法和列主高斯消去法解方程组a*x=d，并比较结果a=[0 1 2 3;9 11 23 34;62.5 23.4 15.5 17.2;192.01 124 25.1 59.3] d=[1;1;1;1]顺序高斯消去法：提示“Warning: Divide by zero.” x =NaN NaN NaN NaN 列主高斯消去法：x =-1.2460 2.8796 5.5206 -4.3069由此可见列主高斯消去法可以解决顺序高斯消去法所不能解决的问题。
matlab中快速求解xa=b的方法
matlab中快速求解xa=b的方法在Matlab中，要快速求解线性方程组xa=b，可以使用以下几种方法：1. 直接求解法（\）：直接使用斜杠操作符（\）可以求解线性方程组。
例如，对于方程组xa=b，可以直接使用x = A\b来解决，其中A是系数矩阵，b是常数向量。
这种方法使用了高效的LU分解算法，并且能够自动适应方程组的类型（如稀疏矩阵或密集矩阵），因此是一种快速求解线性方程组的常用方法。
2. QR分解法：QR分解是一种将矩阵分解为正交矩阵和上三角矩阵的方法。
在Matlab中，可以使用qr函数对系数矩阵进行QR分解，然后使用这个分解求解线性方程组。
具体而言，可以使用[q,r] = qr(A)将系数矩阵A分解为正交矩阵q和上三角矩阵r，然后使用x = r\(q'*b)求解方程组。
这种方法通常适用于方程组的系数矩阵具有较大的条件数或者方程组数目较多的情况。
3. Cholesky分解法：如果线性方程组的系数矩阵是对称正定的，那么可以使用Cholesky分解来求解方程组。
在Matlab中，可以使用chol函数对系数矩阵进行Cholesky分解，然后使用这个分解求解线性方程组。
具体而言，可以使用R = chol(A)将系数矩阵A分解为上三角矩阵R，然后使用x = R'\(R\b)求解方程组。
Cholesky分解法通常适用于系数矩阵具有良好的性质（如对称正定）的情况。
4. 迭代法：如果线性方程组的系数矩阵是稀疏的，那么可以使用迭代法来求解方程组。
迭代法的基本思想是通过迭代改进解的逼近值。
在Matlab中，可以使用pcg函数（预处理共轭梯度法）或者bicg函数（双共轭梯度法）来求解稀疏线性方程组。
这些函数需要提供一个预处理矩阵，用于加速迭代过程。
预处理矩阵可以根据具体问题进行选择，常见的预处理方法包括不完全LU分解（ilu）和代数多重网格（amg）等。
通过使用上述方法，可以在Matlab中快速求解线性方程组xa=b。
matlab实验5 线性方程组的解法
马千里 热动71 970669实验5 线性方程组的解法实验目的1. 用MA TLAB 软件掌握线性方程组的解法，对迭代法的收敛性和解的稳定性作初步分析。
2. 通过实例练习用线性方程组求解实际问题。
实验内容预备：编写雅可比迭代和高斯—塞得尔迭代的程序 雅可比迭代的程序：function xxx=ykb(A,b) D=diag(diag(A)); U=-triu(A,1); L=-tril(A,-1); B1=inv(D)*(L+U); f=inv(D)*b;x(1:length(b),1)=0; xx=x+1;while norm(x-xx)>1e-6 xx=x;x=B1*x+f; end xxx=x;1. 用MA TLAB 软件的“lu ”(LU 分解)，“ \ ”，以及雅可比迭代和高斯—塞德尔迭代解方程组Ax=b (A 如下，b 任意，比较分析其结果包括迭代法收敛或不收敛的原因）。
a. A=[1,2,-1;1,1,1;2,2,1]; 设b=[1; 2; 1]在MA TLAB 下运行：》A=[1,2,-1;1,1,1;2,2,1];b=[1; 2; 1]; 用LU 分解方法： 》[L U p]=lu(A)L =1.0000 0 00.5000 1.0000 00.5000 0 1.0000》x=(L*U\p*b)’ x =-9 8 3 用除号“\”： 》x=(A\b)’ x =-9 8 3 用雅可比迭代法： 》x=ykb(A,b) x =高斯—塞得尔迭代的程序：function xxx=guass(A,b) D=diag(diag(A)); U=-triu(A,1); L=-tril(A,-1); B1=inv(D-L)*U; f=inv(D-L)*b;x(1:length(b),1)=0; xx=x+1;while norm(x-xx)>1e-8 xx=x;x=B1*x+f; end xxx=x;U = 2.0000 2.0000 1.0000 0 1.0000 -2.5000 0 0 0.5000 p = 0 0 1 1 0 0 0 1 0-9 8 3迭代公式中B1 = f=0 -2 2 1-1 0 -1 2-2 -2 0 1 因为ρ(B1)=0.126<1所以是收敛的用高斯—塞德尔迭代法：因为迭代公式中B1= f=0 -2 2 10 2 -3 10 0 2 3因为ρ(B1)=2>1所以迭代是不收敛的b.A为n阶Hilbert矩阵，n=3~10。
matlab线性方程组求解
0.9739 -0.0047 1.0010
n= 5 Jacobi 迭代法： function [x,n]=jacobi(A,b,x0,eps,varargin) if nargin==3 eps= 1.0e-6; M = 200; elseif nargin<3 error return elseif nargin ==5 M = varargin{1}; end D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); B=D/(L+U); f=D/b; x=B*x0+f; n=1; % 迭代次数 % 求 A 的对角矩阵 % 求 A 的上三角阵
n= 5 Gauss-Seidel 迭代法： function [x,n]=gauseidel(A,b,x0,eps,M) if nargin==3 eps= 1.0e-6; M = 200; elseif nargin == 4 M = 200; elseif nargin<3 error return; end D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); G=(D-L)/U; f=(D-L)/b; % 求 A 的对角矩阵 % 求 A 的上三角阵 % 求 A 的下三角阵
批注本地保存成功开通会员云端永久保存去开通
线性方程组求解 1. 直接法 Gauss 消元法： function x=DelGauss(a,b) % Gauss 消去法 [n,m]=size(a); nb=length(b); det=1;% 存储行列式值 x=zeros(n,1); for k=1:n-1 for i=k+1:n if a(k,k)==0 return end m=a(i,k)/a(k,k); for j=k+1:n a(i,j)=a(i,j)-m*a(k,j); end b(i)=b(i)-m*b(k); end det=det*a(k,k); end det=det*a(n,n);
matlab解方程组方法
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函数求解非线性方程的根。
MATLAB编程与工程应用——第7章 MATLAB解方程与函数极值
MATLAB解方程与函数极值
7.1 线性方程组求解
二、迭代解法
迭代解法非常适合求解大型系数矩阵的方程组。在数值 分析中，迭代解法主要包括 Jacobi迭代法、GaussSerdel迭代法、超松弛迭代法和两步迭代法。 Jacobi迭代法 1．Jacobi迭代法 对于线性方程组Ax=b，如果A为非奇异方阵，即aii≠0， 则A可分解为A=D-L-U，其中D为对角阵，其元素为A的对 角元素，L与U为A的下三角阵和上三角阵，于是Ax=b化为： x=D-1(L+U)x+D-1b 与之对应的迭代公式为： x(k+1)=D-1(L+U)x(k)+D-1b 利用jacobi jacobi迭代法求方程的解 例7.4 利用jacobi迭代法求方程的解 exp7_4.m jacobi.m
如s=2时
j =1 j =1
i 1
i 1
1 ym +1 = ym + h[ f ( xm , ym ) + f ( xm + h, ym + h)] 2
matlab求线性方程组的解
matlab求线性方程组的解求解线性方程分为两种方法–直接法和迭代法常见的方法一共有8种直接法Gauss消去法Cholesky分解法迭代法Jacobi迭代法Gauss-Seidel迭代法超松弛迭代法共轭梯度法Bicg迭代法Bicgstab迭代法这里我从计算代码的角度来解释一下，代码按以下顺序给出。
把方程组直接带入已知条件，就可以得到答案。
适用条件Gauss消去法：求解中小规模线性方程（阶数不过1000），一般用于求系数矩阵稠密而且没有任何特殊结构的线性方程组Cholesky分解法：对称正定方程优先使用，系数矩阵A是n 阶对称正定矩阵Jacobi迭代法非奇异线性方程组，分量的计算顺序没有关系Gauss-Seidel迭代法与Jacobi迭代法相似，但计算的分量不能改变超松弛迭代法Jacobi迭代法和Gauss-Seidel迭代法的加速版，由Gauss-Seidel迭代法改进而来，速度较快共轭梯度法需要确定松弛参数w,只有系数矩阵具有较好的性质时才可以找到最佳松弛因子。
但好处是不用确定任何参数，他是对称正定线性方程组的方法也是求解大型稀疏线性方程组最热门的方法Bicg迭代法本质是用双共轭梯度求解线性方程组的方法，对求解的方程没有正定性要求Bicgstab迭代法本质是用稳定双共轭梯度求解线性方程组的方法，对求解的方程没有正定性要求Gauss消去法第一、二个函数ltri、utri是一定要掌握的，后面的几乎每个函数都要用到ltri简单来说，当Ly=bb,L(非奇异下三角矩阵）已知求yfunction y =ltri(L,b)n=size(b,1);y=zeros(n,1);for j =1:n-1y(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-y(j)*L(j+1:n,j); endy(n)=b(n)/L(n,n);utri简单来说，当Ux=yy,U(非奇异上三角矩阵）已知求xfunction x =utri(U,y)n=size(y,1);x=zeros(n,1);for j = n:-1:2x(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-x(j)*U(1:j-1,j);endx(1)=y(1)/U(1,1);gauss算法，计算时粘贴过去就好function[L,U]=gauss(A)n=size(A,1);for k =1:n-1A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k +1:n)-A(k+1:n,k)*A(k,k+1:n);endL=tril(A,-1)+eye(n);U=triu(A);使用例子已经知道一个线性方程组，这里我就不写出数学形式了，A是系数矩阵，直接把上面写好的函数复制过来在运算就可以。
matlab计算方程组
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迭代法的函数来求解方程组。
MATLAB实验一 解线性方程组的直接法
输出 Ax b 中系数 A LU 分解的矩阵 L 和 U ，解向量 x 和 det(A) ;用列主元法 的行交换次序解向量 x 和求 det(A) ；比较两种方法所得结果。
2、用列主高斯消元法解线性方程组 Ax b 。
3.01 6.03 1.99 x1 1 4.16 1.23 x 2 1 （1） 、 1.27 0.987 4.81 9.34 x 1 3 3.00 6.03 1.99 x1 1 4.16 1.23 x 2 1 （2） 、 1.27 0.990 4.81 9.34 x 1 3
index = 1 3、在 MATLAB 窗口：
A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; b=[32 23 33 31]'; x=A\b b1=[32.1 22.9 33.1 30.9]'; x1=A\b1 A1=[10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98]; x2=A1\b delta_b=norm(b-b1)/norm(b) delta_A=norm(A-A1)/norm(A) delta_x1=norm(x-x1)/norm(x) delta_x2=norm(x-x2)/norm(x)
二. 实验要求 1、按照题目要求完成实验内容； 2、写出相应的 Matlab 程序； 3、给出实验结果(可以用表格展示实验结果)； 4、分析和讨论实验结果并提出可能的优化实验。 5、写出实验报告。 三. 实验步骤 1、用 LU 分解及列主元高斯消去法解线性方程组
7 10 3 2.099999 a) 5 1 2 1 1 x1 8 6 2 x 2 5.900001 ， 5 1 x3 5 0 2 1 x 4 0
MATLAB线性方程组求解方法
——GW318 物联网实验室学术活动
MATLAB 线性方程组求解方法
1．线性方程组的问题 在工程计算中，一个重要的问题是线性方程组的求解。在矩阵表示法中，上述问题可以 表述为给定两个矩阵 A 和 B，是否存在惟一的解 X 使得 AX=B 或 XA=B。 例如，求解方程 3x=6 就可以将矩阵 A 和 B 看成是标量的一种情况，最后得出该方 程的 解为 x=6/3=2。 尽管在标准的数学中并没有矩阵除法的概念，但 MATLAB 7.0 采用了与解标量方程中 类 似的约定，用除号来表示求解线性方程的解。MATLAB 7.0 采用第 2 章介绍过的运算 符斜杠 ’/’和运算符反斜杠’\ ’来表示求线性方程的解，其具体含义如下： • X=A\B 表示求矩阵方程 AX=B 的解； • X=B/A 表示求矩阵方程 XA=B 的解。 对于 X=A\B，要求矩阵 A 和矩阵 B 有相同的行数，X 和 B 有相同的列数，X 的行 数等于 矩阵 A 的列数。X=B\A 行和列的性质则与之相反。 在实际情况中，形式 AX=B 的线性方程组比形式 XA=B 的线性方程组要常见得多。 因此 反斜杠’\’用得更多。本小节的内容也主要针对反斜杠’\’除法进行介绍。斜杠’/’ 除法的性质可以 由恒等变换式得到(B/A)’=(A’\B’)。 系数矩阵 A 不一定要求是方阵，矩阵 A 可以是 m×n 的矩阵，有如下 3 种情况： • m=n 恰定方程组，MATLAB7.0 会寻求精确解； • m>n 超定方程组，MATLAB7.0 会寻求最小二乘解； • m<n 欠定方程组，MATLAB7.0 会寻求基本解，该解最多有 m 个非零元素。 值得注意的是用 MATLAB 7.0 求解这种问题时，并不采用计算矩阵的逆的方法。针对 不 同的情况，MATLAB7.0 会采用不同的算法来解线性方程组。 2．线性方程组的一般解 线性方程组 AX=B 的一般解给出了满足 AX=B 的所有解。线性方程组的一般解可以通 过 下面的步骤得到。 • 解相应的齐次方程组 AX=0， 求得基础解。 可以使用函数 null()来得到基础解。 语 句 null(A) 返回齐次方程组 AX=0 的一个基础解，其他基础解与 null(A)是线性关系。 • 求非齐次线性方程组 AX=B，得到一个特殊解。 • 非齐次线性方程组 AX=B 的一般解等于基础解的线性组合加上特殊解。 在后面的章节将介绍求非齐次线性方程组 AX=B 特殊解的方法。 3．恰定方程组的求解
matlab中解方程组
MATLAB中解方程组1. 引言在科学计算和工程领域，解方程组是一个常见的任务。
MATLAB作为一种强大的数值计算软件，提供了多种方法来解决方程组问题。
本文将介绍MATLAB中解方程组的基本方法和技巧。
2. 方程组的表示在MATLAB中，我们可以使用矩阵和向量的形式表示线性方程组。
例如，考虑以下线性方程组：2x + 3y = 74x - y = -1可以将其表示为矩阵和向量的形式：A = [2, 3; 4, -1]B = [7; -1]其中A是系数矩阵，B是常数向量。
3. 使用反斜杠运算符求解方程组MATLAB提供了一个简单而强大的运算符\来求解线性方程组。
例如，我们可以使用以下代码求解上述方程组：A = [2, 3; 4, -1];B = [7; -1];X = A \ B;运行以上代码后，变量X将包含方程组的解。
通过命令disp(X)可以打印出结果。
4. 解非线性方程组除了线性方程组外，MATLAB还可以用于求解非线性方程组。
非线性方程组的求解更加复杂，通常需要使用数值方法来逼近解。
MATLAB提供了多种函数和工具箱来求解非线性方程组。
其中最常用的是fsolve函数，它可以通过迭代方法求解非线性方程组。
例如，考虑以下非线性方程组：x^2 + y^2 = 1x + y = 1我们可以使用fsolve函数求解该方程组：fun = @(x) [x(1)^2 + x(2)^2 - 1; x(1) + x(2) - 1];x0 = [0; 0];options = optimoptions('fsolve', 'Display', 'iter');[x, fval] = fsolve(fun, x0, options);在以上代码中，fun是一个匿名函数，表示要求解的非线性方程组。
x0是初始猜测值，options是优化选项。
运行以上代码后，变量x将包含方程组的解，fval将包含目标函数的值。
10.matlab求解线性方程组
10. 如何求解线性方程组？
当一个方程组Ax=b 有唯一解时可用指令A\b 直接求解，其中，A 是线性方程组的系数矩阵，b是方程组的右端向量。
例如给定线性方程组
[121−2 253−2
−2−235 1323][
x1
x2
x3
x4
]=[
4
7
−1
]
可以先输入系数矩阵和右端向量，然后直接求解，在MATLAB中键入
A=[1 2 1 -2； 2 5 3 -2；-2 -2 3 5； 1 3 2 3]；
b=[4；7；-1；0]；
x=A\b
计算机执行后，将显示数据结果
x= 2 -1 2 -1
由此得知方程组的解为
x1 = 2，x2 = -1，x3 = 2，x4 = -1
上面这一种方法非常适用于方程组有唯一解的情形，在处理实际问题时，有些方程组有无穷多组解，这时可用另一条指令rref([A b]) 化简方程组的增广矩阵，然后利用线性代数的方法得出方程组的通解。
当方程组有唯一解时，仍然可以用这一指令，例如求解上面例子
在MATLAB中键入
A=[1 2 1 -2； 2 5 3 -2；-2 -2 3 5； 1 3 2 3]；
b=[4；7；-1；0]；
rref([A b])
计算机执行后，屏幕将显示数据结果
这一矩阵与增广矩阵等价(4 个非零行5 列)，由此写出与原方程等价的方程组，即
x1 = 2，x2 = -1，x3 = 2，x4 = -1
如果用rref([A b])命令得到最后的矩阵中非零行数小于列数减1，则可求出线性方程组的通解。
有关线性方程组的通解和基础解系的概念请参考线性代数教材。
MATLAB解方程组(线性与非线性方程组)
例7-9 求下列非线性方程组在(0.5,0.5) 附近的数值解。 (1) 建立函数文件myfun.m。 function q=myfun(p) x=p(1); y=p(2); q(1)=x-0.6*sin(x)-0.3*cos(y); q(2)=y-0.6*cos(x)+0.3*sin(y); (2) 在给定的初值x0=0.5,y0=0.5下，调用fsolve函数求方程的根。 x=fsolve('myfun',[0.5,0.5]',optimset('Display','off')) x= 0.6354 0.3734
2．Gauss-Serdel迭代法 在Jacobi迭代过程中，计算时，已经得到，不必再用，即原来的迭代
公式Dx(k+1)=(L+U)x(k)+b可以改进为Dx(k+1)=Lx(k+1)+Ux(k)+b， 于是得到：
x(k+1)=(D-L)-1Ux(k)+(D-L)-1b 该式即为Gauss-Serdel迭代公式。和Jacobi迭代相比，Gauss-Serdel
7.1.2 迭代解法 迭代解法非常适合求解大型系数矩阵的方程组。在数值分析中，迭代
解法主要包括 Jacobi迭代法、Gauss-Serdel迭代法、超松弛迭代法 和两步迭代法。
1．Jacobi迭代法 对于线性方程组Ax=b，如果A为非奇异方阵，即aii≠0(i=1,2,…,n)，则
可将A分解为A=D-L-U，其中D为对角阵，其元素为A的对角元素， L与U为A的下三角阵和上三角阵，于是Ax=b化为： x=D-1(L+U)x+D-1b 与之对应的迭代公式为：
(2) QR分解
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
2.1 方程组地逆矩阵解法及其MATLAB 程序2.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例 2.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 =32.2 三角形方程组地解法及其MATLAB程序2.2.2 解三角形方程组地MATLAB程序解上三角形线性方程组bAX 地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例2.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 03];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]’2.3 高斯（Gauss ）消元法和列主元消元法及其MATLAB 程序2.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例2.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 =42.3.2 列主元消元法及其MATLAB 程序用列主元消元法解线性方程组b AX =地MATLAB 程序function [RA,RB,n,X]=liezhu(A,b)X = 0 -0.5000 0.5000 0B=[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-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例2.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]’2.4 LU 分解法及其MATLAB 程序2.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); return endif 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,return endif h(1,i)~=0disp('请注意：因为A 地各阶主子式都不等于零，所以A 能进行LU 分解.A 地秩RA 和各阶顺序主子式值hl 依次如下：')hl;RA end end例2.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 = 02.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;RA return end endif h(1,i)~=0disp('请注意：因为A 地各阶主子式都不等于零，所以A 能进行LU 分解.A地秩RA 和各阶顺序主子式值hl 依次如下：')for j=1:nU(1,j)=A(1,j);for 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例2.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 2.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,return end endif h(1,i)>0disp('请注意：因为A 地各阶顺序主子式hl 都大于零，所以A 是正定地.A 地转置矩阵zA 和各阶顺序主子式值hl 依次如下：') hl;zA end例2.4.5 判断下列矩阵是否是正定对称矩阵：L = 1 0 0 0 0 1 0 0 1 2 1 0 0 1 0 1 hl = 1 1 2 4（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 是对称矩阵.2.5 求解线性方程组地LU 方法及其MATLAB 程序2.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 00.3333 -0.0833 0 -0.4441 0 03.6667 0.2220 0 0 0 02.0000 0.1332 0 0 0 02.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)在工作窗口输入L = 1 0 0 0 0 1 0 0 1 2 1 0 0 1 0 1 hl = 1 1 2 4>> 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;12 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.500000000000002.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);endendY(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);endend[RA,RB,n,X,Y]’;2.6 误差分析及其两种MATLAB 程序2.6.1 用MATLAB 软件作误差分析例2.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 地扰动是敏感地.2.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例2.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.0000 2.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例2.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+004xAb = 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倍.版权申明本文部分内容，包括文字、图片、以及设计等在网上搜集整理.版权为个人所有This article includes some parts, including text, pictures, and design. Copyright is personal ownership.用户可将本文地内容或服务用于个人学习、研究或欣赏，以及其他非商业性或非盈利性用途，但同时应遵守著作权法及其他相关法律地规定，不得侵犯本网站及相关权利人地合法权利.除此以外，将本文任何内容或服务用于其他用途时，须征得本人及相关权利人地书面许可，并支付报酬.Users may use the contents or services of this article for personal study, research or appreciation, and other non-commercial or non-profit purposes, but at the same time, they shall abide by the provisions of copyright law and other relevant laws, and shall not infringe upon the legitimate rights of this website and its relevant obligees. In addition, when any content or service of this article is used for other purposes, written permission and remuneration shall be obtained from the person concerned and the relevant obligee.转载或引用本文内容必须是以新闻性或资料性公共免费信息为使用目地地合理、善意引用，不得对本文内容原意进行曲解、修改，并自负版权等法律责任.Reproduction or quotation of the content of this article must be reasonable and good-faith citation for the use of news or informative public free information. It shall not misinterpret or modify the original intention of the content of this article, and shall bear legal liability such as copyright.。