|
|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
求解线性方程组
5 t: z9 L7 [% @3 T# T; U% P. Ysolve,linsolve* F0 {2 f; Y: |: \0 N5 Q
例:5 m4 z- T' R+ k. ? D* W
A=[5 0 4 2;1 -1 2 1;4 1 2 0;1 1 1 1];6 M% i! k m- H; W1 }
%矩阵的行之间用分号隔开,元素之间用逗号或空格1 {1 t( I8 {$ D" e
B=[3;1;1;0]4 U# g. }7 U4 z6 ]/ s
X=zeros(4,1);%建立一个4元列向量( z' z8 S9 F$ v/ V
X=linsolve(A,B) k( C/ E2 P e
diff(fun,var,n):对表达式fun中的变量var求n阶导数。( }) v9 c( r9 y
* Q, Y+ `, }0 h& k: G
例如:F=sym('u(x,y)*v(x,y)'); %sym()用来定义一个符号表达式
$ V; Z1 ^8 f+ ^( {+ {# c$ adiff(F); %matlab区分大小写
- E+ q7 i6 ^- I4 A3 cpretty(ans) %pretty():用习惯书写方式显示变量;ans是答案表达式
3 X' Q. `& a: D: x" n9 E; w; q8 O* E9 z
2 x, {; L8 A* ~/ w
非线性方程求解
+ w& ?7 |5 R7 m2 r( D+ Pfsolve(fun,x0,options)
7 j# q# B" }9 m! ?, L0 t! a( T其中fun为待解方程或方程组的文件名;
/ |8 y+ X) [9 `8 t1 r% C- Px0位求解方程的初始向量或矩阵;' ?* B; X! h& u% n+ C4 f
option为设置命令参数! F1 [! G. r0 t n( S& w+ z3 F M0 t- R# c9 `
建立文件fun.m:
# _2 }0 D. D; ~ T; `function y=fun(x)
3 g l; U* b5 p) h6 Ey=[x(1)-0.5*sin(x(1))-0.3*cos(x(2)), ...# T; {) j- a5 a
x(2) - 0.5*cos(x(1))+0.3*sin(x(2))];
5 ^, ^* j2 G. f& e
+ ]0 H7 l/ O( C. a>>clear;x0=[0.1,0.1];fsolve(@fun,x0,optimset('fsolve'))5 @6 V w5 T# D" m# |" P, f4 b
注:
+ U. e2 G4 m5 x$ J P9 |...为续行符( E J& u3 H+ `: b" \
m文件必须以function为文件头,调用符为@;文件名必须与定义的函数名相同;fsolve()主要求解复杂非线性方程和方程组,求解过程是一个逼近过程。0 x& J4 t; Z( Y; R9 d
8 U$ X7 @; s, w N1 r4 E
: U9 E3 P; a- ^! S1 @% j5 ?4 wMatlab求解线性方程组
% x: P5 r+ a; \: N/ I2 \AX=B或XA=B * P0 P0 K. e, P
在MATLAB中,求解线性方程组时,主要采用前面章节介绍的除法运算符“/”和“\”。如:
, |- [- w: t& C% o8 ~X=A\B表示求矩阵方程AX=B的解; 8 j8 b" b, i9 k1 s8 y
X=B/A表示矩阵方程XA=B的解。 ) Q. r- j* |: V' N$ S& y7 T, O
对方程组X=A\B,要求A和B用相同的行数,X和B有相同的列数,它的行数等于矩阵A的列数,方程X=B/A同理。 ! q X/ p4 c Y: e' y/ j. y) H
2 b! C, R& @% U$ D
如果矩阵A不是方阵,其维数是m×n,则有: * S& j' Z6 P5 a
m=n 恰定方程,求解精确解; " J. `$ h- L5 r$ e: ~0 |
m>n 超定方程,寻求最小二乘解; 3 t) A( P" `+ Z6 j& `; p+ W7 |
m<n 不定方程,寻求基本解,其中至多有m个非零元素。
D2 D2 Y. R' g* ^3 g9 C6 Y& P% f针对不同的情况,MATLAB将采用不同的算法来求解。
4 i& O2 [6 N2 l7 T0 I9 L1 q. w# h. K+ h
! e4 g4 n+ a( a3 e一.恰定方程组 8 O5 F. @# k0 d5 `& P& ]' C
恰定方程组由n个未知数的n个方程构成,方程有唯一的一组解,其一般形式可用矩阵,向量写成如下形式:
" N1 a# e( n6 B# F+ [& p6 n: V nAx=b 其中A是方阵,b是一个列向量;
9 z1 ?& ?1 \3 B; g4 S( q4 [7 H在线性代数教科书中,最常用的方程组解法有:
5 ]7 ? ]* j, G F2 i& f' w(1)利用cramer公式来求解法; ( a. z- L9 k6 @- k: q0 H: w
(2)利用矩阵求逆解法,即x=A-1b;
9 i" V, A4 x! l5 P(3)利用gaussian消去法; 0 Z+ U( \. R. l" j. G6 S. A. |: L
(4)利用lu法求解。
; C8 d0 Q- R. D% }" }一般来说,对维数不高,条件数不大的矩阵,上面四种解法所得的结果差别不大。前三种解法的真正意义是在其理论上,而不是实际的数值计算。MATLAB中,出于对算法稳定性的考虑,行列式及逆的计算大都在lu分解的基础上进行。
+ S& }6 P% a- E( k# d! `在MATLAB中,求解这类方程组的命令十分简单,直接采用表达式:x=A\b。 2 [8 H1 |9 h+ y' \
在MATLAB的指令解释器在确认变量A非奇异后,就对它进行lu分解,并最终给出解x;若矩阵A的条件数很大,MATLAB会提醒用户注意所得解的可靠性。
: h+ `) U# Z1 X: b如果矩阵A是奇异的,则Ax=b的解不存在,或者存在但不唯一;如果矩阵A接近奇异时,MATLAB将给出警告信息;如果发现A是奇异的,则计算结果为inf,并且给出警告信息;如果矩阵A是病态矩阵,也会给出警告信息。 " m- y7 ^, t1 k, I8 I& x+ M: x: W
注意:在求解方程时,尽量不要用inv(A)*b命令,而应采用A\b的解法。因为后者的计算速度比前者快、精度高,尤其当矩阵A的维数比较大时。另外,除法命令的适用行较强,对于非方阵A,也能给出最小二乘解。 9 A0 u) b8 \0 B0 G$ _
. q6 |0 b5 P% `7 x" L' c9 L
二.超定方程组 , r4 F# N$ P3 U; ^; u
对于方程组Ax=b,A为n×m矩阵,如果A列满秩,且n>m。则方程组没有精确解,此时称方程组为超定方程组。线性超定方程组经常遇到的问题是数据的曲线拟合。对于超定方程,在MATLAB中,利用左除命令(x=A\b)来寻求它的最小二乘解;还可以用广义逆来求,即x=pinv(A),所得的解不一定满足Ax=b,x只是最小二乘意义上的解。左除的方法是建立在奇异值分解基础之上,由此获得的解最可靠;广义逆法是建立在对原超定方程直接进行householder变换的基础上,其算法可靠性稍逊与奇异值求解,但速度较快; 6 e) `2 q5 ]2 L2 a0 X, B- C
【例7】
' t" c& q! q5 B. f& ~- {, m求解超定方程组 . h# n$ i7 y: c" F
A=[2 -1 3;3 1 -5;4 -1 1;1 3 -13]
( J+ {% o. t, B4 i# B! B. k0 W: sA=
( \8 F5 T8 C& v' u2 -1 3 # X# `# H' Z+ s# h& ?! B6 o
3 1 -5
+ l3 g Y0 T/ g) X( J4 -1 1 ' ^! ]1 ~: R- [' d! C
1 3 -13
( q1 _5 n1 K6 B( Ob=[3 0 3 -6]’;
$ ^/ t& V# n: M! y q/ urank(A)
8 i; R! c( X; m8 a2 ^' Uans=
/ ^& ?; P+ h& Q& \3
6 n F; o5 u1 p' Vx1=A\b + t3 U0 M& P$ u% b# i
x1= 7 m" d0 w* I! k( {
1.0000
1 R: T: V; w2 e4 `0 Y& f' B8 M2.0000 ! n- ~: Q4 {; g# ]' J" b
1.0000 & K5 A4 F. e- ?4 n; v9 A1 K
x2=pinv(A)*b- I z$ a. S$ J0 [" m
x2= # u- B8 n( c3 s) j5 @
1.0000 3 R. J1 Y. s# B: G1 |
2.0000 0 u9 N4 }- j3 a' _) A& S
1.0000 . i+ Y2 d; b2 N
A*x1-b
/ h4 k" Y+ T$ R! ^- c: H3 z% e* x) Ians= ) I7 s0 z* i: B8 ~& z
1.0e-014
4 N1 J8 z7 f; a) K; c. s) _-0.0888 * a( ~: W& }' K( k
-0.0888 1 |1 m& v" N7 V: _2 X k( Q4 I
-0.1332 ( h( q9 H+ G3 @3 k6 N9 [5 j
0 ; j) S$ o, ~5 X# d4 y, Q' {
可见x1并不是方程Ax=b的精确解,用x2=pinv(A)*b所得的解与x1相同。 0 R/ @! |* J' K# A" s: {# I. t2 W
+ J( ?3 I0 K; i3 p) D
三.欠定方程组 - t$ _/ q" f- x$ i% z A0 g! W2 F7 p0 o
欠定方程组未知量个数多于方程个数,但理论上有无穷个解。MATLAB将寻求一个基本解,其中最多只能有m个非零元素。特解由列主元qr分解求得。 & f0 ]* m* ?1 Z5 J- D
【例8】 5 }0 z1 C! y& X2 y0 W& c( `0 E
解欠定方程组 & b! A. T5 M/ s7 k9 B, Y" _
A=[1 -2 1 1;1 -2 1 -1;1 -2 1 5] ; N4 F9 s! i# f4 g& J( H( y
A= * N$ Q6 P4 j: z: ?! l
1 -2 1 1
6 f! }( t: L7 d& T1 -2 1 -1 " ]) D' S9 J% T
1 -2 1 -1 p8 Y* L' v1 V/ Y$ C& E' c
1 -2 1 5
a! M: C: B& E% k0 H* Ob=[1 -1 5]’ $ z6 Y1 j+ \9 {) i9 F, L. q, ^
x1=A\b
5 \! C5 S% _5 |Warning:Rank deficient,rank=2 tol=4.6151e-015
; `) U- I# ]+ N/ o# o+ sx1=
, T0 I0 _% A; M0 ^0
, s" R R; U7 I6 ^4 b7 u R& a-0.0000 1 w7 u3 W }$ ]- k2 R( j
0 7 g/ ?% `# P7 Y. t1 o- C
1.0000
% V" H' E/ b9 i. O jx2=pinv(A)*b
7 T/ q) G" F4 ^. r. u Bx2=
6 P7 F& _/ X& U: @' e" i3 \0
9 d0 P! y+ v! |/ {" v3 a-0.0000 " M$ {+ x( D) o) L
0.0000 ) C4 [ m: |7 v) g4 |( m, I
1.0000 ! V ?$ w) I* \) n
) g; f6 n7 l# @ V2 [* |1 z& y7 V8 c
四.方程组的非负最小二乘解 7 [% r# k5 @, M9 `
在某些条件下,所求的线性方程组的解出现负数是没有意义的。虽然方程组可以得到精确解,但却不能取负值解。在这种情况下,其非负最小二乘解比方程的精确解更有意义。在MATLAB中,求非负最小二乘解常用函数nnls,其调用格式为:
$ O) l. Y' K8 q9 v(1)X=nnls(A,b)返回方程Ax=b的最小二乘解,方程的求解过程被限制在x 的条件下; 0 C3 V0 w/ j7 e+ m
(2)X=nnls(A,b,TOL)指定误差TOL来求解,TOL的默认值为TOL=max(size(A))*norm(A,1)*eps,矩阵的-1范数越大,求解的误差越大; # r: u# l+ B" B5 x+ V; k7 K- y& ~
(3)[X,W]=nnls(A,b) 当x(i)=0时,w(i)<0;当下x(i)>0时,w(i)0,同时返回一个双向量w。 1 D) c \9 {! v# M+ a4 t* c
【例9】求方程组的非负最小二乘解
, ]+ ~: o4 T2 n8 ]* w" R: WA=[3.4336 -0.5238 0.6710 ) `! @. e6 k- ~! ]" o7 d! v7 M
-0.5238 3.2833 -0.7302
3 k4 U/ `1 j: F7 a- Z6 E: j9 b0.6710 -0.7302 4.0261];
* ]1 V& K" |9 Z n2 cb=[-1.000 1.5000 2.5000]; # l3 L4 u; U0 F( N" D) F8 `
[X,W]=nnls(A,b) $ ^# {: D! g d$ p" H& O: ]
X=
9 S6 Q1 I# X3 `0
. G1 g" z2 _& o3 @$ F3 Y0.6563 5 Y. R1 X$ f3 I9 y u
0.6998
: ]' w/ b$ C) a* hW=
8 H( q8 T6 z8 q s* `-3.6820 5 E. p) c, G# W7 J8 V/ m8 @. @
-0.0000 " }5 S8 c: J% g+ J U; k
-0.0000 2 \4 d; M: U" T( b9 ]+ ?
x1=A\b 9 T& c/ \2 V( a* v
x1= 9 t) e8 f: }, c. S; Q2 ?+ A
-0.3569
+ q$ ]- N% W3 m+ z0 X0.5744 $ b$ _& i) y# ?0 ^) I8 t
0.7846
( \* y$ E+ k; G* u4 c* H2 F4 lA*X-b
" U/ B; N6 H9 P* F. u9 o+ L3 bans= & ~' c$ ~4 E# Z! z( x* j& A
1.1258 K$ G: J& e" Z, U, u
0.1437
5 c0 }; o0 [3 t+ ?. r5 ^1 ~6 W-0.1616 ; W" o1 Z* K) c! k) M( i) w
A*x1-b
5 q% U! B8 d9 yans=
m/ _. r0 r/ |$ p$ W* E, i1.0e-0.15 % J: i% ^! @+ C8 t2 L w. u7 n
-0.2220
. J8 W# x3 P- F1 \0.4441 * n% j$ A8 k1 w" }1 h2 n
0 |
|