|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
求解线性方程组
, E$ v Q, R4 asolve,linsolve
3 I/ Q. E1 m0 u6 @$ N例:* l9 m- m+ n' J$ O
A=[5 0 4 2;1 -1 2 1;4 1 2 0;1 1 1 1];; W, }/ ~5 y2 N2 D+ j6 ?
%矩阵的行之间用分号隔开,元素之间用逗号或空格! n+ u: }: z9 c: _) i# ^
B=[3;1;1;0]1 h: T) V4 K4 O8 ~" G; K% Q" {& r' F
X=zeros(4,1);%建立一个4元列向量
+ V5 M j9 ?5 |( L! H/ f5 YX=linsolve(A,B): G# ^8 X0 |3 f! H
diff(fun,var,n):对表达式fun中的变量var求n阶导数。) x" f: A7 R J ]) ]8 K9 O
( s8 I) I* l; Y: g5 d2 e R {5 F
例如:F=sym('u(x,y)*v(x,y)'); %sym()用来定义一个符号表达式# o2 b3 I0 k) \# w( O2 C p
diff(F); %matlab区分大小写
5 x. u* O0 ^ r# r& g9 Xpretty(ans) %pretty():用习惯书写方式显示变量;ans是答案表达式, h ~! A. D0 [: O0 n& f* G# D
: {+ z: F/ D3 o- c9 B) N$ _/ s- |/ F, W
非线性方程求解 9 B. O( d& {6 `$ [* }$ t* |/ m
fsolve(fun,x0,options)
& S% c) `+ J( A, `其中fun为待解方程或方程组的文件名;
6 b4 n" [' t& T$ k: s: |( _x0位求解方程的初始向量或矩阵;) d5 n4 s7 }4 b7 ^" U) C
option为设置命令参数/ c+ c! D. A3 e" r2 q
建立文件fun.m: m/ ~ f' M7 Y( ~& x2 D& }
function y=fun(x)3 [4 w' O* E* V3 z" R0 A
y=[x(1)-0.5*sin(x(1))-0.3*cos(x(2)), .../ h* D' {! h6 E1 {" P1 r' i
x(2) - 0.5*cos(x(1))+0.3*sin(x(2))];7 a1 b+ y" s0 a& U2 ] d7 B L+ D
/ f9 H1 W* C' W4 A/ H8 A" V>>clear;x0=[0.1,0.1];fsolve(@fun,x0,optimset('fsolve')): h; m- }' Z. n4 q- A( r
注:- N6 }3 `3 Y4 b( _9 S
...为续行符8 u$ Z) I# f- d) d/ M& M
m文件必须以function为文件头,调用符为@;文件名必须与定义的函数名相同;fsolve()主要求解复杂非线性方程和方程组,求解过程是一个逼近过程。
G4 q* Z) U8 S0 C: V
$ `' h# a; o% q8 b3 N! F% F2 W& A
0 S# d9 `. a# }, lMatlab求解线性方程组
+ @2 l, n3 \; C4 y1 S# o: F: YAX=B或XA=B % i- ^; q/ c: i; X* `1 A" D
在MATLAB中,求解线性方程组时,主要采用前面章节介绍的除法运算符“/”和“\”。如:
: R; A% `3 e- U9 b3 ]' LX=A\B表示求矩阵方程AX=B的解;
% m+ H1 J6 p- @! \; s2 f) nX=B/A表示矩阵方程XA=B的解。 7 y; \3 E' A& Q; u& ?
对方程组X=A\B,要求A和B用相同的行数,X和B有相同的列数,它的行数等于矩阵A的列数,方程X=B/A同理。
: M% I/ Y6 J* V6 V o, s, E: A, N; P% h2 p
如果矩阵A不是方阵,其维数是m×n,则有: 9 s9 |) [1 c D$ j
m=n 恰定方程,求解精确解;
4 j7 P' P" K8 p* a& `. H' rm>n 超定方程,寻求最小二乘解; , r2 ?/ T" C) q, m2 B" i3 Z! r7 c
m<n 不定方程,寻求基本解,其中至多有m个非零元素。 \: |0 W3 a5 p/ h" n
针对不同的情况,MATLAB将采用不同的算法来求解。
( S# v3 P3 ]# f* z3 K% \3 O: m! k7 |8 n" _0 q8 d2 \
5 l4 a+ ?7 d' o7 W( f2 }" h/ H+ J8 z! |一.恰定方程组 V8 e7 `) j% i- y5 l5 j
恰定方程组由n个未知数的n个方程构成,方程有唯一的一组解,其一般形式可用矩阵,向量写成如下形式: * _* O" J* y( y0 y
Ax=b 其中A是方阵,b是一个列向量;
+ |$ v5 B+ ~ B' x5 O在线性代数教科书中,最常用的方程组解法有: 5 |5 z+ ?7 f+ o/ X5 C6 j" Q0 \5 [
(1)利用cramer公式来求解法;
4 d X7 _0 }1 W/ e; P, j(2)利用矩阵求逆解法,即x=A-1b; 7 l7 C8 ^" _8 s( D# Z
(3)利用gaussian消去法; 5 M, Y1 y D1 V8 o; m/ E
(4)利用lu法求解。 6 p4 e& Z1 g) s! [% h1 R
一般来说,对维数不高,条件数不大的矩阵,上面四种解法所得的结果差别不大。前三种解法的真正意义是在其理论上,而不是实际的数值计算。MATLAB中,出于对算法稳定性的考虑,行列式及逆的计算大都在lu分解的基础上进行。 ' g5 N; L0 t b# W
在MATLAB中,求解这类方程组的命令十分简单,直接采用表达式:x=A\b。 * u+ f9 X0 K2 f* f2 [2 \
在MATLAB的指令解释器在确认变量A非奇异后,就对它进行lu分解,并最终给出解x;若矩阵A的条件数很大,MATLAB会提醒用户注意所得解的可靠性。
: K: d6 B6 y! ~% B如果矩阵A是奇异的,则Ax=b的解不存在,或者存在但不唯一;如果矩阵A接近奇异时,MATLAB将给出警告信息;如果发现A是奇异的,则计算结果为inf,并且给出警告信息;如果矩阵A是病态矩阵,也会给出警告信息。 ' {/ y5 u. j3 \- E
注意:在求解方程时,尽量不要用inv(A)*b命令,而应采用A\b的解法。因为后者的计算速度比前者快、精度高,尤其当矩阵A的维数比较大时。另外,除法命令的适用行较强,对于非方阵A,也能给出最小二乘解。 7 d' s3 H. N: x( l$ R& }* @
. d2 u% W- k' ^/ t, K8 T
二.超定方程组 6 n$ d. D! e1 o/ ?' R+ f: `
对于方程组Ax=b,A为n×m矩阵,如果A列满秩,且n>m。则方程组没有精确解,此时称方程组为超定方程组。线性超定方程组经常遇到的问题是数据的曲线拟合。对于超定方程,在MATLAB中,利用左除命令(x=A\b)来寻求它的最小二乘解;还可以用广义逆来求,即x=pinv(A),所得的解不一定满足Ax=b,x只是最小二乘意义上的解。左除的方法是建立在奇异值分解基础之上,由此获得的解最可靠;广义逆法是建立在对原超定方程直接进行householder变换的基础上,其算法可靠性稍逊与奇异值求解,但速度较快; 5 P( S( s7 T' _
【例7】 # k+ L, h. V! C6 K! y! p4 Y
求解超定方程组
! W8 V* s9 ~3 Q, K- j' WA=[2 -1 3;3 1 -5;4 -1 1;1 3 -13] " S% {2 b5 \/ f, ~
A=
3 F6 g4 c$ W1 K" k5 [' X5 I2 -1 3
) @% ~* l6 |# i: K( K" j3 1 -5 4 D4 S- d1 \5 b' D& q% L& f4 l
4 -1 1
; w) C9 a3 G' W1 k6 j8 Y2 S; ?1 3 -13 5 O6 q- ]# F/ c' C: y
b=[3 0 3 -6]’; : Z/ R# f! ^, o$ p w
rank(A)
$ u& p. O0 V1 H1 |- D" wans= 7 @2 G7 J& t# W' b6 \- U
3
- c9 ], D( ^0 ^4 w: rx1=A\b . J; a9 h. n: \1 K/ a3 W
x1= : j6 l t. E+ a
1.0000 " o2 q, \: a' z7 L" Z" c7 @
2.0000 ; U5 A7 n) }2 _4 S! m* Q; a ?
1.0000 & K; P, V2 x! F8 w# |) C+ }4 {" k7 t
x2=pinv(A)*b# v: r: L: _+ Y5 g
x2=
: I' P: q8 F! b1.0000
. F7 W: ]# @( w9 F2 i/ k( u) e2.0000
1 y+ Z2 j: {5 d# z* m# @1.0000 ; l1 u3 H. X$ u \# q/ I$ T
A*x1-b
( N# [' ?4 e0 z; z0 Aans=
+ \) S7 x4 W H3 [& H; f. K7 p% j1.0e-014 ! ]# [; k8 v+ c5 e% _" k6 M$ \% W
-0.0888
. r* I5 y1 ^( F-0.0888
( R9 W/ X/ D( K-0.1332
9 v7 }& U! ~) ?2 {/ D0
% J9 g7 K1 }# P可见x1并不是方程Ax=b的精确解,用x2=pinv(A)*b所得的解与x1相同。
' x( l5 \' }" F5 J4 ]/ Q! S9 a3 w6 K( H1 g% ]: Z; L
三.欠定方程组 / E" x) s7 x9 Z2 [. m
欠定方程组未知量个数多于方程个数,但理论上有无穷个解。MATLAB将寻求一个基本解,其中最多只能有m个非零元素。特解由列主元qr分解求得。 3 d" @4 ]; t/ q/ R4 w
【例8】
: A. F7 m( p" Y" K9 n7 ?6 G, M4 i解欠定方程组 9 m2 E* k+ H+ _' y: x, x; `
A=[1 -2 1 1;1 -2 1 -1;1 -2 1 5]
% c. U* a* Z' b2 }; R3 wA= 5 [& [ ~- a& m+ p( x8 p. l- b
1 -2 1 1
+ B( X) N; \5 y# b; H$ Q& o& ]1 -2 1 -1 / G; d, ?( \8 e% r* T& G3 v
1 -2 1 -1
8 R" L5 n0 I; k0 x6 E! V1 B0 {1 -2 1 5 ! A$ l1 [, q+ ^& I
b=[1 -1 5]’ 8 y/ b1 y( H. l+ X
x1=A\b j2 Y: H. a9 J! `" m4 ?
Warning:Rank deficient,rank=2 tol=4.6151e-015 1 U. T) m' G( D+ ^9 `4 ?4 Y
x1= 2 R( ?# `* _% ~& |! @
0 ) Y6 K' G3 x# u* o3 E7 U
-0.0000 + r3 ]! y+ t# U2 J
0
5 R; m2 A. N4 v+ [; d/ ]1.0000
( d- J" B# I8 z3 {6 Ux2=pinv(A)*b
. T: I1 H, w3 I/ t9 _- ]x2=
+ }( D: L) j) \4 O, F% Z0
6 C+ _: ]$ ]: x! G' ?% i-0.0000 ' B7 Z1 s$ Y2 _. E7 {9 `! u4 p
0.0000 % O4 |. J0 }/ l8 K; h; Z% s! J
1.0000 - d: {* J8 O8 A0 B
! x; r- a. Y0 k7 K8 c9 a
四.方程组的非负最小二乘解 ! b& E4 u3 W" y
在某些条件下,所求的线性方程组的解出现负数是没有意义的。虽然方程组可以得到精确解,但却不能取负值解。在这种情况下,其非负最小二乘解比方程的精确解更有意义。在MATLAB中,求非负最小二乘解常用函数nnls,其调用格式为: ! s- Q) t: a. k+ t6 Q6 ]
(1)X=nnls(A,b)返回方程Ax=b的最小二乘解,方程的求解过程被限制在x 的条件下;
$ M# `* `% y: h# K+ n3 X" t# p(2)X=nnls(A,b,TOL)指定误差TOL来求解,TOL的默认值为TOL=max(size(A))*norm(A,1)*eps,矩阵的-1范数越大,求解的误差越大;
' I" q( X! }) _ q4 \' @! V5 Z(3)[X,W]=nnls(A,b) 当x(i)=0时,w(i)<0;当下x(i)>0时,w(i)0,同时返回一个双向量w。
4 Q' _1 d1 w, G/ J【例9】求方程组的非负最小二乘解
, |! N' j: E& ]7 o5 oA=[3.4336 -0.5238 0.6710
. H# ~( N1 T& K8 d& r-0.5238 3.2833 -0.7302 0 U+ Q0 f4 G3 y8 c9 N
0.6710 -0.7302 4.0261]; 1 b- j2 B9 ^9 `/ Z
b=[-1.000 1.5000 2.5000];
7 f+ R8 x& T" R) H0 ?& `# t9 ^8 P* q[X,W]=nnls(A,b) - \; A- N. O+ ]" z b! N7 U
X=
2 [( H; i$ J, `. r; [. H5 ]; y0 ! x1 x9 E, @* X+ y
0.6563 , @. }. `- Y8 V) M! M+ Q( j
0.6998 , q$ z5 r* M1 ?8 i- S& F
W= J6 ^: Y+ H: v
-3.6820
8 r9 b( w* t k3 B-0.0000
% Y/ x3 ]/ l! i, ]! A-0.0000 % T$ u: u; L( s' _8 B/ d& ~
x1=A\b
+ E/ P T. a# Ex1= ' o0 N+ Z4 C& V4 ^5 c9 O+ {
-0.3569
" v& ^" Q) Q3 `' A, U0.5744 # p1 J2 l( H. X2 P
0.7846 , n6 ?6 _1 b: J& t4 B/ V
A*X-b
& F7 a) W) m; s$ a/ }% X7 Q/ S& |" qans= 1 n4 ?$ Y$ j4 R5 G& g; `0 K7 F
1.1258
7 E6 z: U, w" d: U0.1437 ; _! h7 ~6 q" d
-0.1616 6 u3 m, x' d; k0 g8 f
A*x1-b
0 H/ l- u$ b$ T* vans= 4 ?4 z* l t+ V' _
1.0e-0.15
. a" ~* F% v# O3 D' t: ?( C-0.2220
7 X4 I4 p, {, H( i- k0 o* H; V0.4441 ) U% g: _) P7 a7 `. u9 s
0 |
|