线性方程组数值解法
第二章 线性方程组的数值解法
第二章 线性方程组的数值解法在科技、工程技术、社会经济等各个领域中很多问题常常归结到求解线性方程组。
例如电学中的网络问题，样条函数问题，构造求解微分方程的差分格式和工程力学中用有限元方法解连续介质力学问题，以及经济学中求解投入产出模型等都导致求解线性方程组。
n 阶线性方程组的一般形式为⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++nn nn n n n n n n b x a x a x a b x a x a x a b x a x a x a L K K K K L L 22112222212********* （1.1） 其矩阵形式为b Ax = （1.2） 其中⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=n n nn n n n n b b b b x x x x a a a a a a a a a A M M L K K K K L L 2121212222111211),,2,1,(n j i a ij L =，),,2,1(n i b i L =均为实数，i b 不全为0，且A 为非奇异。
关于线性方程组的数值解法一般分为两类：1．直接法 就是不考虑计算机过程中的舍入误差时，经有限次的四则运算得到方程组准确解的方法。
而实际中由于计算机字长的限制，舍入误差的存在和影响，这种算法也只能求得线性方程组的近似解。
本章将阐述这类算法中最基本的消去法及其某些变形。
这些方法主要用于求解低阶稠密系数矩阵方程组。
2．迭代法 从某个解的近似值出发，通过构造一个无穷序列，用某种极限过程去逐步逼近线性方程组的精确解的方法。
本章主要介绍迭代法与迭代法。
迭代法是解大型稀疏矩阵（矩阵阶数高而且零元素较多）的线性方程组的重要方法。
§1 高斯)(Gauss 消去法1.1 Gauss 消去法Gauss 消去法是将线性方程组化成等价的三角形方程组求解。
首先举例说明Gauss消去法的基本思想和过程。
第4章线性方程组的数值解法
。
（3）消元过程： 对 k 1,2,, n 1 ,计算：
k ai(k ) ( a k kk)
ci k
( ai( kj 1) ai( kj ) ci k a kkj)
( i , j k 1 , k 2 , , n )
bi( k 1) bi( k ) ci k bk( k )
a n1 x1 a n 2 x2 ann xn bn
(4.1)
( i, j 1, 2 ,3 , , n )
其中，a i j 、bi 为常数， x i 为待求的未知量
用矩阵形式表示是
其中：
Ax b
a11 a 21 A a n 1
a12 a 22 an 2
( i , j k 1 , , n )
bi( k 1) bi( k ) ci k bk( k )
( 只要设 a k kk) 0 就可以继续进行消元，直到经过 n 1 次消元后， 将线性方程组(4.1)化为(4.2)所示上三角方程组（以上计算过程称为 消元过程）。 (1 (1 a11) a12) a1(1) x1 b1(1) n ( (2 (2 ( a 22) a 23) a 22 ) x 2 b22 ) n (k ) (k ) (k ) (4.2) a kk a kn x k bk (n) (n) a n n x n bn
但其绝对值很小时，用它作除数时，根据数值运算中“用绝对值很小
的数做除数，舍入误差会增大，而且严重影响计算结果的精度”的原 则，这种方法在一定程度上也具有局限性。为了克服这一局限性，产
数值计算08-线性方程组数值解法(优选.)
0
(k=1,2,…,n) ,则可通过高斯消元法求出Ax=b 的解。
引理
A的主元素
a(k) kk
0
(k=1,2,…,n) 的充要条件
是矩阵A的各阶顺序主子式不为零，即
a11
a1k
D1 a11 0 Dk
0, k 2, 3, , n
ak1
akk
定理2 Ax=b 可用高 斯消元法求解的充分必要条件是： 系数矩阵 A 的各阶顺序主子式均不为零。
Page 5
线性代数方程组的计算机解法常用方法：
直接法 迭代法
消去法 矩阵三角分解法
Page 6
直接法：经过有限步算术运算，可求得方程组
的精确解的方法（若在计算过程中没有舍入误差）
迭代法：用某种极限过程去逐步逼近线性方程
组精确解的方法 迭代法具有占存储单元少，程序设计简单，原
始系数矩阵在迭代过程中不变等优点，但存在收 敛性及收敛速度等问题
a(k) ik
a(k) kk
aijk
mik
a
k
kj
bik1 bik mikbkk
xn
bnn annn
bii
n
a
i
ij
x
j
,
xi
ji1
aiii
i, j k 1, k 2,, n
i n 1,,2,1
高斯消元法的条件
Page 20
定理1
如果在消元过程中A的主元素
a(k) kk
即:
a111
a112 a222
a11n a22n
x1 x2
bb1212
an22
an2n
xn
bn2
其中：
实验五(线性方程组的数值解法和非线性方程求解)
1大学数学实验 实验报告 | 2014/4/5一、 实验目的1、学习用Matlab 软件数值求解线性代数方程组，对迭代法的收敛性和解的稳定性作初步分析；2、通过实例学习用线性代数方程组解决简化问题。
二、 实验内容项目一：种群的繁殖与稳定收获：种群的数量因繁殖而增加，因自然死亡而减少，对于人工饲养的种群（比如家畜）而言，为了保证稳定的收获，各个年龄的种群数量应维持不变。
种群因雌性个体的繁殖而改变，为方便起见以下种群数量均指其中的雌性。
种群年龄记作k=1,2,…,n ,当年年龄k 的种群数量记作x k ，繁殖率记作b k （每个雌性个体1年的繁殖的数量），自然存活率记作s k （s k =1−d k ，d k 为1年的死亡率），收获量记作ℎk ，则来年年龄k 的种群数量x ̌k 应该为x ̌k =∑b k n k=1x k , x ̌k+1=s k x k −ℎk , (k=1,2,…,n -1)。
要求各个年龄的种群数量每年维持不变就是要求使得x ̌k =x k , (k=1,2,…,n -1).(1) 如果b k , s k 已知，给定收获量ℎk ，建立求各个年龄的稳定种群数量x k 的模型（用矩阵、向量表示）.(2) 设n =5,b 1=b 2=b 5=0,b 3=5,b 4=3,s 1=s 4=0.4,s 2=s 3=0.6,如要求ℎ1~ℎ5为500,400,200,100,100，求x 1~x 5.(3) 要使ℎ1~ℎ5均为500，如何达到？问题分析：该问题属于简单的种群数量增长模型，在一定的条件（存活率，繁殖率等）下为使各年龄阶段的种群数量保持不变，各个年龄段的种群数量将会满足一定的要求，只要找到种群数量与各个参量之间的关系，建立起种群数量恒定的方程就可以求解出各年龄阶段的种群数量。
模型建立：根据题目中的信息，令x ̌k =x k ，得到方程组如下：{x ̌1=∑b k nk=1x k =x 1x ̌k+1=s k x k −ℎk =x k+1整理得到：{−x 1∑b k nk=1x k =0−x k+1+s k x k =ℎk2 大学数学实验 实验报告 | 2014/4/52写成系数矩阵的形式如下：A =[b 1−1b 2b 3s 1−100s 2−1…b n−1b n0000⋮⋱⋮000000000⋯00−10s n−1−1]令h =[0, ℎ1,ℎ2,ℎ3,…,ℎn−2,ℎn−1]Tx =[x n , x n−1,…,x 1]T则方程组化为矩阵形式：Ax =h ，即为所求模型。
第2章 线性方程组的数值解法
第2章 线性方程组的数值解法2.1 引言在自然科学研究和工程技术的应用中，许多问题的解决，诸如非线性问题线性化、求微分方程的数值解最终都归结为线性方程组的求解问题. 我们在后面章节中的样条插值、曲线拟合、数值代数等，也需要求解线性方程组。
一般地，设n 阶线性方程组(linear system of equations of order n )为11112211211222221122,,,n n n n n n nn n n a x a x a x b a x a x a x b a x a x a x b +++=⎧⎪+++=⎪⎨⎪⎪+++=⎩ (2.1.1) 表示成矩阵形式=Ax b , (2.1.2)其中()111212122212n n ij n nn n nn a a a a a a a a a a ⨯⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎣⎦ A ，12n x x x ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦ x ，12n b b b ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦b , (2.1.3) A 为系数矩阵(coefficient matrix).目前在计算机上经常使用的、简单有效的线性方程组的数值解法大致分为两类：直接法(direct method)和迭代法(iterative method). 其中直接法适用于以稠密矩阵为系数矩阵的中低阶线性方程组，而迭代法主要用于求解以稀疏矩阵为系数矩阵的高阶线性方程组。
本章首先介绍解线性方程组的两种常用的直接法：Gauss 消去法与矩阵三角分解法；然后介绍解线性方程组的三种常用的迭代法：Jacobi 迭代法、Gauss-Seidel 迭代法、超松弛法(SOR 法)，并讨论它们的收敛性。
最后，讨论了线性方程组的性态。
2.2 Gauss 消去法Gauss 消去法（Gaussian elimination method ）的基本思想是使用初等行变换将方程组转化为一个同解的上三角形方程组，再通过回代，求出该三角形方程组的解.2.2.1 Gauss 消去法Gauss 消去法包括消元和回代两个过程. 下面先举例说明Gauss 消去法求解线性方程组的主要过程.例2.2.1 求解线性方程组123123123471,2581,3611 1.x x x x x x x x x ++=⎧⎪++=⎨⎪++=⎩ 解 将该线性方程组写成增广矩阵(augmented matrix)的形式1471258136111⎡⎤⎢⎥⎢⎥⎢⎥⎣⎦用Gauss 消去法求解过程如下：1．消元过程12213323323214711471147125810361036136111061020020r r r r r r r r r -+→-+→-+→⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥−−−−→---−−−−→---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥---⎣⎦⎣⎦⎣⎦,从而原方程组等价地变为上三角形方程组123233471,361,20.x x x x x x ++=⎧⎪--=-⎨⎪=⎩2．回代过程从第3 个方程解出30x =，将其代入第2 个方程得()2216/31/3x x =--+=，再将30x =及21/3x =回代到第1个方程，解出1231471/3x x x =--=-. 从而得到原方程组的解123113,0.x x x =-==对于一般线性方程组(2.1.1)，使用Gauss 消去法求解分为以下两步：1．消元过程为方便起见，记()(0)(0),ijn na ⨯==A A ()T(0)(0)(0)(0)1,12,1,1,,,n n n n a a a +++== b b ，则方程组(2.1.1)为()()()()()()()()()()()()000011112211,1000021122222,100001122,1,,.n n n n n n n n nn n n n a x a x a x a a x a x a x a a x a x a x a +++⎧+++=⎪+++=⎪⎨⎪⎪+++=⎩ (2.2.1) 第1次消元：若()0110a ≠，对方程组(2.2.1) 执行初等行变换11i i i r l r r -→, 2,3,,i n = ,得第1个导出方程组——————————————————————————高斯 (Carl Friedrich Gauss 1777年4月30 日 – 1855年2月23 ) 是德国数学家、天文学家，在许多科学领域都做出了杰出的贡献，他为现代数论、微分几何(曲面论)、误差理论等许多数学分支奠定了基础. 他的数学研究以简明、严谨、完美而著称于世. 他在数学上与阿基米德、牛顿和欧拉齐名，被称为“数学王子”，被公认为有史以来最伟大的数学家之一.()()()()()()()()()()000011112211,111122222,111122,1,,,n nn n n n n nn n n n a x a x a x a a x a x a a x a x a +++⎧+++=⎪⎪++=⎨⎪⎪++=⎩(2.2.2)其中()()1111/i i l a a =，()()()10011,2,3,,;2,3,, 1.ij ij i j a a l a i n j n =-==+第2次消元：若()1220a ≠，对方程组(2.2.2)执行初等行变换22,i i i r l r r -→ 3,4,,i n = ,得第2个导出方程组()()()()()()()()()()()()()()()0000011112213311,1111122223322,122233333,122233,1,,,,n nn n n n n nn n nn nn n a x a x a x a x a a x a x a x a a x a x a a x a x a ++++⎧++++=⎪⎪+++=⎪⎪++=⎨⎪⎪⎪++=⎪⎩(2.2.3)其中()()112222/i i l a a =，()()()21122,3,4,,;3,4,, 1.ij ij i j a a l a i n j n =-==+第k 次消元：若()10k kka -≠，对第1k -个导出方程组执行初等行变换i ik k i r l r r -→，1,2,,i k k n =++ , 得第k 个导出方程组()()()()()()()()()()()()()()()000001111221,1111,110112221,1122,11,111,1,1,11,1,,,,k k n nn k k n n n k k k k k k k n nk n k k k n k k nn n n n a x a x a x a x a a x a x a x a a x a x a a x a x a +++++++++++++++⎧+++++=⎪⎪++++=⎪⎪⎨++=⎪⎪⎪⎪++=⎩(2.2.4)其中()()11k k ik ikkkl a a --=，()()()11,k k k ij ij ik kj a a l a --=- 1,,;1,, 1.i k n j k n =+=++ 重复上述过程1n -次，得到第1n -个导出方程组()()()()()()()()()()()()()()0000011112213311,1111122223322,122233333,111,1,,,.n n n n n n n nn n n nn nn n a x a x a x a x a a x a x a x a a x a x a a x a +++--+⎧++++=⎪⎪+++=⎪⎪++=⎨⎪⎪⎪=⎪⎩(2.2.5)其中()()()()()1111,,1,2,,1;1,,;1,, 1.k k kk k ik ikkkij ijik kjl a a a a l a k n i k n j k n ----==-=-=+=++ (2.2.6)这样，通过消元过程就将方程组(2.1.1)化成了等价的上三角形方程组(2.2.5).2．回代过程回代过程就是求上三角形方程组(2.2.5)的解. 若()10n nna -≠，则从最后一个方程开始，先求出()()11,1,/n n n n n n n x a a --+=，再由第1n -个方程解出1n x -，依此类推可解出221,,,n x x x - . 一般(2.2.7)定义 2.2.1 由式(2.2.2)-(2.2.7)确定的求解线性方程组的算法称为Gauss 消去法(Gaussian elimination method)，包括消元(elimination)和回代(backward substitution)两个过程。
数值计算方法第3章解线性方程组的数值解法1
,i
2 ,3 ,...,
n
a
(1 11
)
A（ 1） A ( 2 )
a (1) 11
a (2) 22
...... ......
......
a (2) n2
......
a a
(1) 1n
(2) 2n
a
(2 nn
)
b (1)
b (2)
[
b
( 1
1
)
b (2) 2
a(k) kk
...
a(k) kn
... ... ...
...
...
a(n) nn
b1(1) b2(2)
...
bk(k)
...
bn(n)
21
高斯顺序消去法
也就是对于方程组AX=b系数矩阵做：
ai(jkl1i)k
a(k) ik
a(k) ij
/
a(k) kk
3)顺序消元
31
高斯列主元消去法
第k步
从A ( k ) 的第
k
列
a (k) kk
，a (k) k 1k
，...a
(k) nk
中选取绝对值
最大项，记录所在行，即
|a(k) ikk
|m kina|axi(kk)
|
记 lik
若 l k 交换第k行与l行的所有对应元素，再 进行顺序消元。
32
其中， lii 0, i 1,2,..., n
（1）
10
高斯顺序消元法
线性方程组的四种数值解法
线性方程组的四种数值解法(电子科技大学物理电子学院，四川 成都 610054)摘要：本文介绍了四种求解线性方程组的数值解法: 雅克比迭代法、高斯赛德尔迭代法、高斯消去法和改进的平方根法的基本原理和算法流程，通过求解具体方程，对四种求解方法进行了对比。
对于雅克比迭代法和高斯赛德尔迭代法，研究了两种算法对求解同一方程组的迭代效率差异，结果表明高斯赛德尔迭代法达到同样精度所需迭代次数较少。
对于高斯消去法，通过选择列主元的方法提高算法的准确度，计算结果表明高斯消去法计算精确，且运算复杂度也不是很高。
对于改进的平方根法，其运算复杂度低，但对于给定的方程组有着严苛的要求。
关键词：雅克比迭代法；高斯赛德尔迭代法；高斯消去法；改进的平方根法；线性方程组引言线性方程组的求解在日常生活和科研中有着极其重要的应用，但在实际运算中，当矩阵的维数较高时，用初等方法求解的计算复杂度随维数的增长非常快，因此，用数值方法求解线性方程组的重要性便显现出来。
经典的求解线性方程组的方法一般分为两类：直接法和迭代法。
前者例如高斯消去法,改进的平方根法等，后者的例子包括雅克比迭代法，高斯赛德尔迭代法等。
这些方法的计算复杂度在可以接受的范围内，因此被广泛采用。
一般来说，直接法对于阶数比较低的方程组比较有效；而后者对于比较大的方程组更有效。
在实际计算中，几十万甚至几百万个未知数的方程组并不少见。
在这些情况下，迭代法有无可比拟的优势。
另外，使用迭代法可以根据不同的精度要求选择终止时间，因此比较灵活。
在问题特别大的时候，计算机内存可能无法容纳被操作的矩阵，这给直接法带来很大的挑战。
而对于迭代法，则可以将矩阵的某一部分读入内存进行操作，然后再操作另外部分。
本文使用上述四种算法求解对应的方程组，验证各种算法的精确度和计算速度。
1 算法介绍1.1 雅克比迭代法 1.1.1 算法理论设线性方程组(1)b Ax的系数矩阵A 可逆且主对角元素 均不为零,令并将A 分解成(2)从而(1)可写成令其中. (3)以B 1为迭代矩阵的迭代法(公式)(4)称为雅克比(Jacobi)迭代法(公式),用向量的分量来表示,(4)为(5)其中为初始向量.1.1.2 算法描述 1给定迭代初始向量X 0以及误差要求delta 2根据雅克比迭代公式计算出下一组向量3判断X 是否满足误差要求，即||X k+1 – X k || < delta4若误差满足要求，则停止迭代返回结果；若否，则返回第二步进行下一轮迭代1.2 高斯赛德尔迭代法nna ,...,a ,a 2211()nna ,...,a ,a diag D 2211=()D D A A +-=()b x A D Dx +-=11f x B x +=b D f ,A D I B 1111--=-=()()111f x B x k k +=+⎩⎨⎧[],...,,k ,n ,...,i x a ba xnij j )k (j j i iii)k (i21021111==∑-=≠=+()()()()()Tn x ,...x ,x x 002010=1.2.1 算法理论由雅克比迭代公式可知,在迭代的每一步计算过程中是用的全部分量来计算的所有分量,显然在计算第i 个分量时,已经计算出的最新分量没有被利用，从直观上看,最新计算出的分量可能比旧的分量要好些.因此，对这些最新计算出来的第次近似的分量加以利用,就得到所谓解方程组的高斯—塞德尔（Gauss-Seidel ）迭代法.把矩阵A 分解成(6)其中,分别为的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成即其中(7)以为迭代矩阵构成的迭代法(公式)(8)称为高斯—塞德尔迭代法(公式),用变量表示的形式为(9)1.2.2 算法描述 1给定迭代初始向量X 0以及误差要求delta2根据高斯赛德尔迭代公式计算出下一组向量()k x ()1+k x ()1+k ix ()()1111+-+k i k x ,...,x 1+k()1+k x()1+k jx U L D A --=()nna ,...,a ,a diag D 2211=U ,L --A ()b Ux x L D +=-22f x B x +=()()b L D f ,U L D B 1212---=-=2B ()()221f x B x k k +=+⎩⎨⎧[],...,,k ,n ,,i x a x a b a xi j n i j )k (j ij )k (j ij i ii)k (i21021111111==∑∑--=-=+=++3判断X是否满足误差要求，即||X k+1– X k|| < delta4若误差满足要求，则停止迭代返回结果；若否，则返回第二步进行下一轮迭代1.3 高斯消去法1.3.1 算法理论下面三种变换称为初等行变换：1.对调两行；2.以数k≠0乘某一行中的所有元素；3.把某一行所有元素的k倍加到另一行对应的元素上去。
线性代数方程组的数值解法_百度文库
线性代数方程组的数值解法【实验目的】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后的解。
第二章 线性方程组的数值解法
0
a(k1) n,k1
an (kn1)
bn (k1)
12/114 §2.1 Gauss消去法
1. 消去过程
(2) 第k次消元。
ai(k j1)ai(k j)a ak i((k k k ))kai(k j) bi(k1)bi(k)a ak i((k k k ))kbk (k)
jk1,k2,.n ..— — , 减 减bk (k 去 去 k )的 行 a ak i((k k k ))a a 第 k倍 k i((的 k k k ))k倍 ik1,k2,.n ..
13/114 §2.1 Gauss消去法
1. 消去过程
பைடு நூலகம்
(3) 当k = n – 1时得
a1(11) [A(n) |b(n)]
a(1) 12
a(2) 22
a(1) 1n
a2(2n)
bb12((12))
a(n) nn
bn(n)
完成第n – 1次消元后得到与原方程组等价的三角形方程
组
A(n)x = b(n)
令 lik
a(k) ik
a(k) kk
for (i = k+1 to n) A(i，k) = A(i，k)/ A(k，k) for (j = k+1 to n+1) A(i,j) = A(i,j) - A(i,k)* A(k,j)
a i ( k j 1 ) a i ( k ) j l ia k k ( k ) ,( j j k 1 , ,n )
求解三角形方程组A(n)x = b(n)，
an (n)nxn bn (n)
得到求解公式
xn
abnn((nnn))
bi(i)
数值分析线性方程组求解的数值方法
记
Lk
1
k 1, 2 , , n 1
m ik
a ik
(k ) (k )
， ( i , j k 1, k 2 , , n )
用 4 位浮点数计算精确解， 然后舍入到 4 位有效数字， 解出原方程的组的解为： （MATLAB 求解程序如下） A=[0.001 2.0 3.0;-1 3.712 4.623;-2 1.072 5.643]; >> b=[1 2 3]'; >> X=A\b
x1 0.4904 ， x2 0.05104 ， x3 0.3675 ，
a r i m r k u r i m rk u k i u ri ， 有 ：
k 1 r 1 k 1
n
r1
u r i a r i m r k u k i ， i r , r 1, , n ） （
k 1 r 1
a i r m ik u k r
比较得知，用高斯消元法求解结果误差较大，不能作 为原方程组的近似解，其原因就是在消元过程中使用了小 主元 0.001，至使方程组的系数数量级增加，使误差扩散。
避免方法：高斯主元消元法
高斯主元消去法的MATLAB实现
程 序 8- 2 用 高 斯 列 主 元 消 去 法 求 解 线 性 方 程 组 AX b ， 首 先 将 矩 阵 A 化 为上三角矩阵，再执行回代过程。
用回代法求解上三角线性方程组AX=B，其中A为非奇异。
function X=backsub(A,b) %A是一个n阶上三角非奇异阵。 %b是一个n维向量。 %X是线性方程组AX=b的解。 n=length(b); X=zeros(n,1); X(n)=b(n)/A(n,n); for k=n-1:-1:1 X(k)=(b(k)-A(k,k+1:n)*X(k+1:n))/A(k,k); End
线性方程组的数值解法LU分解法市公开课金奖市赛课一等奖课件
推论 2 D 并入 L，则
A (L D)R L U
此时， L 是下三角阵， U 是单位上三角阵，称之为
Crout 分解.
第8页
矩阵分解理论
推论 3 如果 A AT ,则A LDLT
其中，L 是单位下三角阵，D 是对角阵.
由a2 j l21u1 j 1 u2 j 得u2 j a2 j l2iu1 j ( j 2,3,..., n)；
再由ai2 li1u12 li2u22
得li 2
ai 2
li1u12 u22
(i 3,4,..., n)。
第14页
Doolittle分解
第k步时：计算ukk , ukk1,n j k
i 1
yi bi lij y j i 1,2,..., n j 1 n
xi ( yi uij x j ) / uii i n, n 1,...1 j i 1
x 可获解 (x1, x2 ,..., xn )T。
第18页
例题
例1.试用Doolittle分解求解方程组.
2 5 6 x1 10
3.5 LU分解法 我们知道对矩阵进行一次初等变换，就相
称于用相应初等矩阵去左乘本来矩阵。因 此我们从这个观点来考察Gauss消元法并 用矩阵乘法来表示，即可得到求解线性方 程组另一个直接法：矩阵三角分解。
第1页
高斯消元过程矩阵表示
第1步等价于
:
a (1) 11
0时，将a（211）, a（311）,..., a（n11）消零, 令li1
d1
令
U1
d2
d
n
第29页
《应用数值分析》课件数值分析5.3线性方程组的数值解法
Step k：设ak(kk) ，0计算因子
mik
a(k) ik
/
a(k kk
)
(i k 1, ..., n)
且计算
a ( k 1) ij
b( k 1) i
a(k) ij
m
ik
a
(k kj
)
b(k ) i
mik bk(k )
(i, j k 1, ..., n)
n
bi (bi
aij * b j ) / aii
j i 1
2024/11/23
线性方程组的直接解法
11
计算量 /* Amount of Computation */
由于计算机中乘除 /* multiplications / divisions */ 运算的时 间远远超过加减 /* additions / subtractions */ 运算的时间，故 估计某种算法的运算量时，往往只估计乘除的次数，而且通 常以乘除次数的最高次幂为运算量的数量级。 (n k) 次
(k)
kk
k ,k1
0
a ( k 1) k 1,k 1
a(1) 1n
a(2) 2n
a(k) kn
a ( k 1) k 1,n
0
a ( k 1) n,k 1
a ( k 1) nn
第 6 章 不动点理论及应用 第 1 页 共 1 页
b(1) 1
b(2) 2
b( k ) k
b( k 1) k 1
b( k 1) n
xn
b(n) n
/
a(n) nn
n
b( i ) i
a
(i ij
计算方法线性方程组数值解法
d
2
a3b3c3
x3
d3
an
1bn1cn
1
xn
1
d
n
1
anbn xn dn
其系数矩阵为三对角形，元素满足以下条件：
|b1|>|c1|>0
|bi|≥|ai|+|ci|，且aici≠0 i=2，3，……n-1; |bn|≥|an|>0。
可以采用追赶法求解
4
线性代数方面的计算方法就是研究求解线 性方程组的一些数值解法与研究计算矩阵 的特征值及特征向量的数值方法。
5
设有线性方程组
a11x1 a12x2 a1nxn b1 a21x1a22x2a2nxnb2 an1x1 an2x2 annxn bn
式中，aij，bi为已知常数，xi为待求的未知量。记
u
2
2
u 2 n
u n 1,n 1u n 1,n
u n n
10
若uii≠0（i=1,2,……n）,则由下至上依次回代得
xn yn / unn
xn1 ( yn1 xi yi
un1,n xn ) / un1,n1
n
uij x j ) / uii
0
a
( 2
2 2
)
a
( 2
2) ,k 1
a
( 2
2) ,k
a
( 2
2) ,n
a
( 2
2) ,n 1
0 A(k)
0 0
a
( k
k) ,k
a
( k
k) ,k 1
a
k
k ,n
a
( k
k) 1,n
1
线性方程组的数值解法-安振华-2012011837
实验5：线性方程组的数值解法化学工程系分2 安振华2012011837【实验目的】1、掌握线性方程组的常用数值解法，包括高斯消去法、LU分解法以及校正法。
2、体验数值计算的时间复杂度和计算规模的关系。
3、加深对数值计算误差的理解。
4、学习使用迭代法等算法，求解非线性方程。
5、学习如何使用MATLAB解非线性方程组和方程组。
【实验容】【实验五：习题9】种群的繁殖与稳定收获：种群的数量因繁殖而增加，因自然死亡而减少，对于人工饲养的种群（比如家畜）而言，为了保证稳定的收获，各个年龄的种群数量应保持不变，种群因雌性个体的繁殖而改变，为方便起见以下种群数量均指其中的雌性。
种群年龄记作k=1,2,…,n，当年年龄k的种群数量记作x k，繁殖率记作b k（每个雌性个体在1年繁殖的数量），自然存活率记作s k（s k=1-d k，d k为1年的死亡率），收获量记作h k，则来年年龄k的种群数量k x应为：111,(1,2,,1)n k k k k k k k x b x x s x h k n +===-=⋅⋅⋅-∑要求各个年龄的种群数量每年维持不变就是要使(1,2,,)k k x x k n ==⋅⋅⋅（1） 若b k ，s k 已知，给定收获量h k ，建立求各年龄的稳定种群数量x k 的模型（用矩阵向量表示）（2） 设n=5，b 1=b 2=b 5=0，b 3=5，b 4=3，s 1=s 4=0.4，s 2=s 3=0.6，如果要求h 1~h 5为500,400,200,100,100,求x 1~x 5 （3） 要使h 1~h 5均为500，如何达到？ 【分析】为方便起见以下种群数量均指其中的雌性。
我们并且有以下的假设：(1)雌性个体的繁殖率和存活率在特定的时间是不变的。
(2)人工饲养的种群在质量和数量上是不受外界环境和资源的限制的。
(3)模型中不考虑人为的或是自然的灾害所造成的种群数量、繁殖率和存活率的变动。
数学建模线性方程组的数值解法
直接法 经过有限次算术运算求出精确解（实际上 由 于 有 舍 入 误 差 只 能 得 到 近 似 解 ） ---- 高 斯 （Gauss）消元法及与它密切相关的矩阵LU分解 迭代法 从初始解出发，根据设计好的步骤用逐次 求出的近似解逼近精确解 ---- 雅可比（Jacobi） 迭代法和高斯—塞德尔（Gauss—Seidel）迭代法
(k )
0.1x1
( k 1)
0.3x2
( k 1)
1.4
Gauss-Seideil迭代公式 Dx ( k 1) Lx ( k 1) Ux ( k ) b
用它作除数会导致舍入误 差的很大增加 解决 办法 选
(k ) aik
(i k , n) 最大的一个(列主元)
将列主元所在行与第k行交换后, 再按上面的 高斯消元法进行下去，称为列主元消元法。
直接法 - 高斯消元法的矩阵表示
高斯消元法的第一次消元
a11 x1 a12 x2 a1n xn b1 a21 x1 a22 x2 a2 n xn b2 an1 x1 an 2 x2 ann xn bn
数值解法(迭代解法)的收敛性
实验5的主要内容
1. 两类数值解法: 直接方法；迭代方法
2. 超定线性方程组的最小二乘解 3. 线性方程组数值解法的MATLAB实现 4. 实际问题中方程组的数值解
线性方程组的一般形式、两类解法
a11 x1 a12 x 2 a1n x n b1 a 21 x1 a 22 x 2 a 2 n x n b2 a n1 x1 a n 2 x 2 a nn x n bn
大学数学实验
Mathematical Experiments 实验5 线性代数方程组的数值解法
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
.计算法实验题目：班级：学号：：目录计算法实验 (1)1 实验目的 (3)2 实验步骤 (3)2.1环境配置： (3)2.2添加头文件 (3)2.3主要模块 (3)3 代码 (3)3.1主程序部分 (3)3.2多项式程部分 (3)3.3核心算法部分 (3)3.4数据结构部分 (3)4运行结果 (3)4.1列主元高斯消去法运行结果 (3)4.2LU三角分解法运行结果 (3)4.3雅克比迭代法运行结果 (3)边界情况调试 (3)5总结 (3)输入输出 (3)列主元高斯消元法 (3)雅克比迭代法 (3)6参考资料 (3)1 实验目的1.通过编程加深对列主元高斯消去法、LU三角分解法和雅克比迭代法等求解多项式程法的理解2.观察上述三种法的计算稳定性和求解精度并比较各种法利弊2 实验步骤2.1环境配置：VS2013，C++控制台程序2.2添加头文件#include "stdio.h"#include "stdlib.h"#include "stdafx.h"#include<iostream>2.3主要模块程序一共分成三层，最底层是数据结构部分，负责存储数据，第二层是交互部分，即多项式程部分，负责输入输出获得数据，最上层是核心的算法部分，负责处理已获得的数据。
具体功能如下：●数据结构部分数据结构部分是整个程序的最底层，负责存储部分。
因数组作为数据元素插入和删除操作较少，而顺序表空间利用率大且查看便，故此程序选用二维顺序表保存系数。
数据结构文件中写的是有关其的所有基本操作以供其他文件调用。
●多项式程部分多项式程部分是程序的第二层，容是有关程组的所有函数、构建程、输出程等等，同时在此文件中获得程系数并储存，同时此文件还负责显示菜单部分。
算法部分此文件负责核心算法，处于整个程序最上层部分，负责列主元高斯消去法、LU三角分解法和雅克比迭代法的具体实现过程。
通过调用程文件的函数，将获得的数据进行处理运算，可以得到结果返回给程主函数和输出的第二层。
总结：主函数负责获取程系数并显示，算法和程作为后台程序，顺序表作为存储手段。
3 代码3.1主程序部分// Solutionoflinearquations.cpp : 定义控制台应用程序的入口点。
//#include "stdio.h"#include "stdlib.h"#include "stdafx.h"#include "squencelist.h"#include "equation.h"#include "algorithm.h"#include<iostream>int _tmain(int argc, _TCHAR* argv[])while (Exflag){GetEquation();ShowMenu();}return 0;}3.2多项式程部分 程部分头文件#ifndef _EQUATION_H #define _EQUATION_H#include "stdio.h"#include "stdlib.h"#include "squencelist.h"extern int Xnumbers; extern int Fnumber; extern int Exflag; extern datacoa *A;void GetEquation(void);void ShowMenu(void);void printfunction(datacoa *A); void printresx(datacoa *A); void Tip(void);#endif程部分CPP文件#include "stdafx.h"#include "equation.h"#include "math.h"#include "algorithm.h"#include "stdio.h"#include "stdlib.h"#include<iostream>#include <iomanip>using namespace std;//全局变量int Xnumbers = 0;int Fnumber = 0;int Exflag = 1;datacoa *A;////////////////////////多项式函数系数/////////////////////////void GetEquation(void){int i, j,flag=1;float x;A = InitStruct();while (flag){cout << "程未知量和解总个数：" << endl;cin >> Xnumbers;cout << "程个数：" << endl;cin >> Fnumber;cout << "输入程系数,输入00结束（如输入2 1 5 或者2 1 5 3 4 1 00 " << endl;cout << " 3 4 1 00 ）:" << endl;cin >> x;while (x != 00){for (i = 1; i <= Fnumber; i++){for (j = 1; j <= Xnumbers; j++){if (!Insert(A, x, i, j))exit(0);cin >> x;}}}j = 1;printfunction(A);if (Xnumbers == Fnumber+1)flag = 0;else{cout << "程可能无解" << endl;A->m = 0;A->n = 0;}}}//////////////////////////显示交互/////////////////////////////void ShowMenu(void){int x;cout << "选择求解程根的法：" << endl;cout << "1.列主元高斯消去法" << endl;cout << "2.三角分解法" << endl;cout << "3.迭代法" << endl;cin >> x;switch (x){case 1:ColumnGaussmethod(A);Tip();break;case 2:LUmethod(A);Tip();break;case 3:ISODATAmethod(A);Tip();break;default:break;}}////////////////////////打印输出函数/////////////////////////// void printfunction(datacoa *A){int i,j;cout << "矩阵=" << endl;for (i = 1; i <= A->m; i++){for (j = 1; j <= A->n; j++){cout <<setw(12)<< A->data[i][j];}cout << endl;}}////////////////////////打印结果/////////////////////////// void printresx(datacoa *A){int i;for (i = 1; i <= A->m; i++){cout << "X" << i << "=" << setw(12) << A->data[i][Xnumbers+1];cout << endl;}}////////////////////////返回提示///////////////////////////void Tip(void){int flag;cout << "输入000退出,其余返回：" << endl;cin >> flag;if (flag == 000)Exflag = 0;}3.3核心算法部分算法部分头文件#ifndef _ALGORITHM_H#define _ALGORITHM_H#include "stdio.h"#include "stdlib.h"#include "stdafx.h"int Judge(double x1, double x0, double e); void ColumnGaussmethod(datacoa *A); void LUmethod(datacoa *A);void ISODATAmethod(datacoa *A);#endif算法部分CPP文件#include "algorithm.h"#include "stdafx.h"#include "squencelist.h"#include "equation.h"#define max 100//////////////////////////误差判别////////////////////////////inline int Judge(double x1, double x0, double e){if (e == 0)return 0;if (x1 == 0){if (fabs(x0) < e)return 1;else return 0;}if (x0 == 0){if (fabs(x1) < e)return 1;else return 0;}if (fabs(x1 - x0) < e)return 1;else return 0;}//////////////////////////列主元高斯消元法//////////////////////////// void ColumnGaussmethod(datacoa *A){cout << "///////////////////列主元高斯消元法///////////////////" << endl; int i, j, i2, flagc, k, j2;float temp, res;for (i = 1; i < Fnumber; i++){flagc = i;for (i2 = i + 1; i2 <= Fnumber; i2++)if ((fabs(A->data[i2][i]))>(fabs(A->data[flagc][i])))flagc = i2;if (flagc != i)for (k = i; k <= Xnumbers; k++){temp = A->data[i][k];A->data[i][k] = A->data[flagc][k];A->data[flagc][k] = temp;}for (i2 = i + 1; i2 <= Fnumber; i2++){temp = A->data[i2][i]/A->data[i][i];for (j2 = i; j2 <= Xnumbers; j2++)A->data[i2][j2] = A->data[i2][j2] - temp*A->data[i][j2];}}for (i = Fnumber; i >= 1; i--){for (j = Fnumber; j >= i + 1; j--)A->data[i][Xnumbers] = A->data[i][Xnumbers] - A->data[i][j]*A->data[j][Xnumbers + 1];res = A->data[i][Xnumbers] / A->data[i][i];Insert(A, res, i, Xnumbers+1);}printresx(A);}//////////////////////////直接三角分解法////////////////////////////void LUmethod(datacoa *A){datacoa *L, *U;int i, j, k,q;float temp=0.0;L = InitStruct();U = InitStruct();for (i = 1; i <= Fnumber; i++){for (j = 1; j <= Xnumbers; j++){Insert(L, 0, i, j);Insert(U, 0, i, j);}}for (j = 1; j <= Xnumbers; j++)Insert(U, A->data[1][j], 1, j);for (i = 2; i <= Fnumber; i++)Insert(L, A->data[i][1] / A->data[1][1], i, 1);for (k = 2; k <= Fnumber; k++){for (j = k; j <= Xnumbers; j++){for (q = 1; q < k; q++)temp = temp + L->data[k][q] * U->data[q][j];Insert(U, A->data[k][j] - temp, k, j);temp = 0;}for (i = k + 1; i <= Fnumber; i++){for (q = 1; q < k; q++)temp = temp + L->data[i][q] * U->data[q][k];temp = A->data[i][k] - temp;Insert(L, temp / U->data[k][k], i, k);temp = 0;}}Insert(U, U->data[Fnumber][Xnumbers] / U->data[Fnumber][Xnumbers - 1], Fnumber, Xnumbers + 1);for (k = Fnumber - 1; k >= 1; k--){for (q = k + 1; q < Xnumbers; q++)temp = temp + U->data[k][q] * U->data[q][Xnumbers + 1];temp = U->data[k][Xnumbers] - temp;Insert(U, temp / U->data[k][k], k, Xnumbers + 1);temp = 0;}printresx(U);}//////////////////////////迭代法////////////////////////////void ISODATAmethod(datacoa *A){int i=1, j, k=0;float x0 = 0, x1 = 1, temp = 0;for (i = 1; i <= Fnumber; i++)Insert(A, 0, i, Xnumbers + 1);while (1){for (i = 1; i <= Fnumber; i++){for (j = 1; j <= Fnumber; j++){if (j == i)continue;temp = temp + A->data[i][j] * A->data[j][Xnumbers + 1];}temp = A->data[i][Xnumbers] - temp;temp = temp / A->data[i][i];Insert(A, A->data[i][Xnumbers + 1], i, Xnumbers + 2);Insert(A, temp, i, Xnumbers + 1);temp = 0;}k++;for (i = 1; i <= Fnumber; i++)temp = temp + fabs(A->data[i][Xnumbers + 1] - A->data[i][Xnumbers + 2]);if ((temp < 0.000001) || k > max)break;temp = 0;}DeleteLie(A, Xnumbers + 2);printfunction(A);cout << endl;printresx(A);cout << k << endl;}3.4数据结构部分数据结构头文件#ifndef _SQUENCELIST_H#define _SQUENCELIST_H#include "stdio.h"#include "stdlib.h"#include "stdafx.h"#include<iostream>using namespace std;#define maxsize 1024/***sequenlist*/typedef float datatype;typedef struct{datatype data[maxsize][maxsize];int m, n;}datacoa;datacoa *InitStruct();int Length(datacoa*);int Insert(datacoa*, datatype, int ,int);void DeleteLie(datacoa*L, int j);void DeleteLine(datacoa*L, int i);#endif数据结构CPP文件#include "stdafx.h"#include "squencelist.h"///////////////////////////////////数据结构部分//////////////////////////////////////////////////////////////////////////////sequenlist/////////////////////////////////////////// datacoa *InitStruct(){datacoa*L = (datacoa*)malloc(sizeof(datacoa));L->m = 0;L->n = 0;return L;// datacoa*L = new datacoa;}int Length(datacoa*L){return L->m*L->n;}int Insert(datacoa*L, datatype x, int i,int j){int k;if ((L->m >= maxsize - 1) || (L->n >= maxsize - 1)) {cout << "表已满" << endl;return 0;}for (k = L->n; k >= j; k--)L->data[i][k + 1] = L->data[i][k];L->data[i][j] = x;if (i > L->m)L->m++;if (j > L->n)L->n++;return 1;}void DeleteLie(datacoa*L,int j){int k,i;if ((j<1) || (j>L->n)){cout << "非法删除位置" << endl;}for (i = 1; i <= L->m; i++){for (k = j; k <= L->n; k++)L->data[i][j] = L->data[i][j + 1];}L->n--;}void DeleteLine(datacoa*L, int i){int k,j;if ((i<1) || (i>L->m)){cout << "非法删除位置" << endl;}for (j = 1; j <= L->n; j++){for (k = i; k <= L->m; k++)L->data[i][j] = L->data[i+1][j];}L->m--;}4运行结果4.1列主元高斯消去法运行结果4.2LU三角分解法运行结果4.3雅克比迭代法运行结果边界情况调试超定程等可能无解的情况迭代法迭代次数超出100次的情况5总结这次的程序设计题看似简单实则临界代码区很多，调试时很容易出错。