|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
求解线性方程组
$ [+ s$ d1 ]& [& Z' tsolve,linsolve$ Q8 j* K5 ]: ^2 F# @" L0 b
例:
; H& H7 z2 k N1 P6 ]A=[5 0 4 2;1 -1 2 1;4 1 2 0;1 1 1 1];
9 n9 A( E$ ]1 ?$ i+ D%矩阵的行之间用分号隔开,元素之间用逗号或空格
* K0 @ W0 x' i8 Y lB=[3;1;1;0]
5 O& T& B+ K7 z. J9 |2 X5 e: }X=zeros(4,1);%建立一个4元列向量
/ b: S% ^( S' a0 Z9 sX=linsolve(A,B)/ Q3 _7 o5 t) s
diff(fun,var,n):对表达式fun中的变量var求n阶导数。- r; b% ]' P! b/ v6 E; M
" ~7 l! I- |, N+ _( p+ f例如:F=sym('u(x,y)*v(x,y)'); %sym()用来定义一个符号表达式
- {4 w/ i6 I4 c* F: k4 w F+ Xdiff(F); %matlab区分大小写# k6 k0 s7 J7 t! J
pretty(ans) %pretty():用习惯书写方式显示变量;ans是答案表达式$ r9 D. l7 s v: Y
& C1 n; _7 Q* h! G, k5 S5 y
. d) F' h2 s' X- P非线性方程求解 2 l/ F. B( l9 W: ^4 d" Y- X
fsolve(fun,x0,options). l3 n9 d; C1 P `
其中fun为待解方程或方程组的文件名;" z, y3 B N! n4 d6 T2 j0 c
x0位求解方程的初始向量或矩阵;1 N- ^$ |2 I+ C \$ h/ K
option为设置命令参数
. J0 N2 ?& A: t+ |1 c) z; |; y P- Y' G7 o7 m建立文件fun.m:
% S e0 p& B# @+ B7 }: L' ?; qfunction y=fun(x)
: q: D# Z* N3 i! f2 A3 jy=[x(1)-0.5*sin(x(1))-0.3*cos(x(2)), ...4 g- A) x5 K' y9 C8 b J n
x(2) - 0.5*cos(x(1))+0.3*sin(x(2))];8 ~# B# R1 v9 m0 a
, {0 T# W/ M" I% o! n5 t8 [>>clear;x0=[0.1,0.1];fsolve(@fun,x0,optimset('fsolve'))& e6 J w2 U7 H/ b X0 J4 ?
注:4 Z" Y# R8 D$ e. p
...为续行符4 w+ c R1 G, \, u3 k9 g; s
m文件必须以function为文件头,调用符为@;文件名必须与定义的函数名相同;fsolve()主要求解复杂非线性方程和方程组,求解过程是一个逼近过程。3 z# O8 m; j0 T/ P5 L9 X
" B0 s7 H% z ~( K: c1 Z9 N- F% N6 K
, p% ^! }1 i2 [Matlab求解线性方程组
i6 c5 b1 S+ H2 F4 W! F4 y; yAX=B或XA=B
* v% }( |$ Y- p# _# i在MATLAB中,求解线性方程组时,主要采用前面章节介绍的除法运算符“/”和“\”。如: ! h: w: s7 }: h" c2 z
X=A\B表示求矩阵方程AX=B的解; 0 a$ K0 h$ D9 R# |5 j, T' c) x+ |% e& Y
X=B/A表示矩阵方程XA=B的解。 / B$ b h4 A9 d( g. f+ `9 p" X; F
对方程组X=A\B,要求A和B用相同的行数,X和B有相同的列数,它的行数等于矩阵A的列数,方程X=B/A同理。 9 \+ J' ~0 F4 Z* w l
* w0 f! G5 `9 F; C+ ~
如果矩阵A不是方阵,其维数是m×n,则有: $ D) u2 C, \7 B0 O4 i) u
m=n 恰定方程,求解精确解; 2 z7 R0 a8 d+ k) |, R) b2 B
m>n 超定方程,寻求最小二乘解; 5 u4 I; }! `5 ~, X0 L
m<n 不定方程,寻求基本解,其中至多有m个非零元素。 6 o% S+ }9 l2 Z0 e+ H$ W
针对不同的情况,MATLAB将采用不同的算法来求解。 0 y% Z- ]+ A3 b1 A0 R" |
, q9 h3 N, p3 n0 U P( F+ Q4 [/ Z7 d' v3 Z; Q9 h
一.恰定方程组 8 k" V- n7 r* w* |% `. j
恰定方程组由n个未知数的n个方程构成,方程有唯一的一组解,其一般形式可用矩阵,向量写成如下形式: 3 f/ n5 F r3 T; e) f" S
Ax=b 其中A是方阵,b是一个列向量;
8 d" W3 V! Y. A. d" x( z在线性代数教科书中,最常用的方程组解法有:
7 ?, q( p, a) ^3 }5 }# [2 S(1)利用cramer公式来求解法;
2 j7 E/ j- P+ z0 N x/ N# J(2)利用矩阵求逆解法,即x=A-1b;
/ i, Z* r' W4 n. u(3)利用gaussian消去法; - T4 P% e0 ~; P3 Q' v. D+ V# J% _
(4)利用lu法求解。
+ h2 N7 Y" d9 A& a- z# R一般来说,对维数不高,条件数不大的矩阵,上面四种解法所得的结果差别不大。前三种解法的真正意义是在其理论上,而不是实际的数值计算。MATLAB中,出于对算法稳定性的考虑,行列式及逆的计算大都在lu分解的基础上进行。
/ l9 ^% }4 R3 s) E2 {) B在MATLAB中,求解这类方程组的命令十分简单,直接采用表达式:x=A\b。
4 s* ~/ K B7 w; v" k' q4 l在MATLAB的指令解释器在确认变量A非奇异后,就对它进行lu分解,并最终给出解x;若矩阵A的条件数很大,MATLAB会提醒用户注意所得解的可靠性。
! d3 `; P. w" P; f如果矩阵A是奇异的,则Ax=b的解不存在,或者存在但不唯一;如果矩阵A接近奇异时,MATLAB将给出警告信息;如果发现A是奇异的,则计算结果为inf,并且给出警告信息;如果矩阵A是病态矩阵,也会给出警告信息。
n' u, Q$ C: F( N$ H. t0 j; p) F2 q) ^* @注意:在求解方程时,尽量不要用inv(A)*b命令,而应采用A\b的解法。因为后者的计算速度比前者快、精度高,尤其当矩阵A的维数比较大时。另外,除法命令的适用行较强,对于非方阵A,也能给出最小二乘解。
7 e/ n3 {) \" G4 M* l' v9 a' C( `) \
二.超定方程组 ; r' E4 G( }6 ]: J
对于方程组Ax=b,A为n×m矩阵,如果A列满秩,且n>m。则方程组没有精确解,此时称方程组为超定方程组。线性超定方程组经常遇到的问题是数据的曲线拟合。对于超定方程,在MATLAB中,利用左除命令(x=A\b)来寻求它的最小二乘解;还可以用广义逆来求,即x=pinv(A),所得的解不一定满足Ax=b,x只是最小二乘意义上的解。左除的方法是建立在奇异值分解基础之上,由此获得的解最可靠;广义逆法是建立在对原超定方程直接进行householder变换的基础上,其算法可靠性稍逊与奇异值求解,但速度较快;
* \6 R* `; Y0 T# j, M【例7】
% ]. y& \ K1 q求解超定方程组
8 A: V% e2 a2 l2 [& t% JA=[2 -1 3;3 1 -5;4 -1 1;1 3 -13] 2 [0 f( w S) h# }" m. i$ ^) l( Q
A= ; B3 t# z0 i6 k# t
2 -1 3
) [; j# j7 N9 \3 1 -5
( v6 r9 p6 E! F. u4 -1 1 6 |( {+ ]3 j' n
1 3 -13 ?6 d" P- Y3 t' E
b=[3 0 3 -6]’; & o6 c: P6 s# g5 G' ]2 y% u
rank(A)
* \5 k9 E1 U1 |2 {5 ?9 @ans=
" R+ V9 ?* { ], C3
, M j6 }3 ^4 ~$ E3 B3 l$ n* v! cx1=A\b 5 b/ I; j% ?2 G3 S- s5 {2 s5 m8 }
x1=
% G3 L+ Z3 C! {( F3 \/ m1.0000
6 p8 q% p; C3 M. C L3 [2.0000
: g3 R7 I- S! m2 ?* v1.0000 , g; r/ n/ z, [8 d9 r, U& t# _
x2=pinv(A)*b
3 a$ w Y. I5 L$ [* g' i* Lx2=
- n l: g h. R9 n* S1.0000 + n$ T/ a0 N% s
2.0000 0 e) {2 G( G0 S! b3 r; G3 _
1.0000
+ D+ W1 w' z( P6 u" DA*x1-b 8 a- m8 u9 ~7 g/ Y) G* | d: L
ans= ) h" F a: p4 M3 Z: I4 S2 j
1.0e-014 5 p! _( ?% V/ C
-0.0888 * i8 ]2 H' a0 Y. X1 f( m) v
-0.0888 ) I4 C8 P0 l$ A8 f
-0.1332
; ~' j5 t2 s; h# _) _% Q0
+ q2 V2 \4 ~, k* R: O9 c可见x1并不是方程Ax=b的精确解,用x2=pinv(A)*b所得的解与x1相同。 " t0 L# ?/ l& F$ C6 y8 L- s
- |$ f5 U! m# K& {5 o三.欠定方程组 # S" o# R# ?3 v' J9 Z6 z
欠定方程组未知量个数多于方程个数,但理论上有无穷个解。MATLAB将寻求一个基本解,其中最多只能有m个非零元素。特解由列主元qr分解求得。 * |% O$ ^4 g% w2 g, ^) Q
【例8】 1 F" n/ p( X$ X6 D% p$ {* r' S
解欠定方程组
/ z0 B) D8 v* iA=[1 -2 1 1;1 -2 1 -1;1 -2 1 5] ) k1 `8 g! H' ]& M+ l) T
A= 8 q, o" e, ?/ c; p5 e/ K/ U9 [6 t
1 -2 1 1 2 q& \4 h1 w" c$ o- H5 I
1 -2 1 -1 " i7 B7 u' b! p+ a; A( c4 x0 A2 r# y
1 -2 1 -1 8 o0 G5 x U( s* R S& L
1 -2 1 5
: P, p8 C3 d, p, V0 P# E1 P) Ub=[1 -1 5]’
. [" q3 O8 n& z4 jx1=A\b
! B V1 T$ Q G* o* G3 lWarning:Rank deficient,rank=2 tol=4.6151e-015
1 \' m, [' Y& m6 V3 h) Lx1= 3 e, b6 f+ d3 u: w$ o5 ~& s, ?' O
0 $ y, v4 f. r" w; W8 X0 H+ a A0 L& ~
-0.0000 8 c4 ^) ?! b! D, y$ j: M( ?
0
9 K# e8 l$ d/ k5 T+ ]1.0000
- n9 v0 o' m' h9 e1 B& ?2 b! Jx2=pinv(A)*b
6 S- E7 X& F0 \& n9 Y& R1 [x2=
, F! @# g7 m# }0
1 k# y$ R1 e" E: M* I, x h2 B Z-0.0000
% d) Z8 w) A a7 T; g+ v0.0000
8 l- x3 T; F9 [# W# A) ]/ b' N1.0000 G9 v* R/ Z$ g' F; {
2 U7 l& u; O3 x2 @+ ^
四.方程组的非负最小二乘解 ! ?7 N, M' B- J, r
在某些条件下,所求的线性方程组的解出现负数是没有意义的。虽然方程组可以得到精确解,但却不能取负值解。在这种情况下,其非负最小二乘解比方程的精确解更有意义。在MATLAB中,求非负最小二乘解常用函数nnls,其调用格式为: 0 u* E, S' |# F! L; i: H
(1)X=nnls(A,b)返回方程Ax=b的最小二乘解,方程的求解过程被限制在x 的条件下; 6 t+ u# i( h# m7 U2 D$ G5 x: h
(2)X=nnls(A,b,TOL)指定误差TOL来求解,TOL的默认值为TOL=max(size(A))*norm(A,1)*eps,矩阵的-1范数越大,求解的误差越大; ( U/ @/ m( y* y# e$ j
(3)[X,W]=nnls(A,b) 当x(i)=0时,w(i)<0;当下x(i)>0时,w(i)0,同时返回一个双向量w。
0 T( Z) N9 V* W2 F+ j$ Q+ P【例9】求方程组的非负最小二乘解 4 x( b; P5 y7 K6 a# r2 H
A=[3.4336 -0.5238 0.6710
2 D5 G% e* K& F7 u9 R-0.5238 3.2833 -0.7302
" Q6 C9 A" a& Q" j, J0.6710 -0.7302 4.0261]; % E# u" {1 p; C2 H- ]) l- _
b=[-1.000 1.5000 2.5000];
6 l# c* @& R) U) T/ N5 m" X* N6 `[X,W]=nnls(A,b) 7 X- O/ ~% F1 h4 J" w9 b
X= 4 c O8 c; k& z3 m/ N
0
; h [, \1 H! `+ _+ N0 ?' E0.6563 8 T6 N- ?, u, [# l; D8 y) Y
0.6998 / Q2 _$ R; P5 j1 t
W= % b# O; O, d: t: B& C3 [+ x
-3.6820 ! z% B2 J7 n! q7 x0 X
-0.0000 # Y. O; Z. w; L" _) ?
-0.0000 $ _+ o+ g2 ~/ J
x1=A\b - l- Q/ E( n+ c: f3 z: Z: v
x1=
1 M" }$ T& s1 X9 k5 {8 I-0.3569 & l- T/ {1 Y9 T k1 o
0.5744
% `; m( V4 u, U8 b! S Q% ~0.7846
$ \* A/ g1 @: R }$ \9 {5 NA*X-b - a1 J, [4 z: G
ans=
0 Z6 f* u0 e! I4 A8 v5 P1.1258
, U# _5 u% Q4 E0 W% }# y0.1437 9 a+ B' a4 {9 W3 b7 g; F, T
-0.1616 ; j4 ~4 E: W. j* }4 p- ~
A*x1-b ( t- W. J; W4 ^: Z$ \% W
ans= . ]/ }) `2 W: v! E/ n; O6 e/ I
1.0e-0.15
i( k5 Q9 R( R$ c# Y4 t( a-0.2220
2 a. M0 Z ^1 B/ b0.4441
5 z1 z* ~2 M4 V& a |0 |
|