线性代数方程组数值解法及MATLAB 实现综述
廖淑芳 20122090 数计学院 12计算机科学与技术1班（职教本科） 一、分析课题
随着科学技术的发展，提出了大量复杂的数值计算问题，在建立电子计算机成为数值计算的主要工具以后，它以数字计算机求解数学问题的理论和方法为研究对象。其数值计算中线性代数方程的求解问题就广泛应用于各种工程技术方面。因此在各种数据处理中，线性代数方程组的求解是最常见的问题之一。关于线性代数方程组的数值解法一般分为两大类：直接法和迭代法。
直接法就是经过有限步算术运算，可求的线性方程组精确解的方法（若计算过程没有舍入误差），但实际犹如舍入误差的存在和影响，这种方法也只能求得近似解，这类方法是解低阶稠密矩阵方程组级某些大型稀疏矩阵方程组的有效方法。直接法包括高斯消元法，矩阵三角分解法、追赶法、平方根法。
迭代法就是利用某种极限过程去逐步逼近线性方程组精确解的方法。迭代法具有需要计算机的存储单元少，程序设计简单，原始系数矩阵在计算过程始终不变等优点，但存在收敛性级收敛速度问题。迭代法是解大型稀疏矩阵方程组（尤其是微分方程离散后得到的大型方程组）的重要方法。迭代法包括Jacobi 法SOR 法、SSOR 法等多种方法。
二、研究课题-线性代数方程组数值解法 一、 直接法 1、 Gauss 消元法
通过一系列的加减消元运算，也就是代数中的加减消去法，以使A 对角线以下的元素化为零，将方程组化为上三角矩阵；然后，再逐一回代求解出x 向量。
1.1消元过程
1. 高斯消元法(加减消元)：首先将A 化为上三角阵，再回代求解。
11121121222212n n n n nn n a a a b a a a b a a a b ?? ? ? ? ???L L M M O M M L (1)(1)(1)(1)(1)11121311(2)(2)(2)(2)222322(3)(3)(3)3333()()000000n
n n n n nn n a a a a b a a a b a a b a b ?? ? ? ? ?
? ???L L L M M M O
M M L 步骤如下：
第一步：1
11
1,2,,i a i i n a -?
+=L 第行第行 11121121222212
n n
n n nn n a a a b a a a b a a a b ??
? ?
?
?
??L L M M O M M L
111211(2)(2)(2)2222
(2)(2)(2)2
00n n
n nn n a a a b a a b a a b ??
?
? ? ???
L L M M O M M L
第二步：(2)2
(2)222,3,,i a i i n a -?+=L 第行第行
111211(2)(2)(2)2222
(2)(2)(2)2
00n n
n nn n a a a b a a b a a b ??
?
? ? ???
L L M M O M M L
11121311(2)(2)(2)(2)222322(3)(3)(3)3333(3)(3)(3)300000n n n
n nn n a a a a b a a a b a a b a a b ??
?
? ?
?
? ??
?
L
L L
M M M O
M M L 类似的做下去，我们有：
第k 步：()
()k ,1,,k ik
k kk
a i i k n a -?+=+L 第行第行。
n －1步以后，我们可以得到变换后的矩阵为：
11121311(2)(2)(2)(2)222322(3)(3)(3)3333()()000000n n n n n nn n a a a a b a a a b a a b a b ?? ?
? ? ? ?
???
L L L M M M O M M L
注意到，计算过程中()
k kk a 处在被除的位置，因此整个计算过程要保证它不为0。 所以，Gauss 消元法的可行条件为：()
0k kk
a ≠。 就是要求A 的所有顺序主子式均不为0，即
1111
det 0,1,,i i ii a a i n a a ??
?
≠= ? ???
L
M O M L L 因此，有些有解的问题，不能用Gauss 消元求解。
另外，如果某个()
k kk a 很小的话，会引入大的误差。
例 用Gauss 消去法解方程组：
（1）1234123341518311151111631112x x x x -?????? ? ? ?---- ? ? ?= ? ? ? ? ? ?-?????? （1）对增广矩阵进行初等变换
122133144
3
()21
()121
()4
123
34153715123341505
22218311155321911116044343111277
70
0444E E E E E E E E E +→-+→-+→-??
??-? ?- ? ?---- ? ???????→ ? ? ? ?
?-?? ?-- ??
?
233244344
5
()6
77()()6
11
123341233
4151537370505151522222211
29
11
29
00
0111136367357910000
003633E E E E E E E E E +→+→-+→--????
?
? ? ?-- ? ? ? ??????→??????→ ? ?
?
? ?
?
? ????
?
得等价方程组
123423434412334153715
52221129
1136
910
33x x x x x x x x x x -++=??
?-++=??
?+=??
?=?? 回代得40x =，33x =，22x =，11x =。
第一步：将(1)/3使x 1的系数化为1,再将(2)、(3)式中x 1的系数都化为零，即由(2)-2×(1)(1)得
)
1(321)1( (231)
32=++x x x )1(32)2( (03)
4
32=+x x
由(3)-4×(1)(1)得
第二步：将(2)(1)除以2/3，使x
2
系数化为1，得
再将(3)(1)式中x
2
系数化为零，由(3)(1)-(-14/3)*(2)(2)，得
第三步：将(3)(2)除以18/3，使x
3
系数化为1，得
经消元后，得到如下三角代数方程组：
1.2回代过程
由(3)(3)得 x
3=1,将x
3
代入(2)(2)得x
2
=-2,将x
2
、x
3
代入(1)(1)得x
2
=1,
所以，本题解为[x]=[1,2,-1]T
1.3高斯消元的公式
综合以上讨论，不难看出，高斯消元法解方程组的公式为第一步，消元
（1）令
a
ij (1) = a
ij
, （i,j=1,2,3,…,n）
b
i (1) =b
i
, （i=1,2,3,…,n）
（2）对k=1到n-1，若a
kk
(k)≠0，进行
l ik = a
ik
(k) / a
kk
(k) , （i=k+1,k+2,…,n）
a ij (k+1) = a
ij
(k) - l
ik
* a
kj
(k), （i,j= k+1,k+2,…,n）
b i (k+1) = b
i
(k) - l
ik
* b
k
(k), （i= k+1,k+2,…,n）
第二步，回代
若a
nn
(n) ≠ 0
x n = b
n
(n) / a
nn
(n)
x i = (b
i
(i) – sgm(a
ij
(i) * x
j
)/- a
ii
(i) ,（i = n-1,n-2,…,1）,( j =
)1(
3
2
)3(
......
6
3
10
3
14
-
=
-
-x
x
)2(
3
2
)2(
......
2=
+x
x
)2(
3
)3(
......
6
3
18
-
=
x
)3(
3
)3(
......
1
-
=
x
i+1,i+2,…,n )
2 、LU分解法
求解线性代数方程组除了高斯消元法外，还常用LU分解法（三角形分解法）。LU分解法的优点是当方程组左端系数矩阵不变，仅仅是方程组右端列向量改变，即外加激励信号变化时，能够方便地求解方程组。矩阵的三角分解法可分为直接三角分解法,列主元三角分解法,平方根法,三对角方程组的追赶法。下面讨论直接三角分解法。
设n阶线性方程组Ax=b 。假设能将方程组左端系数矩阵A，分解成两个三角阵的乘积，即A=LU ，式中，L为主对角线以上的元素均为零的下三角矩阵, 且主对角线元素均为1的上三角矩阵；U为主对角线以下的元素均为零
所以有，LUx=b 令 Ux=y，则 Ly=b
由A=LU,由矩阵的乘法公式： a
1j = u
1j
, j=1,2,…,n
a i1 = l
i1
u
11
, i=1,2,…,n
推出 u
1j = a
1j
, j=1,2,…,n
l i1 = a
i1
/u
11
, i=1,2,…,n
这样就定出了U的第一行元素和L的第一列元素。
设已定出了U的前k-1行和L的前k-1列，现在确定U的第k行和L的第k列。由矩阵乘法：
当r>k时，l
kr =0, 且l
kk
=1，因为
∑
=
=
n
r
rj
kr
kj
u
l
a
1
∑
=
+
=
n
r
rj
kr
kj
kj
u
l
u
a
1
所以，
同理可推出计算L 的第k 列的公式： 因此得到如下算法——杜利特（Doolittle ）算法:
(1)将矩阵分解为A=LU,对k=1,2,…,n ；j=k,k+1,…n ； i=k,k+1,…n ；
公式1
(2)解L y =b
(3)解U x =y
对大规模稀疏问题，如果能够通过调整方程及未知量的顺序使得方程组的系数矩阵成带状结构，则对系数矩阵使用通常的LU 分解，可以保障单位下三角矩阵L 及上三角矩阵U 仍为带状结构.
3、直接三角分解法
Gauss 消去法还有许多变形，有些变形是为了利用特殊技巧减少误差，把Gauss 消去法改写为更紧凑的形式，还有一些变形时根据某类矩阵的特性作一些修正和简化，这些方法可统称为直接三角分解法。
矩阵的三角分解
设A 的顺序主子式0(1,2,,)k k n ?≠=L ，则可建立线性方程组Ax b =的Gauss 消去法与矩阵分解的关系，即矩阵A 的LU 分解。这个问题前面已经讲的比较详
∑=+=-=n
r rj
kr kj kj n
k k j u l a u 1,...,1,∑=+=-=n
r kk
rj kr ik ik
n
k k i u u l a l 1
,...,1,/)(∑-==-=1
1,...,2,1:2k r r
kr k k n
k y l b y 公式∑+=-=-
=n
k r kk
r
kr k k n n k u x u
y x 1
1
,...,1,/)(3：公式???
?
?
?
???
=-=-=∑∑==1/)(u 11
kk kk
n
r rj ky ik jk n
r yj
ki kj ki l u u l a l u l a
细了，此处不再赘述。
Doolittle 分解法
首先假设A 的顺序主子式k ?都不为零，则A 可作Doolittle 分解，即
A LU =，其中L 是单位下三角阵，有1ii l =，i j <时0ij l =；U 是上三角阵，i j >时0ij u =。仔细写出为
111211112121
22
22122212
1211
1n n n n n n nn n n nn a a a u u u a a a l u u A LU a a a l l u ?????? ?
??
? ? ???
==
= ? ???
?
?????????
L L L L L L L L L L O O L L
L （2.11） 在前面逐步推导L 和U 的元素公式都要借助于有关的()
k ij a 来表示。现在强调
指出，只要从给定的A 通过比较（2.11）式的两边就可能逐步地把L 和U 构造出来，而不必利用Gauss 消去法的中间结果，这种方法称为Gauss 消去法的紧凑格式。
根据矩阵的乘法规则，比较（2.11）式的两边可得
min(,)
1
i j ij ik kj k a l u ==
∑
，,1,2,,i j n =L （2.12）
先算出U 的第1行，在（2.12）令1i =，由111l =，10(1)j l j n =<≤可得
11j j u a =，1,2,,j n =L
接着在（2.12）中令1j =，由1111i i a l u =?，从而算出L 的第1列
1111111//i i i l a u a a ==，2,,i n =L
若U 的第1至1k -行和L 的第1至1k -列已经算出，则由（2.12）可计算出
U 的第k 行元素
1
1k kj kj kr rj r u a l u -==-∑，,1,,j k k n =+L （2.13）
接着再计算出L 的第k 列元素
1
1
k ik ir rk
r ik kk
a l u l u -=-=
∑，1,,j k n =+L （2.14）
解方程组Ax b =就化为求解LUx b =，分两步求解。第一步解Ly b =，其中L 为下三角阵，只要用逐次向前代入的方法；第二步解Ux y =，其中U 为上三角阵，用逐次向后代入方法即可解除x 。
例 用Doolittle 方法求解：
12346
2
1
16241011141510135x x x x -?????? ? ? ?
- ? ? ?
=
? ? ?- ? ? ?
---?????? 解：第1步算出L 和U ：
111311
16
5119
1610
37L ?? ? ? ? ?=
? ? ?-
- ???，6211102
133337910
1019174U -??
?
? ? ?=-
? ? ?
??
?
第2步解Ly b =得：
231916,3,,574T
y ?
?=- ??
?
第3步解Ux y =得：
(1,1,1,1)T x =--
二、 迭代法
迭代法具有的特点是速度快。与非线性方程的迭代方法一样，需要我们构造一个等价的方程，从而构造一个收敛序列，序列的极限值就是方程组的根。 对方程组：Ax b =做等价变换：x Gx g =+ 如：令A M N =-，则
11()Ax b M N x b Mx b Nx x M Nx M b --=?-=?=+?=+
则，我们可以构造序列(1)() k k x G x g +=+
若()*k x x →* **x G x g Ax b ?=+?= 则，我们可以构造序列(1)() k k x G x g +=+ 若()*k x x →* **x G x g Ax b ?=+?= 同时：(1)()()**(*)k k k x x Gx Gx G x x +-=-=-
1(0)(*)k G x x +==-L
所以，序列收敛0k G ?→ 与初值的选取无关 1. Jacobi 迭代
1111111
n n n nn n n a x a x b a x a x b
++=??
?
?++=?L M L 11221111
2
211233122211 111()
1() 1()
n n n n n n n n n n nn x a x a x b a x a x a x a x b a x a x a x b a ---?
=++-??
-?=+++-???
??
-?=++-??L L M L (1)
()()11221111
(1)()()()2
2112331222(1)()()11 111
()1
() 1
()k k k n n k k k k n n k k k n n n n n n nn
x a x a x b a x a x a x a x b a x a x a x b a +++---?=++-??
-?=
+++-?∴?
??
-?=
++-??
L L M L 格式很简单：
1
(1)
()()1
1
1()i n
k k k i
ij j ij
j
i j j i ii x a x a x
b a -+==+-=+-∑∑
迭代矩阵
记A D L U =--
1100nn a D a ?? ?= ? ???
O
211
1
0000n nn a L a a -??
?- ? ?= ? ? ?--??M O
O L
121100000n n n a a U a ---??
?
? ?= ?
- ?
???
L O M O 易知，Jacobi 迭代有
()D L U x b --= ()Dx L U x b =++ 11()x D L U x D b --=++
11100 B () , D L U I D A f D b ---∴=+=-=，0B 是简单迭代的迭代矩阵。
看上述公式和过程很抽象，我们来举个简单例子：
1111221331
21122223323113223333
a x a x a x
b a x a x a x b a x a x a x b ++=??
++=??++=? （0ii a ≠） 1
112213311
2221123322
33311322331()1()1
()x b a x a x a x b a x a x a x b a x a x a ?=--???=--???=--??
得
得上述变换的方程后，我们任取一向量(0)(0)(0)(0)123(,,)T
x x x x =作初始近似，代入， (1)(1)(1)(1)()()()()12
3123(,,),,(,,)T m m m m T
x x x x x x x x ==L 假设上述方程的准确解是123(,,)T a a a a =
那么如果lim ,m m x a →∞
=?我们就说上述迭代是收敛的。
2. Gauss －Seidel 迭代
在Gauss －Seidel 迭代中，使用最新计算出的分量值
(1)()()11221111
(1)(1)()()2
2112331222(1)(1)(1)11 111
()1
() 1
()k k k n n k k k k n n k k k n n n n n n nn
x a x a x b a x a x a x a x b a x a x a x b a ++++++---?=++-??
-?=
+++-?∴?
??
-?=
++-??L L M L 1
(1)
(1)()1
1
1()i n
k k k i
ij j ij
j
i j j i ii x a x a x
b a -++==+-=+-∑∑
易知，Gauss －Seidel 迭代有
()D L U x b --= ()Dx L U x b =++ 11()x D L U x D b --=++
1110 () , G D L U I D A g D b ---∴=+=-=，0G 是简单迭代的迭代矩阵。
二种方法都存在收敛性问题。 收敛条件
迭代格式收敛的充要条件是B
G
的谱半径<1。
有例子表明：Gauss-Seidel 法收敛时，Jacobi 法可能不收敛；而Jacobi 法收敛时， Gauss-Seidel 法也可能不收敛。
若为严格对角占优的话则是收敛的，如130.10.51
20.0250.3110.030.142
???? ? ?
? ? ? ? ?????
，
假如方程组不满足收敛，有时候对其进行变换，可以改变收敛性。 如求下述方程组的解：
(0)1807109,8,(0,0,0)9117T A b x -????
? ?
=-== ? ? ? ?--????
9117180,71098A b --???? ? ?=-= ? ? ? ?-????
格式
(1)()()
123(1)
(1)1
1(1)(1)1
11/9(7)1/8(7)1/9(8)
k k k k k k k x x x x x x x +++++?=++?=+??=+? 结果
(1)(2)(3)
(4)(0.7778,0.9722,0.9753)(0.9942,0.9993,0.9994)(0.9999,0.9999,0.9999)
(1.0000,1.0000,1.0000)T T T
T
x x x
x ====
0()1B ρ<，G ρ（）<1也满足收敛。
如122111221A --??
?=-- ? ?-??
Jacobi ，0B 的特征方程为221
102
2
λ
λ
λ
---=-
解得1230λλλ===，所以用Jacobi 迭代法收敛。
Gauss-Seidel ，0G 的特征方程为22
1022λλ
λλλλ
---=--
解得1230,22λλλ==-=-
0()21B ρ=->, 所以Gauss-Seidel 迭代法是不收敛的。
3. 松弛迭代 记()(1)()k k k x x x +?=- 则(1)()()k k k x x x +=+?
可以看作在前一步上加一个修正量。若在修正量前乘以一个因子ω，有
(1)()()k k k x x x ω+=+? (1)()(1)()()k k k k x x x x ω++=+-
对Gauss －Seidel 迭代格式(1)1(1)()()k k k x D Lx Ux b +-+=++
(1)()1(1)1()1()()k k k k k x x D Lx D Ux D b x ω+-+--=+++- (1)()1(1)1()1(1)()k k k k x x D Lx D Ux D b ωω+-+--=-+++
写成分量形式，有
(1)()()()111221111
(1)()()()()222112331222
(1)()()()11 111
(1)()(1)()
1
(1)()k k k k n n k k k k k n n k k k k n n n n n n n nn
x x a x a x b a x x a x a x a x b a x x a x a x b a ωωω
ωωω
+++---?=-+++-??
?=--+++-????
-?=-+++-??
L L M
L
关于直接法和迭代法的例题： 1用选主元三角分解法求解下列方程组
1234123
12341234332462421022431742268x x x x x x x x x x x x x x x +++=??++=??
+++=??+++=? 12343324624201022431742268x x x x ?????? ? ? ? ? ? ?=
? ? ?
? ? ?
?????? 3
324
624201022431742268?? ?
?
?
?
??
列主
4226
84
2268422681131363
136222420108111130131112243172233
3
324
633113
100
13422242??
????
?
?-- ? ?
? ?
? ?→→ ?- ?
? ? ? ??? ? ?
--????
这里100042268110003132
6,,1
1810001112333310001014
2L U y ??
???? ?
? ? ?- ? ? ?=== ?
?- ? ? ? ? ?-?? ?????
再通过Ux y =，得x 的解。
123424233x x x x ???? ? ?- ? ?= ? ? ? ?-??
?? 2用Jacobi 迭代法求解下列方程组，并且使精度为310-。
1234
0.240.0880.0930.1590.040.08420x x x -?????? ??? ?-= ???
? ??? ?-??????以4，3，4分别除3方程两边 12310.060.0220.0310.0530.010.0215x x x -?????? ??? ?-= ??? ? ??? ?-?????? 即11223300.060.0220.03
00.0530.010.0205x x x x x x -???????? ? ??? ?
=-+ ? ??? ? ? ??? ?-????????
有Jacobi 迭代式123()
(1)1(1)()2
(1)()300.060.0220.0300.0530.010.0205m m m m m m x x x x x x +++????-???? ? ? ? ?
=-+ ? ? ? ? ? ? ? ? ?-????????
可证明此迭代格式是收敛的，极限值即为解。
40.240.080.060.020.0930.150.030.050.040.0840.010.02λ
λλλλλ--???? ? ?-→- ? ? ? ?--????22
20.06
0.02
0.0018
0.050.0006
0(0.06)(30.010.0009)00
3(0.0018)λλλλ
λλλλλ?? ?- ? ?-+→- ? ?++- ?
?-?
? 解得1230.06,λλλ=-=
=
由此可知，用Jacobi 迭代法是收敛的。 取(0)(2,3,5)T x =为初始值，迭代得，
23(1)1(1)(1)00.060.0222 1.920.03
00.0533 3.190.010.02055 5.04x x x ??-????????
?
??? ? ?=-+= ? ??? ? ? ? ??? ? ? ?-??????????
23(2)1(2)(2)00.060.02 1.922 1.9090.03
00.05 3.193 3.1940.010.020 5.045 5.045x x x ??-???????? ?
??? ? ?=-+= ? ??? ? ? ? ??? ? ? ?-?
???????
??…
其实23(3)1(3)(3)00.060.02 1.9092 1.9090.03
00.05 3.1943 3.1940.010.020 5.0455 5.045x x x ??-????????
? ??? ? ?
=-+= ? ??? ? ? ? ??? ? ? ?-?????????? 即在310-内，(3)(2)x x =，可以停止计算了。所以(3)x 即为近似解。
分析：此题题目可以直接说是解方程组，然后求解过程我们可以先验算Jacobi 和Gauss －Seidel 迭代的收敛情况，再选择算法，一般来说Gauss －Seidel 迭代要比Jacobi 迭代快。
三、用Matlab 软件求解线性方程组 一。高斯消去法 1.顺序高斯消去法 直接编写命令文件 a=[] d=[]'
[n,n]=size(a); c=n+1 a(:,c)=d; for k=1:n-1
a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); 消去 end
x=[0 0 0 0]' 回带 x(n)=a(n,c)/a(n,n); for g=n-1:-1:1
x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g) end
2.列主高斯消去法
*由于“[r,m]=max(abs(a(k:n,k)))”返回的行是“k:n,k”内的第几行，所以要通过修正来把m 改成真正的行的值。该程序只是演示程序，真正机器计算不需要算主元素所在列以下各行应为零的值。
直接编写命令文件
a=[]
d=[] '
[n,n]=size(a);
c=n+1
a(:,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],:); 换行end
a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k,
k:c); %消去
end
end
x=[0 0 0 0]' 回带
x(n)=a(n,c)/a(n,n);
for g=n-1:-1:1
x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g)
end
二迭代法
J迭代公式
a=[5 2 1;-1 4 2;2 -3 10]
d=[-12;20;3]
x=[0;0;0]; 初始向量
stop=1.0e-4 迭代的精度
L=-tril(a,-1)
U=-triu(a,1)
D=inv(diag(diag(a)))
X=D*(L+U)*x+D*d; J迭代公式
n=1;
while norm(X-x,inf)>=stop
x=X;
X=D*(L+U)*x+D*d;
n=n+1;
end
X
n
G-S迭代公式
a=[5 2 1;-1 4 2;2 -3 10]
x=[0;0;0]; 初始向量
d=[-12;20;3]
stop=1.0e-4
L=-tril(a,-1)
U=-triu(a,1)
D=(diag(diag(a)))
X=inv(D-L)*U*x+inv(D-L)*d; G-S迭代公式
n=1;
while norm(X-x,inf)>= stop 时迭代中止否则继续x=X;
X=inv(D-L)*U*x+inv(D-L)*d;
n=n+1;
end
X
n
SOR迭代公式
a=[5 2 1;-1 4 2;2 -3 10]
x=[0;0;0]; 初始向量
d=[-12;20;3]
stop=1.0e-4
w=1.44; 松弛因子
L=-tril(a,-1)
U=-triu(a,1)
D=(diag(diag(a)))
X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*d SOR迭代公式
n=1;
while norm(X-x,inf)>=stop 时迭代中止否则继续
x=X;
X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*d;
n=n+1;end
X
n
3. a*x=d a=[5 2 1;-1 4 2;2 -3 10] d=[-12;20;3]
(1)考察用J、G-S及SOR解此方程组的收敛性.
(2)分别用雅可比迭代法、高斯-赛德尔迭代法及逐次超松弛迭代法解此方程组。要求时迭代中止，观察收敛速度。
(1) J迭代矩阵的谱半径：max(abs(eig(D*(L+U))))= 0.506
G-S迭代矩阵的谱半径：max(abs(eig(inv(D-L)*U)))= 0.2000 SOR迭代矩阵的谱半径：
max(abs(eig(inv(D-w*L)*((1-w)*D+w*U))))=0.9756
所以用J、G-S及SOR均收敛（且有）。
(2)
J： X =-4.0000 3.0000 2.0000 n =18
G-S： X =-4.0000 3.0000 2.0000 n =8
SOR()： X =-4.0000 3.0000 2.0000 n =482
分析：可见高斯-赛德尔迭代法要比雅可比迭代公式的收敛速度快。同时，如果超松弛迭代法的选取不当不但不会加速收敛还会对收敛速度造成很恶劣的结果。
四、结论
在科学研究和工程技术中有许多问题可归结为求解线性代数方程组的问题。本文主要讨论了直接法解线性方程组及迭代法解线性方程，讲解了各种直接法，如高斯消元法，三角分解法，追赶法等的基本思想和原理，介绍了病态方程和误差分析，了解了它们各自的优缺点及适用范围。可以通过计算方法进行一种比较完善的构造用计算机进行计算，使之更普遍化，能够有举一反三的思想，利用直接法解线性方程组解决一些实际中难解的问题。
一、上机目的 1、培养学生运用线性代数的知识解决实际问题的意识、兴趣和能力； 2、掌握常用计算方法和处理问题的方法; 二、上机内容 1、求向量组的最大无关组； 2、解线性方程组； 三、上机作业 1、设A=[2 1 2 4; 1 2 0 2; 4 5 2 0; 0 1 1 7]； 求矩阵A列向量组的一个最大无关组. >> A=[2 1 2 4;1 2 0 2;4 5 2 0;0 1 1 7] A = 2 1 2 4 1 2 0 2 4 5 2 0 0 1 1 7 >> rref(A) ans = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 所以矩阵A的列向量组的一个最大无关组就是它本身； 2、用Matlab解线性方程组 （1） >> A=[2 4 -6;1 5 3;1 3 2] A = 2 4 -6 1 5 3 1 3 2 >> b=[-4;10;5]
b = -4 10 5 >> x=inv(A)*b x = -3.0000 2.0000 1.0000 >> B=[3 41 -62;4 50 3;11 38 25] B = 3 41 -62 4 50 3 11 38 25 >> c=[-41;100;50] c = -41 100 50 >> x=inv(B)*c x = -8.8221 2.5890 1.9465 3、（选作）减肥配方的实现 设三种食物每100克中蛋白质、碳水化合物和脂肪的含量如下表，表中还给出了20世纪80年代美国流行的剑桥大学医学院的简捷营养处方。现在的问题是：如果用这三种食物作为每天的主要食物，那么它们的用量应各取多少才能全面准确地实现这个营养要求？ 四、上机心得体会
D =) () ()(0)()() (001 11112 132 3122211331221 1 312a a a a a a a a a a a a a a a a a a a a a a a a k k k k k k k k ------------ 按第一列展开，再将各列的公因子提出来 D = ) ()()() () () (121323122211331221131 2a a a a a a a a a a a a a a a a a a a a a a a a k k k k k k k k ------------ =(a 2－a 1)(a 3－a 1)…(a k －a 1) 22322 32 111 ---k k k k k a a a a a a 得到的k －1阶范德蒙德行列式，由归纳假设知其值为 ∏≤<≤-k i j j i a a 2)( 于是 D =(a 2－a 1)(a 3－a 1)…(a k －a 1) ∏≤<≤-k i j j i a a 2)(= ∏≤<≤-k i j j i a a 1)( 因此，对于任意正整数n ≥2，范德蒙德行列式的展开式都成立。 证毕 例1.14 计算n 阶三对角行列式： D n = 2 1 120000 021000 12 1000 12------ 解 由行列式的性质1.4，将D n 的第一列的每个元看成两个元之和，得
D n = 2100 12000002100 120 00011----- +2 1 1200000 21000 12 1 00011------ 第一个行列式按第一列展开；第二个行列式从第一行开始依次加到下一行，得 D n =D n －1+ 1 110000 01000 110 00011 ---=D n －1+1 反复利用上面的递推公式，得到 D n =D n －1+1=D n －2+2=…=D 1+n －1=2+n －1=n +1 例1.15 计算n 阶行列式 D n = n a b b b a b b b a 21 (a i ≠b , i =1,2,…,n ) 解 对于这个行列式，采用一种“加边”的技巧。 D n =n a b b b a b b b a b b b 000121 第一行乘以（－1）加到其他各行上去，得
Matlab线性代数实验指导书 理学院线性代数课程组 二零零七年十月
目录 一、基础知识 (1) 1.1、常见数学函数 (1) 1.2、系统在线帮助 (1) 1.3、常量与变量 (2) 1.4、数组（矩阵）的点运算 (3) 1.5、矩阵的运算 (3) 二、编程 (4) 2.1、无条件循环 (4) 2.2、条件循环 (5) 2.3、分支结构 (5) 2.4、建立M文件 (6) 2.5、建立函数文件 (6) 三、矩阵及其运算 (7) 3.1、矩阵的创建 (7) 3.2、符号矩阵的运算 (11) 四、秩与线性相关性 (14) 4.1、矩阵和向量组的秩以及向量组的线性相关性 (14) 4.2、向量组的最大无关组 (14) 五、线性方程的组的求解 (16) 5.1、求线性方程组的唯一解或特解（第一类问题） (16) 5.2、求线性齐次方程组的通解 (18) 5.3、求非齐次线性方程组的通解 (19) 六、特征值与二次型 (22) 6.1、方阵的特征值特征向量 (22) 6.2、正交矩阵及二次型 (23)
一、基础知识 1.1常见数学函数 函数数学计算功能函数数学计算功能 abs(x) 实数的绝对值或复数的幅值floor(x) 对x朝-∞方向取整acos(x) 反余弦arcsinx gcd(m,n) 求正整数m和n的最大公约数acosh(x) 反双曲余弦arccoshx imag(x) 求复数x的虚部angle(x) 在四象限内求复数x的相角lcm(m,n)求正整数m和n的最小公倍 自然对数(以e为底数) asin(x) 反正弦arcsinx log(x) 常用对数(以 10 为底数) asinh(x) 反双曲正弦arcsinhx log10(x) atan(x) 反正切arctanx real(x) 求复数 x 的实部atan2(x,y) 在四象限内求反正切rem(m,n) 求正整数m和n的m/n之余数atanh(x) 反双曲正切arctanhx round(x) 对x四舍五入到最接近的整数 符号函数：求出 x 的符号ceil(x) 对x朝+∞方向取整 sign(x) conj(x) 求复数x的共轭复数 sin(x) 正弦sinx 反双曲正弦sinhx cos(x) 余弦cosx sinh(x) cosh(x) 双曲余弦coshx sqrt(x) 求实数x的平方根exp(x) 指数函数e x tan(x) 正切tanx fix(x) 对 x 朝原点方向取整 tanh(x) 双曲正切tanhx 如：输入 x=[-4.85 -2.3 -0.2 1.3 4.56 6.75]，则： ceil(x)= -4 -2 0 2 5 7 fix(x) = -4 -2 0 1 4 6 floor(x) =-5 -3 -1 1 4 6 round(x) = -5 -2 0 1 5 7 1.2 系统的在线帮助 1.2.1 help 命令： 1.当不知系统有何帮助内容时，可直接输入 help以寻求帮助: >> help(回车) 2.当想了解某一主题的内容时，如输入： >> help syntax (了解Matlab的语法规定) 3.当想了解某一具体的函数或命令的帮助信息时，如输入： >> help sqrt (了解函数sqrt的相关信息) 1.2.2 lookfor 命令 现需要完成某一具体操作，不知有何命令或函数可以完成，如输入： >> lookfor line (查找与直线、线性问题有关的函数) 1.3 常量与变量
数理方程数值解法与其在matlab软件上的实现张体强1026222 廖荣发1026226 [摘要] 数学物理方程的数值解在实际生活中越来越使用，首先基于偏微分数值解的思想上，通过matlab软件的功能，研究其数学物理方程的数值解，并通过对精确解和数值解进行对比，追究其数值解的可行性，在此，给出相关例子和程序代码，利于以后的再次研究和直接使用。 [关键字] 偏微分方程数值解matlab 数学物理方程的可视化 一：研究意义 在我们解数学物理方程，理论上求数学物理方程的定解有着多种解法，但是有许多定解问题却不能严格求解，只能用数值方法求出满足实际需要的近似解。而且实际问题往往很复杂，这时即便要解出精确解就很困难，有时甚至不可能，另一方面，在建立数学模型时，我们已作了很多近似，所以求出的精确解也知识推导出的数学问题的精确解，并非真正实际问题的精确解。因此，我们有必要研究近似解法，只要使所求得的近似解与精确解之间的误差在规定的范围内，则仍能满足实际的需要，有限差分法和有限元法是两种最常用的
求解数学物理方程的数值解法,而MATLAB 在这一方面具有超强的数学功能，可以用来求其解。 二：数值解法思想和步骤 2.1:网格剖分 为了用差分方法求解上述问题，将求解区域 {}(,)|01,01x t x t Ω=≤≤≤≤作剖分。将空间区间[0,1]作m 等分，将时 间[0,1]区间作n 等分，并记 1/,1/,,0,,0j k h m n x jh j m t k k n ττ===≤≤=≤≤。分别称h 和τ 为空间和 时间步长。用两簇平行直线,0,,0j k x x j m t t k n =≤≤=≤≤将Ω分割成矩形网格。 2.2：差分格式的建立 0u u t x ??-=??………………………………（1） 设G 是,x t 平面任一有界域，据Green 公式(参考数学物理方程第五章)： ( )()G u u dxdt udt udx t x Γ??-=--??? ? 其中G Γ=?。于是可将（1）式写成积分守恒形式： ()0udt udx Γ --=? (2) 我们先从（2）式出发构造熟知的Lax 格式设网格如下图所示
用高斯消元法求解线性代数方程组 1234111 5 -413-2823113-2104151 3-21719x x x x ??????????????????=?????? ?????? ?????? 1111X *??????=?????? （X*是方程组的精确解） 1 高斯消去法 1.1 基本思想及计算过程 高斯（Gauss ）消去法是解线性方程组最常用的方法之一，它的基本思想是通过逐步消元，把方程组化为系数矩阵为三角形矩阵的同解方程组，然后用回代法解此三角形方程组得原方程组的解。 为便于叙述，先以一个三阶线性方程组为例来说明高斯消去法的基本思想。 ??? ??=++II =++I =++III) (323034)(5 253)(6432321 321321x x x x x x x x x 把方程（I ）乘（2 3 - ）后加到方程（II ）上去，把方程（I ）乘（2 4- ）后加到方程（III ）上 去，即可消去方程（II ）、（III ）中的x 1，得同解方程组 ?? ? ??=+-II -=-I =++III) (20 223)(445.0)(6 4323232321x x x x x x x 将方程（II ）乘（ 5 .03 ）后加于方程（III ），得同解方程组： ?? ? ??-=-II -=-I =++III) (42)(445.0)(6432332321x x x x x x 由回代公式（3.5）得x 3 = 2，x 2 = 8，x 1 = -13。 下面考察一般形式的线性方程组的解法，为叙述问题方便，将b i 写成a i , n +1，i = 1, 2,…,n 。
线性代数方程组数值解法及MATLAB实现综述廖淑芳20122090 数计学院12计算机科学与技术1班（职教本科）一、分析课题 随着科学技术的发展，提出了大量复杂的数值计算问题，在建立电子计算机成为数值计算的主要工具以后，它以数字计算机求解数学问题的理论和方法为研究对象。其数值计算中线性代数方程的求解问题就广泛应用于各种工程技术方面。因此在各种数据处理中，线性代数方程组的求解是最常见的问题之一。关于线性代数方程组的数值解法一般分为两大类：直接法和迭代法。 直接法就是经过有限步算术运算，可求的线性方程组精确解的方法（若计算过程没有舍入误差），但实际犹如舍入误差的存在和影响，这种方法也只能求得近似解，这类方法是解低阶稠密矩阵方程组级某些大型稀疏矩阵方程组的有效方法。直接法包括高斯消元法，矩阵三角分解法、追赶法、平方根法。 迭代法就是利用某种极限过程去逐步逼近线性方程组精确解的方法。迭代法具有需要计算机的存储单元少，程序设计简单，原始系数矩阵在计算过程始终不变等优点，但存在收敛性级收敛速度问题。迭代法是解大型稀疏矩阵方程组（尤其是微分方程离散后得到的大型方程组）的重要方法。迭代法包括Jacobi法SOR法、SSOR法等多种方法。 二、研究课题-线性代数方程组数值解法 一、直接法 1、Gauss消元法 通过一系列的加减消元运算，也就是代数中的加减消去法，以使A对角线以下的元素化为零，将方程组化为上三角矩阵；然后，再逐一回代求解出x向量。
1.1消元过程 1. 高斯消元法(加减消元)：首先将A 化为上三角阵，再回代求解。 11121121222212n n n n nn n a a a b a a a b a a a b ?? ? ? ? ???(1)(1)(1)(1)(1)11121311(2)(2)(2)(2)222322 (3)(3)(3)3333()()000 00 n n n n n nn n a a a a b a a a b a a b a b ?? ? ? ? ? ? ??? 步骤如下： 第一步：1 11 1,2,,i a i i n a -? +=第行第行 11121121222212 n n n n nn n a a a b a a a b a a a b ?? ? ? ? ???1112 11(2)(2)(2)22 22 (2)(2)(2)2 00n n n nn n a a a b a a b a a b ?? ? ? ? ??? 第二步：(2)2 (2)222,3, ,i a i i n a -?+=第行第行 111211(2)(2)(2)2222(2)(2)(2)2 00n n n nn n a a a b a a b a a b ?? ? ? ? ???11 12 1311(2)(2)(2)(2)222322 (3)(3)(3)33 33(3)(3)(3)3 0000 0n n n n nn n a a a a b a a a b a a b a a b ?? ? ? ? ? ? ??? 类似的做下去，我们有： 第k 步：() ()k ,1, ,k ik k kk a i i k n a -?+=+第行第行。 n －1步以后，我们可以得到变换后的矩阵为： 11121311(2)(2)(2)(2)222322 (3)(3)(3)3333()()00000 n n n n n nn n a a a a b a a a b a a b a b ?? ? ? ? ? ? ?? ?
合肥学院 2018—2019学年第2学期 线性代数及应用 (模块) 实验报告 实验名称：线性代数MATLAB实验 实验类别：综合性 设计性□验证性 专业班级： 17通信工程（2）班 实验时间： 9-12周 组别：第组人数 3人 指导教师：牛欣成绩： 完成时间： 2019年 5 月9日
一. 小组成员 姓名学号具体分工 汪蔚蔚(组长) 1705022025 A报告最后的整合，编写，案例四的计算与应用 以及案例一的计算与证明 陶乐 1 1705022009 C案例二，化学方程式配平问题 程赢妹1505022036 A案例三，应用题灰度值的计算问题 二. 实验目的 1、案例一利用MATLAB进行线性代数计算，求出矩阵B 2、案例二利用MATLAB计算出每一个网格数据的值，然后每一个网格数据的值乘以256以后进行归一化处理，根据每个网格中的灰度值，绘制出灰度图像。 3、案例三利用MATLAB完成对化学方程式进行配平的应用 4、案例四利用MATLAB求极大线性无关组，并表示出其余向量 三. 实验内容 1、案例一： 0,1,0 ,=1,0,0, 0,0,0 A B AB BA A B ?? ?? =?? ?? ?? 已知矩阵和矩阵满足乘法交换律，即且求矩阵。 2、案例二 配平下列化学方程式： 3、案例三： 3*32 0.81.21.70.20.3 0.6021.61.20.6. 1MATLAB 2256MATLAB 给定一个图像的个方向上的灰度叠加值：沿左上方到右 下方的灰度叠加值依次为，，，，；沿右上方到左下 方的灰度叠加值依次为，。，，， ）建立可以确定网络数据的线性方程组，并用求解 ）将网络数据乘以，再取整，用绘制该灰度图像
一维抛物线偏微分方程数值解法（3） 上一篇参看一维抛物线偏微分方程数值解法（2）（附图及matlab程序） 解一维抛物线型方程(理论书籍可以参看孙志忠：偏微分方程数值解法） Ut-Uxx=0, 0 r=h2/(h1*h1); for(i=1:n) %外循环，先固定每一时间层，每一时间层上解一线性方程组% a(1)=0;b(1)=1+r;c(1)=-r/2;d(1)=r/2* (u(i+1,1)+u(i,1))+h2*f(i,j)... +(1-r)*u(i,2)+r/2*u(i,3); for(k=2:m-2) a(k)=-r/2;b(k)=1+r;c(k)=- r/2;d(k)=h2*f(i,j)+r/2*u(i,k)+(1-r)... *u(i,k+1)+r/2*u(i,k+2); %输入部分系数矩阵，为0的矩阵元素不输入% end a(m-1)=-r/2;b(m-1)=1+r;d(m-1)=h2*f(i,j)+r/2* (u(i,m+1)+u(i+1,m+1)... )+r/2*u(i,m-1)+(1-r)*u(i,m); for(k=1:m-2) %开始解线性方程组消元过程 a(k+1)=-a(k+1)/b(k); b(k+1)=b(k+1)+a(k+1)*c(k); d(k+1)=d(k+1)+a(k+1)*d(k); end u(i+1,m)=d(m-1)/b(m-1); %回代过程% for(k=m-2:-1:1) u(i+1,k+1)=(d(k)-c(k)*u(i+1,k+2))/b(k); end end for(i=1:n+1) for(j=1:m+1) p(i,j)=exp(x(j)+t(i)); %p为精确解 e(i,j)=abs(u(i,j)-p(i,j));%e为误差 end end [u p e x t]=CN(0.1,0.005,10,200);surf(x,t,e); shading interp; >> xlabel('x');ylabel('t');zlabel('e'); >> title('误差曲面') 第二章线性代数方程组的直接解法 教学目标： 1.了解线性代数方程组的结构、基本理论以及相关解法的发展历程； 2.掌握高斯消去法的原理和计算步骤，理解顺序消去法能够实现的条件，并在此基础上理解矩阵的三角分解（即LU分解），能应用高斯消去法熟练计算简单的线性代数方程组； 3.在理解高斯消去法的缺点的基础上，掌握有换行步骤的高斯消去法，从而理解和掌握选主元素的高斯消去法，尤其是列主元素消去法的理论和计算步骤，并能灵活的应用于实际中。 教学重点： 1. 高斯消去法的原理和计算步骤； 2. 顺序消去法能够实现的条件； 3. 矩阵的三角分解（即LU分解）； 4. 列主元素消去法的理论和计算步骤。 教学难点： 1. 高斯消去法的原理和计算步骤； 2. 矩阵的三角分解（即LU分解）； 3. 列主元素消去法的理论和计算步骤。 教学方法： 教具： 引言 在自然科学和工程技术中，许多问题的解决常常归结为线性方程组的求解，有的问题的数学模型中虽不直接表现为线性方程组，但它的数值解法中将问题“离散化”或“线性化”为线性方程组。例如，电学中的网络问题、船体数学放样中建立三次样条函数问题、最小二乘法用于求解实验数据的曲线拟合问题、求解非线性方程组问题、用差分法或有限元法求解常微分方程边值问题及偏微分方程的定解问题，都要导致求解一个或若干个线性方程组的问题。 目前，计算机上解线性方程组的数值方法尽管很多，但归纳起来，大致可以分为两大类：一类是直接法（也称精确解法）；另一类是迭代法。例如线性代数中的Cramer法则就是一种直接法，但其对高阶方程组计算量太大，不是一种实用的算法。实用的直接法中具有代表性的算法是高斯（Gauss）消元法，其它算法都是它的变形和应用。 在数值计算历史上，直接法和迭代法交替生辉。一种解法的兴旺与计算机的硬件环境和问题规模是密切相关的。一般说来，对同等规模的线性方程组，直接法对计算机的要求高于迭代法。对于中、低阶（200 n ）以及高阶带形的线性方程组，由于直接法的准确性和可靠性高，一般都用直接法求解。对于一般高阶方程组，特别是系数矩阵为大型稀疏矩阵的线性方程组用迭代法有效。 线性代数MATLAB 实验指导书 MATLAB 是Matrix Laboratory 的缩写，是一个集数值计算、图形处理、符号运算、文字处理、数学建模、实时控制、动态仿真和信号处理等功能为一体的数学应用软件，而且该系统的基本数据结构是矩阵，又具有数量巨大的内部函数和多个工具箱，使得该系统迅速普及到各个领域，尤其在大学校园里，许多学生借助它来学习大学数学和计算方法等课程，并用它做数值计算和图形处理等工作。我们在这里介绍它的基本功能，并用它做与线性代数相关的数学实验。 在正确完成安装MATLAB 软件之后，直接双击系统桌面上的MATLAB 图标，启动MATLAB ，进入MATLAB 默认的用户主界面，界面有三个主要的窗口：命令窗口（Commend Window ）, 当前目录窗口（Current Directory ）,工作间管理窗口（Workspace ）。 命令窗口是和Matlab 编译器连接的主要窗口，“>>”为运算提示符，表示Matlab 处于准备状态，当在提示符后输入一段正确的运算式时，只需按Enter 键，命令窗口中就会直接显示运算结果。 实验1 矩阵的运算，行列式 实验名称：矩阵的运算，行列式 实验目的：学习在matlab 中矩阵的输入方法以及矩阵的相关运算，行列式。 实验原理：介绍相关的实验命令和原理 (1)一般矩阵的输入 (2)特殊矩阵的生成 (3)矩阵的代数运算 (4)矩阵的特征参数运算 (5)数字行列式和符号行列式的计算 实验命令 1 矩阵的输入 Matlab 是以矩阵为基本变量单元的，因此矩阵的输入非常方便。输入时，矩阵的元素用方括号括起来，行内元素用逗号分隔或空格分隔，各行之间用分号分隔或直接回车。 例1 输入矩阵 ???? ? ??--=654301211A ，可以在命令窗口中输入 >>A=[1 1 2;-1 0 3;4 -5 6] A = 1 1 2 -1 0 3 4 - 5 6 2 特殊矩阵的生成 某些特殊矩阵可以直接调用相应的函数得到，例如： zeros(m,n) 生成一个m 行n 列的零矩阵 一维抛物线偏微分方程数值解法（1） 解一维抛物线型方程(理论书籍可以参看孙志忠：偏微分方程数值解法） Ut-Uxx=0, 0 线性代数方程组数值解法及MATLAB 实现综述 廖淑芳 20122090 数计学院 12计算机科学与技术1班（职教本科） 一、分析课题 随着科学技术的发展，提出了大量复杂的数值计算问题，在建立电子计算机成为数值计算的主要工具以后，它以数字计算机求解数学问题的理论和方法为研究对象。其数值计算中线性代数方程的求解问题就广泛应用于各种工程技术方面。因此在各种数据处理中，线性代数方程组的求解是最常见的问题之一。关于线性代数方程组的数值解法一般分为两大类：直接法和迭代法。 直接法就是经过有限步算术运算，可求的线性方程组精确解的方法（若计算过程没有舍入误差），但实际犹如舍入误差的存在和影响，这种方法也只能求得近似解，这类方法是解低阶稠密矩阵方程组级某些大型稀疏矩阵方程组的有效方法。直接法包括高斯消元法，矩阵三角分解法、追赶法、平方根法。 迭代法就是利用某种极限过程去逐步逼近线性方程组精确解的方法。迭代法具有需要计算机的存储单元少，程序设计简单，原始系数矩阵在计算过程始终不变等优点，但存在收敛性级收敛速度问题。迭代法是解大型稀疏矩阵方程组（尤其是微分方程离散后得到的大型方程组）的重要方法。迭代法包括Jacobi 法SOR 法、SSOR 法等多种方法。 二、研究课题-线性代数方程组数值解法 一、 直接法 1、 Gauss 消元法 通过一系列的加减消元运算，也就是代数中的加减消去法，以使A 对角线以下的元素化为零，将方程组化为上三角矩阵；然后，再逐一回代求解出x 向量。 1.1消元过程 1. 高斯消元法(加减消元)：首先将A 化为上三角阵，再回代求解。 11121121222212n n n n nn n a a a b a a a b a a a b ?? ? ? ? ???L L M M O M M L (1)(1)(1)(1)(1)11121311(2)(2)(2)(2)222322(3)(3)(3)3333()()000000n n n n n nn n a a a a b a a a b a a b a b ?? ? ? ? ? ? ???L L L M M M O M M L 步骤如下： 本章主要介绍方程根的有关概念，求方程根的步骤，确定根的初始近似值的方法（作图法，逐步搜索法等），求根的方法（二分法，迭代法，牛顿法，割线法，米勒（M üller ）法和迭代法的加速等）及其MATLAB 程序，求解非线性方程组的方法及其MATLAB 程序. 2.1 方程（组）的根及其MATLAB 命令 2.1.2 求解方程（组）的solve 命令 求方程f (x )=q (x )的根可以用MATLAB 命令: >> x=solve('方程f(x)=q(x)',’待求符号变量x ’) 求方程组f i (x 1,…,x n )=q i (x 1,…,x n ) (i =1,2,…,n )的根可以用MATLAB 命令: >>E1=sym('方程f1(x1,…,xn)=q1(x1,…,xn)'); ……………………………………………………. En=sym('方程fn(x1,…,xn)=qn(x1,…,xn)'); [x1,x2,…,xn]=solve(E1,E2,…,En, x1,…,xn) 2.1.3 求解多项式方程（组）的roots 命令 如果)(x f 为多项式，则可分别用如下命令求方程0)(=x f 的根，或求导数)('x f （见表 2-1）. 2.1.4 求解方程（组）的fsolve 命令 如果非线性方程（组）是多项式形式，求这样方程（组）的数值解可以直接调用上面已经介绍过的roots 命令.如果非线性方程（组）是含有超越函数，则无法使用roots 命令，需要调用MATLAB 系统中提供的另一个程序fsolve 来求解.当然，程序fsolve 也可以用于多项式方程（组），但是它的计算量明显比roots 命令的大. fsolve 命令使用最小二乘法（least squares method ）解非线性方程（组） (F X =)0 的数值解，其中X 和F （X ）可以是向量或矩阵.此种方法需要尝试着输入解X 的初始值（向量或矩阵）X 0，即使程序中的迭代序列收敛，也不一定收敛到(F X =)0的根（见例2.1.8）. fsolve 的调用格式： X=fsolve(F,X0) 输入函数)(x F 的M 文件名和解X 的初始值（向量或矩阵）X 0，尝试着解方程（组） 线性代数方程组求解 一、实验要求 编程求解方程组: 方程组1: 方程组2: 方程组3： 要求: 用C/C++语言实现如下函数： 1.bool lu(double* a， int＊ pivot, int n）； 实现矩阵的LU分解。 pivot为输出参数,pivot[0,n）中存放主元的位置排列. 函数成功时返回false，否则返回true。 2.bool guass（double const* lu， int const* p, double* b, int n）； 求线代数方程组的解 设矩阵Lunxn 为某个矩阵anxn 的LU 分解，在内存中按行优先次序存放。p[0，n)为LU 分解的主元排列.b 为方程组Ax=b 的右端向量.此函数计算方程组Ax=b 的解，并将结果存放在数组b ［0，n ）中.函数成功时返回false ，否则返回true 。 3。 void qr(double* a ， double ＊ d, int n);矩阵的QR 分解 假设数组anxn 在内存中按行优先次序存放。此函数使用HouseHolder 变换将其就地进行QR 分解。 d 为输出参数，d [0,n) 中存放QR 分解的上三角对角线元素。 4。 bool hshld(double const*qr ， double const*d, double*b ， int n); 求线代数方程组的解 设矩阵qrnxn 为某个矩阵anxn 的QR 分解,在内存中按行优先次序存放。d [0,n ） 为QR 分解的上三角对角线元素。b 为方程组Ax=b 的右端向量。 函数计算方程组Ax=b 的解,并将结果存放在数组b[0，n)中。 函数成功时返回false ，否则返回true 。 二、问题分析 求解线性方程组Ax=b ，其实质就是把它的系数矩阵A 通过各种变换成一个下三角或上三角矩阵，从而简化方程组的求解。因此,在求解线性方程组的过程中，把系数矩阵A 变换成上三角或下三角矩阵显得尤为重要,然而矩阵A 的变换通常有两种分解方法：LU 分解法和QR 分解法。 1、LU 分解法： 将A 分解为一个下三角矩阵L 和一个上三角矩阵U,即:A=LU ， 其中 L=??????? ?????1001 00 12121 n n l l l ， U=? ? ??? ? ??????nn n n u u u u u u 000 00222112 11 2、QR 分解法： 将A 分解为一个正交矩阵Q 和一个上三角矩阵R,即：A=QR 三、实验原理 解Ax=b 的问题就等价于要求解两个三角形方程组： ⑴ Ly=b，求y; ⑵ Ux=y,求x 。 设A 为非奇异矩阵，且有分解式A=LU , L 为单位下三角阵，U 为上三角 双曲型方程基于MATLAB 的数值解法 （数学1201,陈晓云,41262022） 一：一阶双曲型微分方程的初边值问题 0,01,0 1.(,0)cos(),0 1. (0,)(1,)cos(),0 1. u u x t t x u x x x u t u t t t ππ??-=≤≤≤≤??=≤≤=-=≤≤ 精确解为 ()t x cos +π 二：数值解法思想和步骤 2.1:网格剖分 为了用差分方法求解上述问题，将求解区域{}(,)|01,01x t x t Ω=≤≤≤≤作剖分。将空间区间[0,1]作m 等分，将时间[0,1]区间作n 等分，并记 1/,1/,,0,,0j k h m n x jh j m t k k n ττ===≤≤=≤≤。分别称h 和τ为空间和时 间步长。用两簇平行直线,0,,0j k x x j m t t k n =≤≤=≤≤将Ω分割成矩形网格。 2.2：差分格式的建立 0u u t x ??-=?? 2.2.1：Lax-Friedrichs 方法 对时间、空间采用中心差分使得 2h 1 1111)(2 1u u x u u u u u t u k j k j k j k j k j k j -+-++-= +=-= ????τ τ 则由上式得到Lax-Friedrichs 格式 1 11111()202k k k k k j j j j j u u u u u h τ+-+-+-+-+= 截断误差为 ()[]k k k j h j j R u L u Lu =- 1 11111()22k k k k k k k j j j j j j j u u u u u u u h t x τ+-+-+-+-??=+-+?? 23222 3 (),(0,0)26k k j j u u h O h j m k n t x ττ??= -=+≤≤≤≤?? 所以Lax-Friedrichs 格式的截断误差的阶式2()O h τ+ 令/s h τ=：则可得差分格式为 1111 11(),(0,0)222 k k k k k j j j j j s s u u u u u j m k n +--++=-+++≤≤≤≤ 0cos()(0)j j u x j m π=≤≤ 0cos(),cos(),(0)k k k m k u t u t k n ππ==-≤≤ 其传播因子为： ()()()e e G h i h i s h i h i σσσστσ---=-+e e 221， 化简可得： ()()()()()h s G h is h G στσσστ σsin 11,sin cos ,2 2 2--=-= 所以当1s ≤时，()1,≤τσG ,格式稳定。 * 2.2.2：LaxWendroff 方法 用牛顿二次插值公式可以得到LaxWendroff 的差分格式，在此不详细分析，它的截断误差为() h 2 2 +O τ ，是二阶精度；当2s ≤时，()1,≤τσG ， 格式稳定。在这里主要用它与上面一阶精度的Lax-Friedrichs 方法进行简单对比。 2.3差分格式的求解 Matlab部分函数名的义源 rand(m,n) random 随机 inv(a) inverse 逆矩阵 root 平方根sqrt(a) squared abs(a) absolute value 绝对值 det(a) determinant 行列式 rank(a) rank 秩 trace(a) trace 迹 rref(a) reduced row echelon form 最简行阶梯形 space 零核空间null(a) null sym(a) symbol 符号 orth(a) orthogonal 正交 norm norm 模 poly(a) polynomial 多项式 roots(p) root 根 eig(a) eigen- 特征的eigensys(a) eigen- system 特征的 线性代数部分词汇英汉对照 adjoint matrix 伴随矩阵 algebraic cofactor 代数余子式 augmented matrix 增广矩阵 block matrix 分块矩阵 basic solution set 基础解系 characteristic equation 特征方程 characteristic polynomial 特征多项式 coefficient matrix 系数矩阵 cofactor 余子式 column vector 列向量 canonical form [二次型的]标准形 cramer’s rule 克莱姆法则 determinant of order n n阶行列式 diagonal matrix 对角矩阵 dimension 维数 echelon form 阶梯形 eigenvalue 特征值 eigenvector 特征向量 elementary matrix 初等矩阵 elementary row operation 行初等变换 full rank 满秩 general solution 通解 gram-schmidt process 施密特正交化过程 identity matrix 单位矩阵 index of inertia 惯性指数 线性代数方程组求解 The Standardization Office was revised on the afternoon of December 13, 2020 线性代数方程组求解 一、实验要求 编程求解方程组： 方程组1： 方程组2： 方程组3： 要求： 用C/C++语言实现如下函数： 1.bool lu(double* a, int* pivot, int n); 实现矩阵的LU分解。 pivot为输出参数，pivot[0,n) 中存放主元的位置排列。 函数成功时返回false ，否则返回true 。 2. bool guass(double const* lu, int const* p, double* b, int n); 求线代数方程组的解 设矩阵Lunxn 为某个矩阵anxn 的LU 分解，在内存中按行优先次序存放。p[0,n)为LU 分解的主元排列。b 为方程组Ax=b 的右端向量。此函数计算方程组Ax=b 的解，并将结果存放在数组b[0,n)中。函数成功时返回false ，否则返回true 。 3. void qr(double* a, double* d, int n);矩阵的QR 分解 假设数组anxn 在内存中按行优先次序存放。此函数使用HouseHolder 变换将其就地进行QR 分解。 d 为输出参数，d [0,n) 中存放QR 分解的上三角对角线元素。 4. bool hshld(double const*qr, double const*d, double*b, int n); 求线代数方程组的解 设矩阵qrnxn 为某个矩阵anxn 的QR 分解，在内存中按行优先次序存放。d [0,n) 为QR 分解的上三角对角线元素。b 为方程组Ax=b 的右端向量。 函数计算方程组Ax=b 的解，并将结果存放在数组b[0,n)中。 函数成功时返回false ，否则返回true 。 二、问题分析 求解线性方程组Ax=b ，其实质就是把它的系数矩阵A 通过各种变换成一个下三角或上三角矩阵，从而简化方程组的求解。因此，在求解线性方程组的过程中，把系数矩阵A 变换成上三角或下三角矩阵显得尤为重要，然而矩阵A 的变换通常有两种分解方法：LU 分解法和QR 分解法。 1、LU 分解法： 将A 分解为一个下三角矩阵L 和一个上三角矩阵U ，即：A=LU, 其中 L=??????? ?????10010012121 n n l l l ， U=? ? ??? ? ??????nn n n u u u u u u 000 00222112 11 2、QR 分解法： 将A 分解为一个正交矩阵Q 和一个上三角矩阵R,即：A=QR 实验 5：线性代数方程组的数值解法 习题3： 已知方程组Ax b =,其中20*20 A R ∈，定义为： 3 1/21/41/231/21/41/41/231/21/41/41/231/21/4 1/23--????---? ???---??-???? ---?? --? ? 试通过迭代法求解此方程组，认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性质对 收敛速度的影响。实验要求： （1） 选取不同的初始向量x0和不同的方程组右端向量b ，给定迭代误差要求，用雅可比 迭代法和高斯-赛德尔迭代法计算，观测得到的迭代向量序列是否均收敛？若收敛，记录迭代次数，分析计算结果并得出结论； （2） 取定右端向量b 和初始向量x0，将A 的主对角线元素成倍的增长若干次，非主对角 元素不变，每次用雅可比迭代法计算，要求迭代误差满足(1)()5 ||||10k k x x +-∞-<,比 较收敛速度，分析现象并得出结论。 1、 程序设计（可直接粘贴运行） 1） Jacobi 迭代法 function y=jacobi(a,b,x0,e,m) %定义jacobi 函数，其中：a,b 为线性方程组Ax b =中的矩阵和右端向量；x0为初始值； %e 和m 分别为人为设定的精度和预计迭代次数；运行结果y 为迭代的结果和所有中间值组成的 %矩阵 y=0; %对y 初始化 d=diag(diag(a)); %按雅可比迭代标准形形式取主对角元素作为矩阵D u=-triu(a,1); %取上三角矩阵u l=-tril(a,-1); %取下三角矩阵l bj=d^-1*(l+u); fj=d^-1*b; x=[x0,zeros(20,m-1)]; %初始化x,其中x1=x0,即初始值 for k=1:m %人为规定迭代次数，防止不收敛迭代导致死循环 x(:,k+1)=bj*x(:,k)+fj; %jacobi 迭代 if norm(x(:,k+1)-x(:,k),inf)线性代数方程组的直接解法_赖志柱
matlab线性代数实验
一维抛物线型方程数值解法(1)(附图及matlab程序)
线性代数方程组数值解法及MATLAB实现综述
第二章非线性方程(组)的数值解法的matlab程序
线性代数方程组求解
双曲方程基于matlab的数值解法
线性代数及matlab英汉对照
线性代数方程组求解
数学实验 5：线性代数方程组的数值解法
相关主题
文本预览