找回密码
 注册
关于网站域名变更的通知
查看: 1884|回复: 1
打印 上一主题 下一主题

Matlab求解线性方程组、非线性方程组

[复制链接]

该用户从未签到

跳转到指定楼层
1#
发表于 2018-10-24 13:30 | 只看该作者 回帖奖励 |正序浏览 |阅读模式

EDA365欢迎您登录!

您需要 登录 才可以下载或查看,没有帐号?注册

x
求解线性方程组
; i; X. U# c. i' I' z- Csolve,linsolve6 a& i( x7 F5 N8 d+ F" T8 K
例:. _1 `+ T. ?! {! A, `
A=[5 0 4 2;1 -1 2 1;4 1 2 0;1 1 1 1];
3 q. N+ X' z- v" Z5 u%矩阵的行之间用分号隔开,元素之间用逗号或空格
# _4 q2 S) z$ M( mB=[3;1;1;0]4 {7 m: X6 ~4 K: O" Q
X=zeros(4,1);%建立一个4元列向量
/ x# H1 |: l9 I  M9 c" C) t8 qX=linsolve(A,B)/ ~, }% I- o) a& ^
diff(fun,var,n):对表达式fun中的变量var求n阶导数。+ s/ b2 H( e0 }+ [
9 d# q, O& x7 Y, k
例如:F=sym('u(x,y)*v(x,y)'); %sym()用来定义一个符号表达式
6 c4 ]3 }0 T1 Z5 x: C5 q3 H! vdiff(F); %matlab区分大小写
0 |5 h- U6 z9 J( E; S- opretty(ans) %pretty():用习惯书写方式显示变量;ans是答案表达式
7 ~; D1 T* |9 k" _3 _; \6 `- Z0 Q' _5 d3 w- P$ [" Z3 U
- h5 [, R' S8 q, M
非线性方程求解 6 e% @2 e, s. O0 W+ P
fsolve(fun,x0,options)
$ N5 k2 \# W7 A& O3 Z. c. e* u其中fun为待解方程或方程组的文件名;2 V2 |/ W3 x4 D, A0 V. R
x0位求解方程的初始向量或矩阵;
# d+ a. {2 e( D* Moption为设置命令参数
1 r- E! W+ c! o! o9 @建立文件fun.m:
* g& J- y! O" j: `3 K- u( ifunction y=fun(x), U8 E: g: i* u) e3 d5 N& |8 k
y=[x(1)-0.5*sin(x(1))-0.3*cos(x(2)), ...
0 J$ Y# r  H1 k( Rx(2) - 0.5*cos(x(1))+0.3*sin(x(2))];
1 C/ A3 z0 H0 D0 n1 i* @3 N$ \/ G
>>clear;x0=[0.1,0.1];fsolve(@fun,x0,optimset('fsolve'))
1 \2 x6 V# [1 y注:
: ?, r) K) I5 U  C- h...为续行符
5 T* Y1 D8 y3 |/ m! Gm文件必须以function为文件头,调用符为@;文件名必须与定义的函数名相同;fsolve()主要求解复杂非线性方程和方程组,求解过程是一个逼近过程
4 `1 d5 Z: c6 f3 t( `* I  x6 H# l
+ @8 p0 D, }6 f$ `$ r* ?; @; J  ?
Matlab求解线性方程组
$ y: x# k* G8 \2 [1 RAX=B或XA=B 5 D. s" X) B# N) r( F# W
在MATLAB中,求解线性方程组时,主要采用前面章节介绍的除法运算符“/”和“\”。如: 7 t) }6 [% V- c7 Z6 p& y0 Z  D
X=A\B表示求矩阵方程AX=B的解;
) J4 A, k0 P9 s" g: |X=B/A表示矩阵方程XA=B的解。 : Y& R; D2 g& R
对方程组X=A\B,要求A和B用相同的行数,X和B有相同的列数,它的行数等于矩阵A的列数,方程X=B/A同理。 9 i% \8 N; `& e) ^# d& n" Z/ i
5 |  ]3 M7 O- r* Y, e1 w
如果矩阵A不是方阵,其维数是m×n,则有:
$ D* ^5 O+ X2 G, t+ _7 Km=n 恰定方程,求解精确解;
" I# N4 l5 Q! ~& H3 c& Q6 qm>n 超定方程,寻求最小二乘解;
' v! k# z  V9 d3 N1 s! Am<n 不定方程,寻求基本解,其中至多有m个非零元素。 ( }( F- {9 g6 e$ J9 ~
针对不同的情况,MATLAB将采用不同的算法来求解。
9 G2 p0 l6 k/ a0 a; q5 P1 J( q; P( A/ U
6 D* s( r  c, t% F
一.恰定方程组
) N+ L2 q3 g+ f& }, P恰定方程组由n个未知数的n个方程构成,方程有唯一的一组解,其一般形式可用矩阵,向量写成如下形式:
/ @9 B( h3 n: U( {Ax=b 其中A是方阵,b是一个列向量; 2 O$ l/ f9 C% t7 o% {8 `: A
在线性代数教科书中,最常用的方程组解法有:
& Q* X/ G( m+ N$ x7 `$ E; P% I(1)利用cramer公式来求解法;
1 M$ w9 n) |* Y: C2 }: @% z! f(2)利用矩阵求逆解法,即x=A-1b;
, z& j7 `# B0 `  a7 J(3)利用gaussian消去法; 9 C; U% E8 w# X5 M8 v8 ?
(4)利用lu法求解。 5 N: s! T# j. P( x: w! E. d
一般来说,对维数不高,条件数不大的矩阵,上面四种解法所得的结果差别不大。前三种解法的真正意义是在其理论上,而不是实际的数值计算。MATLAB中,出于对算法稳定性的考虑,行列式及逆的计算大都在lu分解的基础上进行。 5 ^  |1 S3 M/ \
在MATLAB中,求解这类方程组的命令十分简单,直接采用表达式:x=A\b。
5 g, n2 P! v# C4 i在MATLAB的指令解释器在确认变量A非奇异后,就对它进行lu分解,并最终给出解x;若矩阵A的条件数很大,MATLAB会提醒用户注意所得解的可靠性。 2 ]0 d, \: D0 Q, A+ \$ t# K3 [
如果矩阵A是奇异的,则Ax=b的解不存在,或者存在但不唯一;如果矩阵A接近奇异时,MATLAB将给出警告信息;如果发现A是奇异的,则计算结果为inf,并且给出警告信息;如果矩阵A是病态矩阵,也会给出警告信息。
5 W) T, \% R  H& Q( @, _注意:在求解方程时,尽量不要用inv(A)*b命令,而应采用A\b的解法。因为后者的计算速度比前者快、精度高,尤其当矩阵A的维数比较大时。另外,除法命令的适用行较强,对于非方阵A,也能给出最小二乘解。 # Y) A; h) k# J: M

8 y! J& \2 E& q( S0 m+ S二.超定方程组
) C7 _8 l8 I- x& l6 s0 t' L对于方程组Ax=b,A为n×m矩阵,如果A列满秩,且n>m。则方程组没有精确解,此时称方程组为超定方程组。线性超定方程组经常遇到的问题是数据的曲线拟合。对于超定方程,在MATLAB中,利用左除命令(x=A\b)来寻求它的最小二乘解;还可以用广义逆来求,即x=pinv(A),所得的解不一定满足Ax=b,x只是最小二乘意义上的解。左除的方法是建立在奇异值分解基础之上,由此获得的解最可靠;广义逆法是建立在对原超定方程直接进行householder变换的基础上,其算法可靠性稍逊与奇异值求解,但速度较快; 5 i0 T( p# L! R1 [7 F& ~
【例7】 . I% C: D  N4 z$ K8 p- R4 @
求解超定方程组 # \! U: W' R* s1 y6 ^
A=[2 -1 3;3 1 -5;4 -1 1;1 3 -13] $ J2 m, j% H5 z7 R  L; A8 v
A=   B0 `4 i* K- g& b
2 -1 3 ; N6 q3 ]* N. z2 V1 U# q" Y% k
3 1 -5 3 ~! d: `# i; t: y. H; o) M% G
4 -1 1
$ j! w* W$ Y2 [+ O2 o2 `6 x* P1 3 -13
" z& J% S" r$ c$ a; r& hb=[3 0 3 -6]’; 7 v4 p  @8 O6 z
rank(A)
9 v% W! Z# p( Z( Lans=
) U8 \* f9 [4 F4 l. P6 Z$ S7 u3 ' E4 j. m# {. T- \& {1 G* \0 I
x1=A\b
# X) @2 z7 E3 w2 n0 A, M4 `3 Nx1=
; d5 ]/ B2 F# Q. E: X% S7 @1.0000
! E; T3 j, F6 o& L6 ]6 d6 T2.0000
% H& M+ v, E( R/ ]1.0000 0 K7 [: q0 Q! m; Q3 S; y$ {
x2=pinv(A)*b7 i. M: X$ ?4 G8 |
x2= / U* \- Q. }$ L$ M# [) v
1.0000
/ A2 s; ~8 [1 Z  C2.0000
+ D! L" M& V" t& `- b/ G1.0000
+ P2 s' ]4 @* W* r6 c* wA*x1-b % S2 u* {8 {% B
ans= % Y: N. c/ i% y
1.0e-014
5 I+ K  ?' A& A* r5 u-0.0888 . i6 t+ K/ z5 D# o, Z( }0 j' q
-0.0888 2 A. {1 ~+ v6 ]9 D8 }: o) M- O4 d6 T
-0.1332
4 ~, l( A' G/ a* ^* H0
* l( u' F5 ~5 z可见x1并不是方程Ax=b的精确解,用x2=pinv(A)*b所得的解与x1相同。
5 @1 f6 b- e# W6 h4 U4 \. @9 \
1 `1 t' c  o9 E三.欠定方程组 1 M- D  H! c4 U" L3 \
欠定方程组未知量个数多于方程个数,但理论上有无穷个解。MATLAB将寻求一个基本解,其中最多只能有m个非零元素。特解由列主元qr分解求得。 ) L" M0 E% c, o% A
【例8】 $ {' j9 M, B, j) h; ^8 w
解欠定方程组
/ C! x) b2 `/ OA=[1 -2 1 1;1 -2 1 -1;1 -2 1 5] + O  G6 D- s6 W; H7 G0 }
A= 3 \1 G, `5 w/ T# e, ?$ o
1 -2 1 1 ( ~$ A) J9 K1 J( R
1 -2 1 -1 * b* ^8 O. _8 `) h% y
1 -2 1 -1
( m& d" E) M" j; W$ d1 -2 1 5
# B% x0 d1 F7 {+ E$ }( q6 Sb=[1 -1 5]’ % Z7 ^1 @$ ?8 A8 A
x1=A\b
! z: U! J; P" f0 p0 f  cWarning:Rank deficient,rank=2 tol=4.6151e-015 % C/ k* L5 w; y' ~# J
x1= 5 Q- {* l7 B/ j( a
0 2 B( P8 s$ J, q1 A) ^
-0.0000
2 e, Q/ S0 K+ l3 }9 Z7 Q) C0
% n, i+ u8 ?# A: _" P1.0000 , H2 n7 Q, [" d9 f* Y- \0 |
x2=pinv(A)*b # J( l1 i' k7 T  T% I2 K+ k' |
x2= 9 T% `5 w  a9 D; R/ [6 b# g% K) w1 i
0 7 i5 ]1 O/ {& o1 M4 f
-0.0000 + v3 ^3 Z- `. ^
0.0000 , q8 G! g# `! b* W) O& g, ]: C9 j
1.0000
8 ~/ F$ @1 t2 O, f$ g& R  y8 @
7 g8 u% P, [7 o6 A: m0 N四.方程组的非负最小二乘解 ; T# n8 b" x  }5 L. j- N
在某些条件下,所求的线性方程组的解出现负数是没有意义的。虽然方程组可以得到精确解,但却不能取负值解。在这种情况下,其非负最小二乘解比方程的精确解更有意义。在MATLAB中,求非负最小二乘解常用函数nnls,其调用格式为: ) G2 x7 g) s) M% e% D) p
(1)X=nnls(A,b)返回方程Ax=b的最小二乘解,方程的求解过程被限制在x 的条件下; - K4 P& O; x4 K3 O8 [7 g: E
(2)X=nnls(A,b,TOL)指定误差TOL来求解,TOL的默认值为TOL=max(size(A))*norm(A,1)*eps,矩阵的-1范数越大,求解的误差越大;
6 @6 @2 y6 n6 d% f8 l8 E1 m(3)[X,W]=nnls(A,b) 当x(i)=0时,w(i)<0;当下x(i)>0时,w(i)0,同时返回一个双向量w。 ; Y- m( f2 l7 J( I6 R' b% }% \
【例9】求方程组的非负最小二乘解
6 Q0 ?5 L6 y/ b; L" [" {A=[3.4336 -0.5238 0.6710
( U" z: r/ a  o' m) M1 K8 g-0.5238 3.2833 -0.7302
% Z" s4 \! ]6 Y: |" D1 f$ Q8 l4 ~0.6710 -0.7302 4.0261];
  l8 |% `4 B; K1 L$ fb=[-1.000 1.5000 2.5000];
! f8 o4 Q, o4 p) R( s& o! R" r% t[X,W]=nnls(A,b) - M+ t' ^' B' x) b: [% R; `- A
X= 3 E! {6 z2 b. _8 a9 z; D/ R  i4 D
0 & b8 R; j  K8 z$ I( |, Z
0.6563 3 B7 U2 v) i, K% m; p
0.6998
+ n/ y; ]# v$ C! Q! {W= + M/ e/ H' N& `5 M
-3.6820
, J+ o% u  `+ b$ d, X-0.0000 " B6 U- P# v# d; z/ x
-0.0000
" z$ u6 q8 r3 f$ X# _x1=A\b
9 _9 K/ ^4 K. d8 r& @+ h6 {x1= ( K4 L! a. ~6 e% v' m( K% @
-0.3569
7 S4 ~, L- S0 Q- K0.5744
5 Z( f! P( h% G) x0.7846
: F' u, Y1 g7 E  o; PA*X-b ( j+ o# D% n9 ~  G5 P
ans=
. M2 R3 i5 v; g4 \1 u1.1258
3 n4 w( j0 Y4 w( ?8 k5 p. ~% M0.1437 % u) \5 c2 g! H
-0.1616
  k$ u2 h1 a5 @3 w5 AA*x1-b
. k: T) B" `' E4 N. uans=
! S" Y4 K1 j( M3 g3 i! g* g1.0e-0.15
6 n7 T5 `; n/ A  t' [-0.2220
! V$ k& D7 L: w6 Y% z* `0.4441 % j" D1 }9 R+ D& N9 b
0
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

推荐内容上一条 /1 下一条

EDA365公众号

关于我们|手机版|EDA365电子论坛网 ( 粤ICP备18020198号-1 )

GMT+8, 2025-11-23 15:35 , Processed in 0.156250 second(s), 24 queries , Gzip On.

深圳市墨知创新科技有限公司

地址:深圳市南山区科技生态园2栋A座805 电话:19926409050

快速回复 返回顶部 返回列表