数学建模社区-数学中国
标题:
高斯消元法解线性方程组
[打印本页]
作者:
2744557306
时间:
2023-12-31 17:28
标题:
高斯消元法解线性方程组
这段代码是用于解线性方程组的高斯消元法（Gaussian elimination）。以下是代码的主要步骤：
4 y3 x, m- l- B& h- a
/ _3 Z# M- M. g& h
1.初始化： 定义系数矩阵 a 和常数向量 b，以及一个排列矩阵 L，用于记录行的交换顺序。
' h6 d. U) Y2 \$ D
6 M* R. O) C* |) e1 R
a=[1,-1,2,-1;2,-2,3,-3;1,1,1,0;1,-1,4,3]; % 系数矩阵 a
! Q" V8 U; N" [+ [( s, D
b=[-8,-20,-2,4]'; % 常数向量 b
4 l& @: t7 e# @. y; p
L=[1,2,3,4]; % 排列矩阵 L
3 U; I$ e8 f. o6 [9 C
n=length(b);
7 u }0 _% K# {, M
' W% s# { D) }; y2 r# `
7 I7 m: h' A; L2 D1 f
2.高斯消元： 通过一系列行变换将系数矩阵转化为上三角矩阵，并相应地更新常数向量。这里使用了列主元素法，即每次选取绝对值最大的元素所在的行与当前行进行交换。
( W% s9 o1 M( }/ V7 C0 U, n& M2 d% J
" M( r0 l- g& C$ n) m
for k=1:n-1
$ p) C2 v' N) J+ N( K
[p,q]=find(abs(a)==max(max(abs(a(k:n,k:n)))));
, m8 W/ N* b8 ]- D' u* B
( Y; c& S9 q) k! o
if(p~=k | q~=k)
* ?, a3 x* F$ Q7 y" U( p, O5 C
t=a(k,
;
3 ]1 ?8 m8 j6 o4 C _- ~$ @) U$ v
a(k,
=a(p,
;
+ C% W/ G( c' P; G+ F
a(p,
=t;
. s3 h; _+ ~- g/ \% M- \+ Q
r=a(:,k);
+ e7 L+ S2 Z5 G' e1 H' F/ _
a(:,k)=a(:,q);
0 u, P. u5 {7 X& E* g
a(:,q)=r;
( T4 Z6 [ O4 \9 f5 E- b! g
t=L(k);
0 M+ L+ i4 v+ @
L(k)=L(q);
( Q' b6 j0 Z1 o8 |6 K1 B
L(q)=t;
3 g9 T( ?: Y9 K! ?
u=b(k);
" d% A" u2 x; Y
b(k)=b(p);
& z8 N5 ]: k& K
b(p)=u;
9 K& R. F. c# K$ x$ ?
end
: e0 n9 ^2 s( B/ ]8 D7 f7 c
m(k+1:n,k)=a(k+1:n,k)./a(k,k);
, }6 p: Y9 g- ?/ T+ m4 o$ I0 w
a(k+1:n,k:n)=a(k+1:n,k:n)-m(k+1:n,k)*a(k,k:n);
/ B q, r6 E9 _* g
b(k+1:n)=b(k+1:n)-m(k+1:n,k)*b(k);
% m0 _9 X6 s x2 _" C7 U; E% i
end
; T, V7 A$ I0 G4 e
4 M. d6 e2 |+ e: p6 C4 Q
# u' k6 W+ S, s$ F0 U. g
3.回代： 通过回代过程求解方程组。从最后一行开始，逐步计算未知数的值。
3 D5 z! i2 S. n' g
& T& N+ D q9 [9 t, @- q. S, O- A& C
y(n)=b(n)/a(n,n);
- I% K+ H2 p$ g( V7 n8 C; }' W; C
for i=n-1:-1:1
2 S. z/ }+ t! k @
sum=0;
# A1 G2 A, [' b" X0 l
for j=i+1:n
+ g; g7 J& V, ~) B: U) z: ?
sum=sum+a(i,j)*y(j);
' I7 s# [! Q3 p' C$ t
end
J9 {+ M* c" w+ ^0 s, t
y(i)=(b(i)-sum)/a(i,i);
5 j3 i+ G5 e# N) t/ E; O4 q
end
1 t. j" D* J+ {/ A- S& ^: l, x+ |
9 i5 x) Y# d4 o' k
4 V& n2 Z; `9 a" [" J$ Q% s
4.输出结果： 将解存储在 x 中，并输出结果。
& }" d1 _: i. {' N
3 z$ F* A* k( U6 i$ Q. W' v
x(L(n))=y(n);
9 J3 k5 T" o0 O- Z4 m3 p
x(L(1:n-1))=y(1:n-1);
/ m+ h8 _! f3 g9 z) o6 v6 H
jie=x'
- [3 f$ b# F. L/ U1 W$ E5 E$ }, N0 O
) y" M0 A/ s# p+ G8 r
最后，解向量 jie 包含了线性方程组的解。请注意，这段代码在求解之前进行了列主元素的行交换，以提高数值稳定性。
; s( u! ?2 n) T6 G* p
. T/ O$ ]' \8 V: g6 W$ A
6 X5 I" U: k* p- q3 H
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5
;=a(p,;+ C% W/ G( c' P; G+ F=t;