大学精品课【工程计算基础】工程计算3线性代数方程组的数值解法
工程经济学课件 5.1线性方程组的数值解法
(2)
22
23
(2)
(2)
2n
2
0
0
( 3 )
(3)
(3)
a a b 33
3n
3
0
0
( 3 )
(3)
(3)
a a b n3
nn
n
令mi 2
则
a(2) i2
a(2) 22
,i
2,3,..., n
a a m a ( 3 )
(2)
ij
ij
(2) i2 2j
(i 3,...,n; j 3,...,n)
b b m b ( 3 )
(2)
i
i
(2) i2 2
(i 3,...,n)
• 重复上述过程, 最后得
a1(11 )
a( 1 ) 12
a( 2 ) 22
... ...
... ...
...
a( 1 ) 1n
...
a( 2 ) 2n
b( 1
b( 2
1 2
) )
[ A( n )b( n )]
n
)
令mi1
则
a(1) i1
a(1) 11
,i
2,3,..., n
a a m a (2)
(1)
ij
ij
(1) i1 1 j
(i 2,...,n; j 2,...,n)
b b b m (2)
(1)
(1)
i
i
1 i1
(i 2,...,n)
第二步：当a2(22)
0,第2行
a(2) i2 a(2) 22
a(2) 2n
a(2) nn
线性代数方程组的数值解法讨论
线性代数方程组的数值解法讨论解线性方程组的方法，主要分为直接方法和迭代方法两种。
直接法是在没有舍入误差的假设下能在预定的运算次数内求得精确解。
而实际上，原始数据的误差和运算的舍入误差是不可以避免的，实际上获得的也是近似解。
迭代法是构造一定的递推格式，产生逼近精确解的序列。
对于高阶方程组，如一些偏微分方程数值求解中出现的方程组，采用直接法计算代价比较高，迭代法则简单又实用，因此比较受工程人员青睐。
小组成员本着工程应用，讨论将学习的理论知识转变为matlab 代码。
讨论的成果也以各种代码的形式在下面展现。
1 Jacobi 迭代法使用Jacobi 迭代法，首先必须给定初始值，其计算过程可以用以下步骤描述： 步骤1 输入系数矩阵A ，常熟向量b ，初值(0)x ，误差限ε，正整数N ，令1k =.步骤2 (0)11ni i ij jj ii j i x b a x a =≠⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦∑，(0)j x 代表(0)x 的第j 个分量。
步骤3 计算11ni i ij j j ii j i y b a x a =≠⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦∑，判断1max i i i n x y ε≤≤-<，如果是，则结束迭代，转入步骤5；否则，转入步骤4。
步骤4 判断k N =？如果是，则输出失败标志；否则，置1k k =+，i i x y ⇐，1,2,,i n =，转入步骤2。
步骤5 输出12,,n y y y 。
雅可比迭代代码function [x,k]=Fjacobi(A,b,x0,tol)% jacobi 迭代法 计算线性方程组% tol 为输入误差容限，x0为迭代初值max1= 300; %默认最多迭代300，超过要300次给出警告 D=diag(diag(A)); L=-tril(A,-1);U=-triu(A,1); B=D\(L+U); f=D\b; x=B*x0+f;k=1; %迭代次数while norm(x-x0)>=tol x0=x;x=B*x0+f; k=k+1;if(k>=max1)disp('迭代超过300次，方程组可能不收敛'); return; end%[k x'] %显示每一步迭代的结果 End2 高斯赛德尔迭代由Jacobi 迭代法中，每一次的迭代只用到前一次的迭代值，若每一次迭代充分利用当前最新的迭代值，即在计算第i 个分量(1)k i x +时，用最新分量11()k x +，12()k x +…(1)1k i x +-代替旧分量)1(k x ', )2(k x …)3(k x 就得到高斯赛德尔迭代格式，其数学表达式为：1(1)(1)()111(1,2,,)i n k k k ii ij j ij j j j i ii xb a x a x i n a -++==+⎛⎫=--= ⎪⎝⎭∑∑具体形式如下：()()()(1)()()()11221331111(1)(1)()()22112332222(1)(1)(1)(1)(1)112233,11111k k k k n n k k k k n n k k k k k n n n n n n n n nnx a x a x a x b a x a x a x a x b a x a x a x a x a x b a ++++++++--=----+=----+⋯⋯⋯⋯⋯⋯=-----+矩阵形式表示为：()(1)1(1)()(0,1,2,,),k k k k n +-+=++=x D Lx Ux b将(1)(1)()(0,1,2,,)k k k k n ++=++=Dx Lx Ux b 移项整理得： (1)1()1()()(0,1,2,,))k k x D L Ux D L b k n +--=-+-=记11(),()--=-=-M D L U g D L b ，则(1)()k k x x +=+M g高斯塞德尔迭代function [x,k]=Fgseid(A,b,x0,tol)%高斯-塞德尔迭代法 计算线性方程组 % tol 为误差容限max1= 300; %默认最高迭代300次D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); G=(D-L)\U; f=(D-L)\b; x=G*x0+f;k=1; while norm(x-x0)>=tol x0=x;x=G*x0+f; k=k+1;if(k>=max1)disp('迭代次数太多，可能不收敛'); return; end% [k,x'] %显示每一步迭代结果 End3 超松弛迭代法在工程中最常遇到的问题便是线性代数方程组的求解,而线性代数方程组的求解一般可以分为两类,一类是直接法(精确法),包括克莱姆法则方法、LD 分解法等，另一类是迭代法（近似法），包括雅克比迭代法、高斯迭代法、超松弛迭代法等。
实验5_线性代数方程组的数值解法
实验5 线性代数方程组的数值解法化工系 毕啸天 2010011811【实验目的】1. 学会用MATLAB 软件数值求解线性代数方程组，对迭代法的收敛性和解的稳定性作初步分析；2. 通过实例学习用线性代数方程组解决简化的实际问题。
【实验内容】题目3已知方程组Ax=b ，其中A ，定义为⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--------------=32/14/12/132/14/14/12/132/14/14/12/132/14/12/13A试通过迭代法求解此方程组，认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性质对收敛速度的影响。
实验要求：（1）选取不同的初始向量x (0)和不同的方程组右端项向量b ，给定迭代误差要求，用雅可比迭代法和高斯-赛德尔迭代法计算，观测得到的迭代向量序列是否均收敛？若收敛，记录迭代次数，分析计算结果并得出你的结论；（2）取定右端向量b 和初始向量x (0)，将A 的主对角线元素成倍增长若干次，非主对角线元素不变，每次用雅可比迭代法计算，要求迭代误差满足 ，比较收敛速度，分析现象并得出你的结论。
3.1 模型分析选取初始向量x(0) =(1,1,…,1)T ，b=(1,1,…,1)T ，迭代要求为误差满足 ，编写雅各比、高斯-赛德尔迭代法的函数，迭代求解。
3.2 程序代码function x = Jacobi( x0,A,b,m ) D=diag(diag(A)); U=-triu(A,1); L=-tril(A,-1); B1=D\(L+U); f1=D\b; x(:,1)=x0;x(:,2)=B1*x(:,1)+f1; k=1;while norm((x(:,k+1)-x(:,k)),inf)>m x(:,k+2)=B1*x(:,k+1)+f1; k=k+1; endendfunction x = Gauss( x0,A,b,m )D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B2=(D-L)\U;f2=(D-L)\b;x(:,1)=x0;x(:,2)=B2*x(:,1)+f2;k=1;while norm((x(:,k+1)-x(:,k)),inf)>mx(:,k+2)=B2*x(:,k+1)+f2;k=k+1;endendA1=3.*eye(20,20);A2=sparse(1:19,2:20,-1/2,20,20);A3=sparse(1:18,3:20,-1/4,20,20);AA=A1+A2+A3+A2'+A3';A=full(AA);b=ones(20,1); %输入自选右端项向量b x0=ones(20,1); %输入自选初始向量x0 m=1e-5;x1=Jacob(x0,A,b,m);x2=Gauss(x0,A,b,m);结果输出数据：3.3.1将x0各分量初值置为0【分析】从数据中可以看出，当迭代的初值变化了，达到相同精度所需要的迭代次数也变化了。
第3章 线性代数方程组的数值解法
第k步消元：
丛rk+1,rk+2, ,rn中消去xk项，条件akk(k)≠0，使得
A(k+1)x = b(k+1)A(k)x = b(k)
其中
第3章 线性代数方程组的数值解法Gauss消去法
(1 a11) (1 a12)
(k a kk )
( a11) n ( a 22 ) n
( a nn,1n)1 1
( a nn 1) 1n (n a nn)
b1(1) (1) b1 ( bnn11) (n) bn
(3.2.5)
由(3.2.4)式按倒序可方便的求出解向量x：
xn
( n) bn
( n) ann
( ( ( xn1 bnn11) ann,1n) xn ann,1n)1 1 1
ri(k)likrk(k) ri(k+1)，i= k+1,k+2, ,n
(3.2.2)
以矩阵[A(k)，b(k)]中的第k行乘以-lik加到i行，即 其中第i行 aij(k+1) = aij(k) likakj(k)，i, j= k+1,k+2,,n bi(k+1) = bi (k) likbk (k)，i= k+1,k+2,,n 当完成第k=n1步时， A(1)变为上三角阵A(n) ，Gauss消元过程
b1(1) ( 2) b2 ( bn2 )
第3章 线性代数方程组的数值解法Gauss消去法
具体方法：
-(r1(1)/a11(1))a21(1)加到第2行， -(r1(1)/a11(1))a31(1)加到第3行，，
线性代数方程组的数值解法_百度文库
线性代数方程组的数值解法【实验目的】1. 学会用MATLAB 软件数值求解线性代数方程组，对迭代法的收敛性和解的稳定性作初步分析；2. 通过实例学习用线性代数方程组解决简化的实际问题。
【实验内容】【题目1】通过求解线性方程组A1x=b1和A2x=b2，理解条件数的意义和方程组的性态对解的影响。
其中A1是n阶范德蒙矩阵，即⎡1x0⎢1x1⎢A1=⎢⎢⎢⎣1xn-12x0x12 2xn-1n-1⎤ x0⎥ x1n-1⎥1,...,n-1 ，xk=1+0.1k，k=0,⎥ n-1⎥ xn-1⎥⎦A2是n阶希尔伯特矩阵，b1，b2分别是A1，A2的行和。
（1）编程构造A1（A2可直接用命令产生）和b1，b2；你能预先知道方程组A1x=和A2x=。
b2的解吗？令n=5，用左除命令求解（用预先知道的解可检验程序）b1（2）令n=5,7,9,…，计算A1，A2的条件数。
为观察它们是否病态，做以下试验：b1，b2不变，A1和A2的元素A1(n,n)，A2(n,n)分别加扰动ε后求解；A1和A2不变，b1，b2的分量b1(n)，分析A和b的微小扰动对解的影响。
b2(n)分别加扰动ε求解。
ε取10-1010，-8，10-6。
（3）经扰动得到的解记做x~，计算误差-x~x，与用条件数估计的误差相比较。
1.1构造A1，A2和b1，b2首先令n=5，构造出A1，A2和b1，b2。
首先运行以下程序，输出A1。
运行以下程序对A1,A2求行和：由于b1，b2分别是A1，A2的行和，所以可以预知x1=运行下列程序，用左除命令对b1，b2进行求解：得到以下结果： T。
x2=(1,1, ,1)1.2 计算条件数并观察是否为病态1.不加扰动，计算条件数。
运行以下程序：由此可知，A1，A2的条件数分别是3.574∗10, 4,766∗10。
2.b1，b2不变，A1(n,n)，A2(n,n)分别加扰动（1）n=5时设x11,x12,x13分别为A1添加扰动10−10，10−8,10−6后的解。
【推荐】数值计算方法：第二章 线性方程组的数值解法.ppt
2.4 经典迭代法(Classic Iterative Methods)
迭代法思想：
42
2.4.1 Jacobi迭代法 (以对角元为分母)
43
建立迭代格式
将初值
代入后迭代得
将上述过程一般化
44
以分量表示方程组得
对角元对应的量 移到左边，其它 量在右边便得 ：
从而可建立迭代格式
Jacobi迭代
20
比较第k列, j=k,k+1,…,n
21
比较第k行, i=k+1,…,n
22
先解Ly = b 再解Ux = y
23
例 2.4 试用Crout分解法解线性方程组
解
24
25
26
2.2.2 对称正定矩阵分解法
若A 为对称正定矩阵，则容易证明存在下三角
矩阵L，使得
。这称为矩阵的乔里斯
基（Cholesky）分解。
50
定理 2.1 设方程组Ax=b的精确解为x*。如果存在一 个矩阵范数使得（2.4.9）中的迭代矩阵满足条件
则由(2.4.9)确立的迭代任何初始向量均收敛。且
证 迭代式相减取范数得
进一步递推得
定理得证。
51
利用定理2.1很容易判别迭代法的收敛性。以常用 矩阵范数为例，有下列结论。 推论 2.1 若(2.4.9)迭代矩阵 满足条件
解 相应的迭代公式为
Jacobi迭代
令 由Jacobi迭代得
由Gauss-Seidel迭代得
Gauss-Seidel迭代 取四位小数迭代计算
49
2.4.3 一般 迭代法的收敛性
定义3.2 设
相应的迭代格式为
如果存在某个向量范数使得
第3章_线性代数方程组的数值解法_new
其中biT为A的第i行的n维行向量.
上页 下页
矩阵的基本运算：
(1) 矩阵加法
C A B cij aij bij ( A, B, C R
(2) 矩阵与标量的乘法
mn
).
C A cij aij ( A, C Rmn ,是一个数).
(3) 矩阵与矩阵的乘法
(1) b1 (1) b2
(1) (1) (1) a12 a1 x b n 1 1 (1) (1) (1) a 22 a 2 n x 2 b2 . (1) ( 1 ) (1) am a mn 2 x n bm
(1) bm
上页
下页
将(2.1)记为A(1)x=b(1)，其中
(1) a11 (1) a 21 a (1) m1 (1) (1) a11 a12 a1 n (1) (1) a 22 a2 a21 n (1) (1) a am a 2 mn m1
解为
b1 b2 l 21 y1 y1 , y2 , l11 l 22
, yn
bn l1 j y j
j 1
n 1
l nn
计算量（乘除法的主要部分）都为 n2/2. 因此，我们将一般的线性方程组化成等价的三 角形方程组来求解.
上页 下页
首先举一个简单的例子来说明消去法的基本思想. 例1 用消去法解方程组
上页 下页
x1 x 2 (m维列向量). m x R x x m A a1 a2 an ,
其中aj为A的第j列的m维列向量. 同理
数值计算方法-第3章--线性方程组的解法PPT课件
个顺序主子式
a a (1)
(1)
11
12
Dk
a(1) 21
a(1) 22
a(1) 1k
a(1) 2k
0
(k 1, 2,..., n 1).
a a (1)
(1)
k1
k2
a(1) kk
.
13
顺序Gauss消去法计算过程中的 akk(k) 称为主元素，在 第k步消元时要用它作除数,则可能会出现以下几种情况
.
是原方程组 Ax=b 的解向量。
27
对于
Ly =b
1
由
l21
1
l31
l32 1
y1 b1
y2
b2
y3
b3
ln1 ln2 lnn1 1 yn bn
.
解得
y1 yk
b1 bk
k 1 i 1
lki
yi
,
k 2,3,, n
28
对于 Ux =y
u11 u12 u1n x1 y1
2x3 6
⑤
x1 6 (x2 x3 ) 1
x2 x3 5 / 4 2
x3 (6) / (2) 3
用x3, x2的值求x1 把x3的值代入②求x2
.
8
从下向上逐步求解
对应的增广矩阵的变化
1 1 1 6 1 1 1 6
( A | b) 0
4
1 5 0
4
1
5
2 2 1 1 0 4 1 11
0.8334
5.910
12.10
0.0120 0.0100 0.1670 0.6781
3200
1200
4.200 981.0
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
b11x1 b12 x2
b22 x2
b1,n x 1 n1 b1n xn g1 b2,n1xn1 b2n xn g2
b x n1,n1 n1 bn1,n xn gn1
bnn xn gn
这个过程称为消元过程。 由此可逐个求出xn，xn-1 ，… ，x1 ，这个过程称为回代过程
a(0) 11
l21
a(0) 12
a(1) 22
a(0) 13
a (1) 23
l31
l32
a(2) 33
ln1 ln2 ln3
A的行列式的值为
|
A
|
a a a (0) (1) (2) 11 22 33
2020/9/29
a(0) 1n
a (1) 2n
a(2) 3n
a ( n 1) nn
n
a(n1) nn
在每一次消元前，都进行上述的列选主元的过程 。这样的
解可由线性方程组
得到。
A(n)x=b(n)
2020/9/29
27
3.2 高斯消去法
以上消元过程中，对角线上的
a(k) kk
称为主元素。
这就是列选主元素方法的消元过程。
列选主元素方法中要多次交换系数矩阵的行，若共进行了m 次 ，则行列式
n
| A | (1)m
bk( k
)
b( k ) k 1
bn( k
)
11
3.2 高斯消去法
a(k) kk
0
计算因子
消元得到
a (1) 11
a (1) 12
a(2) 22
(
A(k 1) , b(k 1)
)
2020/9/29
lik
a(k) ik
/
a(k kk
)
(i k 1,
,n)
a (1) 1k
a(2) 1k
3.2 高斯消去法
其中
baii((jkk11))
a(k) ij
b(k ) i
lik
a(k kj
)
likbi(k )
(i, j k 1, k 2,
, n)
按上述做法，完成n-1次消元后，方程组化 成同解的上三角方程
A(n)x = b(n)
a (1) 11
( A(n) , b(n) )
a (1) 12
取五位有效数字，进行高斯消元计算，把第二个方程中的x1 消去，得到
0.0003x1
3.0000x2 9999.0x2
2.0001 6666.0
从第二个方程解得x2=0.6667，再代入第一个方程得x1=0。可 见，直接利用高斯消元法得到的结果与准确解相差较大。
2020/9/29
25
3.2 高斯消去法
a11x1 a12 x2 a21x1 a22 x2
a1n xn b1 a2n xn b2
an1x1 an2 x2 ann xn bn
用矩阵和向量的记号来表示，可写成 Ax b
a11 a12 A a21 a22
an1 an2
a1n
a2n
ann
x1
x
x2
高斯消去法的特点：消元和回代不同步！
2020/9/29
14
3.2 高斯消去法
消元过程可写为
[ A(n) , b(n) ] Ln1[ A(n1) , b(n1) ]
Ln L 1 n2 L2 L1[ A(0) , b(0) ]
Ln L 1 n2 L2L1[ A, b]
分开，有
A(n) Ln1Ln2 L2 L1 A
2020/9/29
21
3.2 高斯消去法
计算量：n阶行列式有n！项，每项计算n-1次乘法， 总乘法次数
Sn=(n+1)×n! ×(n-1)= (n+1)! ×(n-1) 当n=20时，乘法次数为9.7×1020次。若用每秒运 算10亿次乘法的计算机计算，约需3.08万年，这是无 法实现的。
若用Gauss消去法求解，其乘除法运算次数只需 3060次。
可分析两种解法不同结果的原因。
2020/9/29
26
3.2 高斯消去法
选主元 列选主元素方法
第一步时
在A(1)的第一列元素
ai(11)(i=1,2,
…,n)中，选主元素
a(1) i11
，使
a(1) i11
max
1in
|
a(1) i1
|
并交换(A(1)，b(1))的第一行与第i1行。元素记号保持不变。
a (1) 1,k 1
a(k) 1,k 1
a(k) kk
a(k) k ,k 1
a(k 1) k 1,k 1
a(k 1) n,k 1
a (1) 1n
a(2) 2n
a(k) kn
a(k 1) k 1,n
a(k 1) n,n
b(1) 1
b(2) 2
b(k ) k
bk(
k 1) 1
bn( k
1)
12
注意：右端的A(n)是一个上三角矩阵。
由于
Lk I lk ekT
因此
Lk1 I lk ekT
仍为单位下三角矩阵
2020/9/29
15
3.2 高斯消去法
令
L (Ln1Ln2
L2 L1)1 L11L21
L L 1 1 n2 n1
为如下形式的单位下三角矩阵
1
l21
L
l31
ln1,1 ln1
A LA(n) LU
即，一个方阵可以分解为一个单位下三角矩阵与一个
上三角矩阵的乘积。称为矩阵的LU分解，也称为 Doolittle分解。
方程
Ax=b
可以写成
LUx=b
等价于两个方程
Ly=b
Ux=y
2020/9/29
17
3.2 高斯消去法
存储上的特点。由于L是单位下三角矩阵，其对角元 素为 1，是已知量，可以不存储。因此可将L存于U空出的下三 角部分。即利用矩阵A的存储空间 ，有
xn
b1
b
b2
bn
设|A|≠0，则方程组有唯一解
2020/9/29
3
3.1 引言
直接法 将Ax=b化为等价的上三角形方程
Ux=L-1b 其中，U是上三角矩阵，L是下三角形矩阵。可用高斯消去 法得到 直接三角分解法：将A做三角分解
A=LU
然后求解下列两个三角形方程 Ly=b Ux=y
2020/9/29
4
3.1 引言
中小规模代数方程组，常选用直接法 如果右端项有多个，多采用直接法 大型稀疏矩阵，有些用迭代法
2020/9/29
5
3.2 高斯消去法
3.2.1 顺序消去过程和矩阵的LU三角分解 3.2.2 可行性和计算量 3.2.3 数值稳定性：选主元
2020/9/29
6
3.2 高斯消去法
3.2.1 顺序消去过程和矩阵的LU三角分解 高斯消去法把方程组转化为一个等价的三角形方程组
Summer Grass Fade
Arial Font Family
3 线性代数方程组的数值解法
• 3.1 引言 • 3.2 高斯消去法 • 3.3 矩阵的直接三角分解 • 3.4 方程组的形态、条件数 • 3.5 大型方程组的迭代法
2020/9/29
2
3.1 引言
本章研究的对象是n阶线性代数方程组
定理3.2.3 若A严格对角占优，则
a(k) kk
0
(k 1, 2,
,n)
2020/9/29
19
3.2 高斯消去法
高斯消去法的计算量
把A(1)化为等价的A(n)所需要的乘法运算量是
n (k 2 k) n (n2 1)
k 1
3
除法运算量是
n1
n
k (n 1)
k 1
2
回代过程需要的乘除法总运算量是n(n1)/2
b(1) 1
b(2) 2
a(2) nn
b(2) n
9
3.2 高斯消去法
第i行的消元公式为
baii((j22))
a(1) ij
b(1) i
li1a1(1j) li1b1(1)
(i, j 2,3,
, n)
矩阵表示
初等矩阵 其中
L1=I-l1e1T l1=(0,l21, l21, …, ln1)T
a ( i 1) ii
i 1
2020/9/29
28
3.2 高斯消去法
完全选主元方法
如果将选主元的范围扩大到余下的全部元素，即在
0
3
3
9
1 4 5 14 0 3 3 9 0 0 2 2
可以继续计算了
2020/9/29
24
3.2 高斯消去法
2) 主元过小，导致计算过程不稳定
例3.2.2 设有线性方程组
10..00000003xx11
3.0000x2 1.0000x2
2.0001 1.0000
该方程组的准确解为x1=1/3，x2=2/3。
2020/9/29
7
3.2 高斯消去法
记
a11 a12
(
A,
b)
a21
a22
an1
an2
且把(A，b)记为(A(1)，b(1))，即
(
A,
b)
(
A(1)
,
b(1)
)
a (1) 11
a (1) 21
a (1) 12
a (1) 22
a (1) n1
a (1) n2
2020/9/29