|
|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
求解线性方程组
4 E3 @8 ~$ w3 U4 dsolve,linsolve! Q. g' p4 n3 T9 B
例:
. ?4 \- u. ~4 l4 x3 OA=[5 0 4 2;1 -1 2 1;4 1 2 0;1 1 1 1];
2 b$ e, ?. u3 g* i' }%矩阵的行之间用分号隔开,元素之间用逗号或空格
$ p- b( B0 \5 s- k+ u( C. LB=[3;1;1;0]0 @, |) u6 k! o0 S, U6 G* u
X=zeros(4,1);%建立一个4元列向量% B4 `, c. V S7 x* Z
X=linsolve(A,B)
+ j" ]7 |# n7 a! A8 S8 Q- @" sdiff(fun,var,n):对表达式fun中的变量var求n阶导数。; _; N! g; j0 R. z9 U- u
e. R" Y( q' v1 K' s( v2 A/ T例如:F=sym('u(x,y)*v(x,y)'); %sym()用来定义一个符号表达式/ x# N2 y+ Y1 N/ d* n( X
diff(F); %matlab区分大小写8 |! u) O7 ~* J/ G7 H
pretty(ans) %pretty():用习惯书写方式显示变量;ans是答案表达式3 R0 T4 B0 b' ~( @5 x
% Z6 V& d: H6 M" R9 |& |) D9 L6 M
- c# ` z- _. d7 N+ C, o3 d, d2 B非线性方程求解 0 v2 F4 w& v n/ S6 s7 q2 v- a w
fsolve(fun,x0,options)7 K m9 m# X1 h9 {. u; B) Y0 X
其中fun为待解方程或方程组的文件名;
) z2 \; `! _! j* L. [x0位求解方程的初始向量或矩阵;
- D8 N' w9 |5 |( q+ L' f# e* Uoption为设置命令参数1 y" x4 n( S" C5 H, Y+ p
建立文件fun.m:& c7 S. G: T/ v4 e
function y=fun(x)5 N+ H, Y" U" l4 S* D
y=[x(1)-0.5*sin(x(1))-0.3*cos(x(2)), ...
( r6 T: f0 h7 f' Ox(2) - 0.5*cos(x(1))+0.3*sin(x(2))];# f7 D! {5 g* C5 e
- }, R4 Z8 E" b) V$ l>>clear;x0=[0.1,0.1];fsolve(@fun,x0,optimset('fsolve'))
6 }8 y! i3 i: c0 j: M注:
, l, c) x# W" Z0 V$ B...为续行符
# R, O) M6 V0 f6 h+ am文件必须以function为文件头,调用符为@;文件名必须与定义的函数名相同;fsolve()主要求解复杂非线性方程和方程组,求解过程是一个逼近过程。
% J; g5 A# H4 q1 Q+ h, o% v; |2 ]
/ X3 {/ k( d. a4 k/ n3 t3 c% E8 s
. A4 c; b8 t" q1 O0 b2 m1 a1 `7 @Matlab求解线性方程组2 m. m3 U9 F% D5 B. l4 F
AX=B或XA=B
; u3 a8 }$ n1 }) X% F4 m5 D在MATLAB中,求解线性方程组时,主要采用前面章节介绍的除法运算符“/”和“\”。如:
- ~ T& @* \+ J$ B4 e) pX=A\B表示求矩阵方程AX=B的解; 2 i. \$ ?4 J; C
X=B/A表示矩阵方程XA=B的解。
3 E7 K F# P3 |3 j2 h6 L对方程组X=A\B,要求A和B用相同的行数,X和B有相同的列数,它的行数等于矩阵A的列数,方程X=B/A同理。
& Z$ X. c$ }4 ~1 X' U1 v
) c6 |/ K) Z+ C1 @# Y6 H如果矩阵A不是方阵,其维数是m×n,则有: 2 N# x* c' S/ G+ i* U
m=n 恰定方程,求解精确解;
" M7 g6 D/ c, W* Z: Tm>n 超定方程,寻求最小二乘解;
8 U- z6 s( f7 l! gm<n 不定方程,寻求基本解,其中至多有m个非零元素。
* q( z+ l% h, Z- Y针对不同的情况,MATLAB将采用不同的算法来求解。 $ T- |$ u1 @2 C) g; Y l
/ Z4 v9 \$ j/ f1 n5 ]# ]# O( U+ S q3 y) o$ L' _
一.恰定方程组 6 A% q2 g/ W' m( {6 j
恰定方程组由n个未知数的n个方程构成,方程有唯一的一组解,其一般形式可用矩阵,向量写成如下形式:
$ ?9 ^ Y% P' v! n0 }7 j: zAx=b 其中A是方阵,b是一个列向量; % ]& C1 k: J- u* |
在线性代数教科书中,最常用的方程组解法有:
' R- f3 o# e$ y(1)利用cramer公式来求解法; ; x o1 J% R& S4 a/ {8 _
(2)利用矩阵求逆解法,即x=A-1b;
1 }% j- Z0 y3 u! i, A5 z" e(3)利用gaussian消去法;
3 E# S N# i/ P2 P0 h \" f(4)利用lu法求解。 ; c% a7 @9 C, h, D" ]0 Y4 m& N
一般来说,对维数不高,条件数不大的矩阵,上面四种解法所得的结果差别不大。前三种解法的真正意义是在其理论上,而不是实际的数值计算。MATLAB中,出于对算法稳定性的考虑,行列式及逆的计算大都在lu分解的基础上进行。
0 B! h' W- `) W6 i1 H) b在MATLAB中,求解这类方程组的命令十分简单,直接采用表达式:x=A\b。 8 w- j' V& a7 ^% ~' q; |
在MATLAB的指令解释器在确认变量A非奇异后,就对它进行lu分解,并最终给出解x;若矩阵A的条件数很大,MATLAB会提醒用户注意所得解的可靠性。 ; g1 @, e# E0 J; U9 k2 V
如果矩阵A是奇异的,则Ax=b的解不存在,或者存在但不唯一;如果矩阵A接近奇异时,MATLAB将给出警告信息;如果发现A是奇异的,则计算结果为inf,并且给出警告信息;如果矩阵A是病态矩阵,也会给出警告信息。
9 b' F2 y) T4 b注意:在求解方程时,尽量不要用inv(A)*b命令,而应采用A\b的解法。因为后者的计算速度比前者快、精度高,尤其当矩阵A的维数比较大时。另外,除法命令的适用行较强,对于非方阵A,也能给出最小二乘解。 ( u, T' ]4 o: o @/ p" y
w+ I; d* W$ q: o2 z
二.超定方程组 % Z' T8 p! Y3 f+ Z; g3 T/ |
对于方程组Ax=b,A为n×m矩阵,如果A列满秩,且n>m。则方程组没有精确解,此时称方程组为超定方程组。线性超定方程组经常遇到的问题是数据的曲线拟合。对于超定方程,在MATLAB中,利用左除命令(x=A\b)来寻求它的最小二乘解;还可以用广义逆来求,即x=pinv(A),所得的解不一定满足Ax=b,x只是最小二乘意义上的解。左除的方法是建立在奇异值分解基础之上,由此获得的解最可靠;广义逆法是建立在对原超定方程直接进行householder变换的基础上,其算法可靠性稍逊与奇异值求解,但速度较快;
2 E+ A, a+ k5 w5 i【例7】 ; |! H6 a( L, O: g; w
求解超定方程组 + ?* i/ F( ~' l. t7 q
A=[2 -1 3;3 1 -5;4 -1 1;1 3 -13]
: `1 i" Z3 C8 q( [+ {A= # i& t2 H& S/ j( y' ]) K% a
2 -1 3
3 m3 x& `8 r5 d5 ~3 1 -5
- D" O u/ u& n u! ?4 -1 1 + O5 Q8 d& z' p$ S( B" Y: C
1 3 -13 ' w) z' D* ]/ @2 c# O4 J, }
b=[3 0 3 -6]’; ; J- m! q# d, S
rank(A) 6 W% O# ]2 S4 U) h5 n/ h
ans=
. y) L, W" y( e; c8 d3 # t M1 E1 {- z0 V P* g+ X; X
x1=A\b + a: X8 Q( V# B
x1= 0 q2 H5 ]' s7 b& ~) S. q
1.0000
/ w. I8 T( j0 Y) q+ i2.0000
6 D" l9 d2 y( b5 B( X Y8 t1.0000 2 @4 r. S7 ?0 k
x2=pinv(A)*b
1 a5 I) b5 I2 V7 y" f! u2 cx2=
0 S6 B. U+ t$ F4 O. V- L8 W1.0000 % Z, \# e2 K, v; x: J
2.0000 & |9 S( X; _. ?7 k" P3 n
1.0000 8 |: `* B6 \( R5 e+ z. H
A*x1-b
+ K* r# o7 H! K$ p; D& }) {5 kans= ) a$ V1 X y" F ?* u
1.0e-014 0 F9 w' `. g' z" I7 ~+ U
-0.0888 5 O1 v! J2 T& s3 h1 Z
-0.0888 2 R( h7 w" H% a: A, M8 K, l( g6 X" d9 ^# _2 @
-0.1332 : S* W5 i: O" v6 U. m
0
* L$ f* ^, O- _; t7 M$ t! B可见x1并不是方程Ax=b的精确解,用x2=pinv(A)*b所得的解与x1相同。
! |5 H% R" v7 S$ w H. T! T3 F2 G" M( G! C' s% N
三.欠定方程组
+ j3 U1 {# j7 |6 r欠定方程组未知量个数多于方程个数,但理论上有无穷个解。MATLAB将寻求一个基本解,其中最多只能有m个非零元素。特解由列主元qr分解求得。
% T2 T# S% \- ^, ~【例8】
" P. L; \* o( {1 B8 t解欠定方程组
. s/ R" b5 f" N9 ^A=[1 -2 1 1;1 -2 1 -1;1 -2 1 5]
8 h8 M2 ]# M- G: o6 R8 U7 {A=
3 a1 z* I& f8 q- A" W+ r1 -2 1 1
+ ~. t" \% k4 m9 A3 `' f0 W1 -2 1 -1 $ G+ b. d9 S% Y
1 -2 1 -1
" Z+ p* E# B) e r1 K/ i1 -2 1 5
% o' l/ r$ T+ ~+ N; bb=[1 -1 5]’
* B9 m5 F8 j# `2 h8 ?x1=A\b 0 _. D7 g( q9 q4 G6 ]/ ? a5 `
Warning:Rank deficient,rank=2 tol=4.6151e-015 1 F0 E- o- D. @9 E6 f$ U8 f
x1= + r5 n' T, z( C
0
* y. q+ Q# Q2 S6 ^; S) G-0.0000 1 Z3 D/ L4 w$ G& L/ @
0
, @& M; X2 C+ O$ {& A: j1.0000
& k p6 I5 B' c* f) b4 C7 Y1 Ox2=pinv(A)*b
( N _* l& n, z) c# \) Bx2=
/ S* W' t4 Z) g' O+ k0
0 ~% _! _% Q1 a' i" _-0.0000
3 y$ U. d8 w3 ~, W( h* d0.0000 5 p$ E, z$ G+ K4 b
1.0000
" `& F, {' e' T+ g- B2 N: s
. M: E( e8 J+ Y7 i四.方程组的非负最小二乘解 $ X5 w/ s1 V( t) S' M) m( ^
在某些条件下,所求的线性方程组的解出现负数是没有意义的。虽然方程组可以得到精确解,但却不能取负值解。在这种情况下,其非负最小二乘解比方程的精确解更有意义。在MATLAB中,求非负最小二乘解常用函数nnls,其调用格式为:
! Y, u9 s: t* [7 w) X7 F$ q(1)X=nnls(A,b)返回方程Ax=b的最小二乘解,方程的求解过程被限制在x 的条件下;
' o2 \' [0 S/ N2 L& W" D(2)X=nnls(A,b,TOL)指定误差TOL来求解,TOL的默认值为TOL=max(size(A))*norm(A,1)*eps,矩阵的-1范数越大,求解的误差越大;
! J0 m* _9 R) w3 R2 @(3)[X,W]=nnls(A,b) 当x(i)=0时,w(i)<0;当下x(i)>0时,w(i)0,同时返回一个双向量w。
3 ]+ b* ?/ t9 ]3 g【例9】求方程组的非负最小二乘解 ' n, W8 v, O' _- I) i9 G
A=[3.4336 -0.5238 0.6710
" w# d' [! ~0 T' o9 S( u-0.5238 3.2833 -0.7302 3 {; ^; m" ?1 d* w" Q m& e* N
0.6710 -0.7302 4.0261]; 3 y% C0 u; ~: c$ o! [
b=[-1.000 1.5000 2.5000];
' i% }9 y& r# ]: @2 Q[X,W]=nnls(A,b)
5 U- |+ k6 ?, i7 `, A* AX= 9 F; m: b, T' o1 q
0 7 k7 G; d4 O0 I4 k" j
0.6563 2 t! S1 l' P/ T, i8 `; Z
0.6998 - D" o" l2 u f4 h9 ~% q v
W=
: u4 M* p' M' e( O7 h-3.6820
! E+ M* m1 E# N, m9 F" o4 w-0.0000
: H/ U2 N0 ]% ]6 O" ]3 g-0.0000 2 g4 g# O+ t/ i4 m& _1 e) z
x1=A\b
+ R6 U/ `/ a3 K. Y0 cx1= / }6 H- w8 {+ s
-0.3569
/ \ [) w5 D2 T {/ m7 e- {0.5744
) _; h/ z1 U" S/ e7 @* G7 \& c& ]0.7846 + K' u) p4 {3 s, I* [: \
A*X-b
4 S3 {4 j# D: A0 X$ y7 O; yans= ) H# n; a3 A- \+ u- D0 x2 r8 ~9 U
1.1258
$ X4 n' q/ I/ j$ ^1 {) V1 H/ [2 A0.1437
& Q- Z* E d: b: U5 L8 ]# w-0.1616
- O- n* ?/ Z, o4 S8 GA*x1-b $ G1 z/ m* R; V/ G6 ?7 T# P
ans= 0 G* a! J+ X, p( U! u& M$ `& y
1.0e-0.15 " k; k4 E2 p4 @4 I# w! b8 R6 Z
-0.2220
% d! L/ g' b$ q& O0.4441 ) d9 R# N+ d- X! H9 G
0 |
|