概述
MATLAB 求解线性方程组
- 前言
- 1 线性方程组的直接解法
- 1.1利用左除运算符的直接解法
- 1.2 LU分解法
- 1.2.1 LU分解法的基本思想
- 1.2.2 Matlab的LU分解函数
- 1.2.3 用LU分解法求解线性方程组
- 2 线性方程组的迭代解法
- 2.1 基本思想
- 2.2 雅可比(Jacobi)迭代法
- 2.3 高斯-赛德尔(Gauss-Seidel)迭代法
- 总结
前言
线性方程组计算机解法有直接法和迭代法两大类。
- 直接法
用计算公式直接计算求出线性方程组解的方法。包括有高斯(Gauss)消元法、列主元消元法、全主元消元法、LU分解法(矩阵的三角分解法)、追赶法等。 - 迭代法
用迭代公式通过迭代计算求出满足精度要求的线性方程组近似解的方法。包括有Jacobi迭代法、seidel迭代法(即Gauss-seidel迭代法)、sor法(超松弛迭代法)等。(迭代法是一种逐次逼近线性方程组解的方法)
1 线性方程组的直接解法
1.1利用左除运算符的直接解法
A
x
=
b
⇒
d
e
t
(
A
)
≠
0
x
=
A
b
Ax=b ;;xRightarrow{det(A) neq 0} ;;x=Abackslash b
Ax=bdet(A)=0x=Ab
eg:
{
2
x
1
+
x
2
−
5
x
3
+
x
4
=
13
x
1
−
5
x
2
+
7
x
4
=
−
9
2
x
2
+
x
3
−
x
4
=
6
x
1
+
6
x
2
−
x
3
−
4
x
4
=
0
begin{cases} 2x_1+x_2-5x_3+x_4=13\ x_1-5x_2+7x_4 = -9 \ 2x_2+x_3-x_4=6\x_1+6x_2-x_3-4x_4=0 end{cases}
⎩⎪⎪⎪⎨⎪⎪⎪⎧2x1+x2−5x3+x4=13x1−5x2+7x4=−92x2+x3−x4=6x1+6x2−x3−4x4=0
demo:
A=[2 1 -5 1;1 -5 0 7;0 2 1 -1;1 6 -1 -4];
b=[13,-9,6,0]';
x=Ab
运行结果如下:
x =
-66.555555555555543
25.666666666666664
-18.777777777777775
26.555555555555550
1.2 LU分解法
1.2.1 LU分解法的基本思想
A x = b ⇒ A = L U L U x = b ⟹ { L y = b U x = y Ax=b quad xRightarrow{A=LU}LUx=b Longrightarrow left{begin{matrix} Ly=b \ Ux=y end{matrix}right. Ax=bA=LULUx=b⟹{Ly=bUx=y
1.2.2 Matlab的LU分解函数
- [L,U]=lu(A):产生一个上三角矩阵U和一个变换形式的下三角矩阵L,使之满足A=LU。(!注意:矩阵A必须为方阵)
- [L,U,P]=lu(A):产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P,使之满足PA=LU。(!注意:矩阵A也必须为方阵)
1.2.3 用LU分解法求解线性方程组
- A x = b ⟶ L U x = b ⟶ x = U ( L b ) Ax=blongrightarrow LUx=b longrightarrow x=Ubackslash(Lbackslash b) Ax=b⟶LUx=b⟶x=U(Lb)
- A x = b ⟶ P A x = P b ⟶ L U x = P b ⟶ x = U ( L p ∗ b ) Ax=blongrightarrow PAx=Pb longrightarrow LUx=Pblongrightarrow x=Ubackslash(Lbackslash p*b) Ax=b⟶PAx=Pb⟶LUx=Pb⟶x=U(Lp∗b)
eg:
{
2
x
1
+
x
2
−
5
x
3
+
x
4
=
13
x
1
−
5
x
2
+
7
x
4
=
−
9
2
x
2
+
x
3
−
x
4
=
6
x
1
+
6
x
2
−
x
3
−
4
x
4
=
0
begin{cases} 2x_1+x_2-5x_3+x_4=13\ x_1-5x_2+7x_4 = -9 \ 2x_2+x_3-x_4=6\x_1+6x_2-x_3-4x_4=0 end{cases}
⎩⎪⎪⎪⎨⎪⎪⎪⎧2x1+x2−5x3+x4=13x1−5x2+7x4=−92x2+x3−x4=6x1+6x2−x3−4x4=0
demo:
A=[2 1 -5 1;1 -5 0 7;0 2 1 -1;1 6 -1 -4];
b=[13,-9,6,0]';
[L,U]=lu(A);
x=U(Lb)
% [L,U,P]=lu(A);
% x=U(LP*b)
运行结果如下:
x =
-66.555555555555543
25.666666666666664
-18.777777777777775
26.555555555555550
2 线性方程组的迭代解法
2.1 基本思想
将线性方程组 Ax=b 等价变形为 x=Bx+g,然后构造向量迭代公式:
x
(
k
+
1
)
=
B
x
(
k
)
+
g
k
=
0
,
1
,
2...
x^{(k+1)}=Bx^{(k)}+g quad k=0,1,2...
x(k+1)=Bx(k)+gk=0,1,2...
再给定一个初始向量
x
(
0
)
x^{(0)}
x(0),代入迭代公式计算出迭代向量序列:
x
(
1
)
,
x
(
2
)
,
.
.
.
x
(
k
)
,
.
.
.
x^{(1)},x^{(2)},...x^{(k)},...
x(1),x(2),...x(k),...
eg:
{
10
x
1
−
x
2
=
9
−
x
1
+
10
x
2
−
2
x
3
=
7
−
2
x
2
+
10
x
3
=
6
⟹
{
x
1
=
10
x
2
−
2
x
3
−
7
x
2
=
10
x
1
−
9
x
3
=
(
6
+
2
x
2
)
/
10
begin{cases} 10x_1-x_2=9\ -x_1+10x_2-2x_3 = 7\ -2x_2+10x_3=6\ end{cases} quad Longrightarrow quadbegin{cases} x_1=10x_2-2x_3-7\ x_2=10x_1-9\ x_3=(6+2x_2)/10\ end{cases}
⎩⎪⎨⎪⎧10x1−x2=9−x1+10x2−2x3=7−2x2+10x3=6⟹⎩⎪⎨⎪⎧x1=10x2−2x3−7x2=10x1−9x3=(6+2x2)/10
⟹
{
x
1
(
k
+
1
)
=
10
x
2
(
k
)
−
2
x
3
(
k
)
−
7
x
2
(
k
+
1
)
=
10
x
1
(
k
)
−
9
x
3
(
k
+
1
)
=
(
6
+
2
x
2
(
k
)
)
/
10
Longrightarrow begin{cases} x_1^{(k+1)}=10x_2^{(k)}-2x_3^{(k)}-7\ x_2^{(k+1)}=10x_1^{(k)}-9\ x_3^{(k+1)}=(6+2x_2^{(k)})/10\ end{cases}
⟹⎩⎪⎨⎪⎧x1(k+1)=10x2(k)−2x3(k)−7x2(k+1)=10x1(k)−9x3(k+1)=(6+2x2(k))/10
2.2 雅可比(Jacobi)迭代法
A
x
=
b
⇒
(
A
=
D
−
L
−
U
)
(
D
−
L
−
U
)
x
=
b
Ax=b quad xRightarrow{(A=D-L-U)} (D-L-U)x=b
Ax=b(A=D−L−U)(D−L−U)x=b
D
=
[
a
11
a
22
⋱
a
n
n
]
L
=
−
[
0
a
21
0
⋮
⋱
⋱
a
n
1
⋯
a
n
,
n
−
1
0
]
U
=
−
[
0
a
12
⋯
a
1
n
0
⋱
⋮
⋱
a
n
−
1
,
n
0
]
D=begin{bmatrix} a_{11} & \ & a_{22}\ & & ddots & \ &&& a_{nn} end{bmatrix}quad L=-begin{bmatrix} 0 & \ a_{21}&0\ vdots & ddots & ddots & \a_{n1} &cdots &a_{n,n-1}& 0 end{bmatrix}quad U=-begin{bmatrix} 0 & a_{12}& cdots & a_{1n} \ & 0& ddots &vdots\ & &ddots &a_{n-1,n}& \ &&& 0 end{bmatrix}
D=⎣⎢⎢⎡a11a22⋱ann⎦⎥⎥⎤L=−⎣⎢⎢⎢⎡0a21⋮an10⋱⋯⋱an,n−10⎦⎥⎥⎥⎤U=−⎣⎢⎢⎢⎡0a120⋯⋱⋱a1n⋮an−1,n0⎦⎥⎥⎥⎤
求解公式为:
x
=
D
−
1
(
L
+
U
)
+
D
−
1
b
x=D^{-1}(L+U)+D^{-1}b
x=D−1(L+U)+D−1b
迭代公式为:
x
(
k
+
1
)
=
D
−
1
(
L
+
U
)
x
(
k
)
+
D
−
1
b
⇒
B
=
D
−
1
(
L
+
U
)
,
g
=
D
−
1
b
x
(
k
+
1
)
=
B
x
(
k
)
+
g
x^{(k+1)}=D^{-1}(L+U)x^{(k)}+D^{-1}b ;;xRightarrow{B=D^{-1}(L+U),g=D^{-1}b};; x^{(k+1)}=Bx^{(k)} +g
x(k+1)=D−1(L+U)x(k)+D−1bB=D−1(L+U),g=D−1bx(k+1)=Bx(k)+g
代码实现:
(雅可比迭代法的函数文件Jacobi.m)
function [y,n]=jacobi(A,b,x0,ep) %A:系数矩阵;b:常数矩阵;x0:迭代初值;ep:迭代精度;
D=diag(diag(A)); %取A的主对角线上的元素建立对角矩阵;
L=-tril(A,-1); %取A的主对角线以下(不包括主对角线)的元素建立下三角矩阵;
U=-triu(A,1); %取A的主对角线以上(不包括主对角线)的元素建立上三角矩阵;
B=D(L+U);
f=Db;
y=B*x0+g;
n=1;
while norm(y-x0)>=ep %用范数判断误差是否满足精度要求
x0=y;
y=B*x0+g;
n=n+1;
end
2.3 高斯-赛德尔(Gauss-Seidel)迭代法
D
x
(
k
+
1
)
=
(
L
+
U
)
x
(
k
)
+
b
⟹
D
x
(
k
+
1
)
=
L
x
k
+
1
+
U
x
(
k
)
+
b
⟹
(
D
−
L
)
x
(
k
+
1
)
=
U
x
(
k
)
+
b
Dx^{(k+1)}=(L+U)x^{(k)}+b;Longrightarrow ;Dx^{(k+1)}=Lx^{k+1}+Ux^{(k)}+b ;Longrightarrow ;(D-L)x^{(k+1)}=Ux^{(k)}+b
Dx(k+1)=(L+U)x(k)+b⟹Dx(k+1)=Lxk+1+Ux(k)+b⟹(D−L)x(k+1)=Ux(k)+b
x
(
k
+
1
)
=
(
D
−
L
)
−
1
U
x
(
k
)
+
(
D
−
L
)
−
1
b
→
B
=
(
D
−
L
)
−
1
U
,
g
=
(
D
−
L
)
−
1
b
x
(
k
+
1
)
=
B
x
(
k
)
+
g
x^{(k+1)}=(D-L)^{-1}Ux^{(k)}+(D-L)^{-1}b;xrightarrow{B=(D-L)^{-1}U,g=(D-L)^{-1}b};x^{(k+1)}=Bx^{(k)}+g
x(k+1)=(D−L)−1Ux(k)+(D−L)−1bB=(D−L)−1U,g=(D−L)−1bx(k+1)=Bx(k)+g
代码实现:
(Gauss-Seidel迭代法的函数文件gauseidel.m)
function [y,n]=gauseidel(A,b,x0,ep) %A:系数矩阵;b:常数矩阵;x0:迭代初值;ep:迭代精度;
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
B=(D-L)U;
f=(D-L)b;
y=B*x0+g;
n=1;
while norm(y-x0)>=ep %用范数判断误差是否满足精度要求
x0=y;
y=B*x0+g;
n=n+1;
end
eg1:
分别用雅可比迭代法和高斯-赛德尔迭代法求解线性方程组。设迭代初值为0,迭代精度为
1
0
−
6
10^{-6}
10−6。
[
4
−
2
1
−
2
4
3
−
1
−
3
3
]
[
x
1
x
2
x
3
]
=
[
1
5
0
]
begin{bmatrix} 4 quad-2quad1 \ -2quad4quad3 \ -1;-3quad3 end{bmatrix} begin{bmatrix} x_1 \ x_2\x_3 end{bmatrix}=begin{bmatrix} 1 \ 5\0 end{bmatrix}
⎣⎡4−21−243−1−33⎦⎤⎣⎡x1x2x3⎦⎤=⎣⎡150⎦⎤
demo:
clc,clear;
A=[4 -2 -1;-2 4 3;-1 -3 3];
b=[1 5 0]';
[x,n]=jacobi(A,b,[0 0 0]',1.0e-6)
[x1,n1]=gauseidel(A,b,[0 0 0]',1.0e-6)
运行结果如下:
x =
0.970588483091126
0.852941039929380
1.176470992801838
n =
35
x1 =
0.970588105856975
0.852941168882785
1.176470537501776
n1 =
16
eg2:
分别用雅可比迭代法和高斯-赛德尔迭代法求解下列线性方程组,看是否收敛。
[
1
2
−
2
1
1
1
2
2
1
]
[
x
1
x
2
x
3
]
=
[
9
7
6
]
begin{bmatrix} 1 quad2;-2 \ 1quad1quad1\ 2quad2quad1 end{bmatrix} begin{bmatrix} x_1 \ x_2\x_3 end{bmatrix}=begin{bmatrix} 9 \ 7\6 end{bmatrix}
⎣⎡12−2111221⎦⎤⎣⎡x1x2x3⎦⎤=⎣⎡976⎦⎤
demo:
clc,clear;
A1=[1 2 -2;1 1 1;2 2 1];
b1=[9;7;6];
[x11,n11]=jacobi(A1,b1,[0 0 0]',1.0e-6)
[x22,n22]=gauseidel(A1,b1,[0 0 0]',1.0e-6)
运行结果如下:
x11 =
-27
26
8
n11 =
4
x22 =
NaN
NaN
NaN
n22 =
1012
总结
- 直接法:以矩阵初等变换为基础,可以求得方程的精确解;占用的内存空间大、程序实现较复杂;一般适合求解低阶稠密线性方程组。
- 迭代法:从给定的初始值逐步逼近精确解的过程,求解过程占用存储空间小,程序设计简单;适用于求解大型稀疏矩阵线性方程组;要考虑算法的收敛性。
补充
迭代法收敛性判断:
- 迭代收敛定理: x ( k + 1 ) = B x ( k ) + g 对 任 意 x ( 0 ) 都 收 敛 ⟺ ρ ( B ) < 1 ( ρ ( B ) 为 矩 阵 B 的 谱 半 径 ) ;x^{(k+1)}=Bx^{(k)}+g;对任意x^{(0)}都收敛;iff;rho(B)<1quadquad(rho(B)为矩阵B的谱半径) x(k+1)=Bx(k)+g对任意x(0)都收敛⟺ρ(B)<1(ρ(B)为矩阵B的谱半径)
- 迭代收敛判别条件:
判别条件Ⅰ | 若存在迭代矩阵B的某种矩阵范数满足 ∣ ∣ B ∣ ∣ midmid Bmidmid ∣∣B∣∣<1,则 ∀ x ( 0 ) ∈ R n , x ( k + 1 ) = B x ( k ) + g , forall x^{(0)}in R^{n},x^{(k+1)}=Bx^{(k)}+g,; ∀x(0)∈Rn,x(k+1)=Bx(k)+g,产生的迭代序列都收敛到其不动点 x ∗ 。 ( 充 分 条 件 , 非 充 要 条 件 ! ; 矩 阵 范 数 有 一 个 小 于 1 即 可 ) x^*。(充分条件,非充要条件!;矩阵范数有一个小于1即可) x∗。(充分条件,非充要条件!;矩阵范数有一个小于1即可) |
判别条件Ⅱ | 若A为严格对角占优矩阵,则线性方程组Ax=b的Jacobi和Seidel迭代对任何初值 x ( 0 ) x^{(0)} x(0)都收敛。 |
判别条件Ⅲ | 若A为正定矩阵,则线性方程组Ax=b的Seidel迭代对任何初值 x ( 0 ) x^{(0)} x(0)都收敛。 |
如果有用,请记得给我点赞!(^ 0^)(//∇//)
最后
以上就是欢喜康乃馨为你收集整理的MATLAB直接法和迭代法求解线性方程组前言1 线性方程组的直接解法2 线性方程组的迭代解法总结的全部内容,希望文章能够帮你解决MATLAB直接法和迭代法求解线性方程组前言1 线性方程组的直接解法2 线性方程组的迭代解法总结所遇到的程序开发问题。
如果觉得靠谱客网站的内容还不错,欢迎将靠谱客网站推荐给程序员好友。
发表评论 取消回复