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

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

[复制链接]

该用户从未签到

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

EDA365欢迎您登录!

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

x
求解线性方程组: A. ]+ H3 P9 c  V# }: a
solve,linsolve- G& D! V! G) c
例:
2 c& o6 c3 F! |" N) m1 a) OA=[5 0 4 2;1 -1 2 1;4 1 2 0;1 1 1 1];' _) f- @' e& n' D  m2 g% A: u
%矩阵的行之间用分号隔开,元素之间用逗号或空格: C$ F4 n( a# v. P! F
B=[3;1;1;0]
3 z% p2 s( S* G4 t; H3 MX=zeros(4,1);%建立一个4元列向量% u9 K' ]+ \; |3 O+ E1 C
X=linsolve(A,B)
8 F# ]8 P' A: S! Y* C% t* fdiff(fun,var,n):对表达式fun中的变量var求n阶导数。
. ~* Q# [" E5 f' z) _& O8 D' R+ f. J7 E
例如:F=sym('u(x,y)*v(x,y)'); %sym()用来定义一个符号表达式
4 q6 R, |6 I0 J# Cdiff(F); %matlab区分大小写
8 G" U5 i4 j; V  Z9 apretty(ans) %pretty():用习惯书写方式显示变量;ans是答案表达式' A) I2 H, D# @4 ?% }3 g8 e
+ e, L/ A8 Q' p: V

$ H- {* z. o, R1 }& E" a1 q2 s  [非线性方程求解
5 ^$ q: t4 }/ L. ~1 i- y1 E  Afsolve(fun,x0,options)
$ y& d. r; A: S其中fun为待解方程或方程组的文件名;
: ]$ _- z) G+ s1 b9 _1 _x0位求解方程的初始向量或矩阵;
+ j" K9 Q$ g0 {0 ~/ ?option为设置命令参数
- J2 H7 o7 Y  V建立文件fun.m:
3 K! j5 x0 X! P8 qfunction y=fun(x); Z) J, ?+ L. U/ L( c
y=[x(1)-0.5*sin(x(1))-0.3*cos(x(2)), ...6 N% m; x5 N6 b; |1 L- K% f) n! N
x(2) - 0.5*cos(x(1))+0.3*sin(x(2))];
7 J' o& f0 s7 n
1 o7 Y$ [, M; e9 Q" R9 G5 B>>clear;x0=[0.1,0.1];fsolve(@fun,x0,optimset('fsolve'))* O9 Y0 {- r  h. f+ v1 S' f
注:
( X: e9 u: E( v% j& f' ?- j5 r...为续行符: ^8 j6 r$ l; n" B  i* [
m文件必须以function为文件头,调用符为@;文件名必须与定义的函数名相同;fsolve()主要求解复杂非线性方程和方程组,求解过程是一个逼近过程
: D* Q9 b5 P' |; i6 e; u6 A
, c' c8 h2 V( D% Z/ K
" E1 K  y" X& SMatlab求解线性方程组
8 M9 ]& u* i9 t; x+ I! ^; T7 ]AX=B或XA=B
6 {( h- `) }2 Y- v+ M# k在MATLAB中,求解线性方程组时,主要采用前面章节介绍的除法运算符“/”和“\”。如: 5 V# q/ V3 A0 u/ J& U
X=A\B表示求矩阵方程AX=B的解; : V- @% G( P1 f8 d
X=B/A表示矩阵方程XA=B的解。
5 z- P( P# c5 |9 z对方程组X=A\B,要求A和B用相同的行数,X和B有相同的列数,它的行数等于矩阵A的列数,方程X=B/A同理。 ' q2 I3 J/ P" j( d

$ a* D% V2 V$ V# l6 F6 @如果矩阵A不是方阵,其维数是m×n,则有: % f' L+ K" h. E. `! p2 ]. {! M
m=n 恰定方程,求解精确解;
+ [5 c! e" y0 ~; A* D9 P* _6 ym>n 超定方程,寻求最小二乘解; , {& K. `" g$ M/ K, [" }5 K
m<n 不定方程,寻求基本解,其中至多有m个非零元素。
* |: H1 `3 W; r3 i* V: e针对不同的情况,MATLAB将采用不同的算法来求解。
/ V) Z8 f" l* V6 {3 u7 _2 R! B$ M3 m( P9 X8 J5 C! n+ R

8 _/ |3 M* L5 p1 ~) j一.恰定方程组
$ q1 J, F$ z5 _% V: g# ^- h恰定方程组由n个未知数的n个方程构成,方程有唯一的一组解,其一般形式可用矩阵,向量写成如下形式: - ~2 N7 ]7 Q# M" z( J
Ax=b 其中A是方阵,b是一个列向量;
2 @8 e9 `5 b& {. y8 m在线性代数教科书中,最常用的方程组解法有: - ?1 i$ B7 O. E, G
(1)利用cramer公式来求解法;
/ m' w. n# Q) Y1 J3 B7 @! t3 Y% i1 z! u(2)利用矩阵求逆解法,即x=A-1b; . f' p  A7 l9 h1 K: J* Y8 Q2 ]8 q
(3)利用gaussian消去法;
# J( L. |6 b# R+ B% l3 X(4)利用lu法求解。 + m# h# t$ V- `
一般来说,对维数不高,条件数不大的矩阵,上面四种解法所得的结果差别不大。前三种解法的真正意义是在其理论上,而不是实际的数值计算。MATLAB中,出于对算法稳定性的考虑,行列式及逆的计算大都在lu分解的基础上进行。 ; \# \% F1 i$ e& t% V/ P5 g7 L, u
在MATLAB中,求解这类方程组的命令十分简单,直接采用表达式:x=A\b。
& a0 q2 a% x" b) m( o7 H( A2 o: u在MATLAB的指令解释器在确认变量A非奇异后,就对它进行lu分解,并最终给出解x;若矩阵A的条件数很大,MATLAB会提醒用户注意所得解的可靠性。 ) g7 f" g/ X$ a
如果矩阵A是奇异的,则Ax=b的解不存在,或者存在但不唯一;如果矩阵A接近奇异时,MATLAB将给出警告信息;如果发现A是奇异的,则计算结果为inf,并且给出警告信息;如果矩阵A是病态矩阵,也会给出警告信息。 ; B% l3 B& p0 ~! F6 _
注意:在求解方程时,尽量不要用inv(A)*b命令,而应采用A\b的解法。因为后者的计算速度比前者快、精度高,尤其当矩阵A的维数比较大时。另外,除法命令的适用行较强,对于非方阵A,也能给出最小二乘解。 $ k  R; I6 e* n. c0 P( Q4 k

6 m( Q* B$ e- K2 d+ b( i  u" ?二.超定方程组 + [  n+ {6 t3 q  @, u8 F6 i
对于方程组Ax=b,A为n×m矩阵,如果A列满秩,且n>m。则方程组没有精确解,此时称方程组为超定方程组。线性超定方程组经常遇到的问题是数据的曲线拟合。对于超定方程,在MATLAB中,利用左除命令(x=A\b)来寻求它的最小二乘解;还可以用广义逆来求,即x=pinv(A),所得的解不一定满足Ax=b,x只是最小二乘意义上的解。左除的方法是建立在奇异值分解基础之上,由此获得的解最可靠;广义逆法是建立在对原超定方程直接进行householder变换的基础上,其算法可靠性稍逊与奇异值求解,但速度较快;
0 m) e: j' u( k3 `4 m, U【例7】   y/ [! ]$ a  o8 ~% o! D
求解超定方程组 $ S0 y$ |1 z+ i! t! E
A=[2 -1 3;3 1 -5;4 -1 1;1 3 -13]
9 D; m( X) D$ j2 EA= 1 u- D! y5 u# Q3 s$ l
2 -1 3 1 ^1 ^; [' h0 v* S& T4 E
3 1 -5   Q( m1 O' L3 j0 L# ^7 h6 E
4 -1 1
3 D: z( q+ f* ^' D6 ^4 _* _1 3 -13
2 y. ], K3 A& p9 bb=[3 0 3 -6]’;
* S0 B* i& p- c, H2 t# l4 j+ g/ Orank(A) * B1 {! J! y0 d0 S3 U! C
ans= $ s& y) v+ x6 F- l" k' q, S: b
3
2 T. S! E' w) j& K- j2 H, L1 u: Qx1=A\b 9 _, L, w1 r3 j
x1= , o3 `5 L! ^+ t
1.0000   v9 W! l& f3 A7 G+ j% ~
2.0000
: k8 f- d. Y( k1 n1 ]. K3 f1.0000
4 w' ~$ k7 X9 O' }- j  Ux2=pinv(A)*b- m, W5 o. j& _8 J6 F+ j% O
x2=
3 E& [! f, ~# Q1.0000   J( E3 T9 {" r( c8 N  `9 K% t" _- k
2.0000 8 \) R; H7 E5 W1 L
1.0000 7 ^+ Z; G% `2 ?+ u* `9 ~
A*x1-b
) o5 ~$ i; Q- V6 ?1 p. n6 l2 ~5 kans= ; L5 X8 ^+ E. q: O7 r0 l
1.0e-014
9 A; m1 X+ P7 R-0.0888
3 x; \1 F- F# l- l) v-0.0888 0 ?+ S  Y. U% \
-0.1332   Q3 j! A* e) g* L
0
+ k1 g% p7 X6 N) f可见x1并不是方程Ax=b的精确解,用x2=pinv(A)*b所得的解与x1相同。
* m1 ~& ~9 m% Z( M) G0 v: A
8 K# b' y8 ^1 t$ g- Q) n三.欠定方程组
& I+ e! w# |2 b& H# k0 v* n) R欠定方程组未知量个数多于方程个数,但理论上有无穷个解。MATLAB将寻求一个基本解,其中最多只能有m个非零元素。特解由列主元qr分解求得。
. I2 Z* L$ s$ Q2 Q. G【例8】 6 _# p+ v3 @* p/ Q7 Q4 c
解欠定方程组
. f+ S$ ?" a7 {A=[1 -2 1 1;1 -2 1 -1;1 -2 1 5] & [8 R$ i2 ^& d, ^: M. D
A= 0 c0 x3 u. i  K
1 -2 1 1 + b8 _, N# l+ K$ D6 T: ~
1 -2 1 -1 0 V( E  Q& @7 C# |9 o7 @
1 -2 1 -1 6 `6 Y* x1 m0 [4 j: d
1 -2 1 5
) U: \  u9 J3 h4 }b=[1 -1 5]’ 6 ~- J' K( g9 S) l
x1=A\b
& D% \$ P4 J" d  z% bWarning:Rank deficient,rank=2 tol=4.6151e-015 ( D9 N6 b, u( e8 r* K
x1= + w* d" h* G' F& X0 A: l
0
3 e& H; s: P  z& V% ~% C* P9 x" Z-0.0000 0 X: N- ^% i5 Z1 G' B  F, o* r
0 ) ]) l  z7 C# s
1.0000 5 g  D9 h) X. c$ R) z* u; d
x2=pinv(A)*b " D  b' ]; e. U+ X& h! f7 o/ P
x2=
6 N3 ]6 X4 e. _  c0 s0
+ N( k9 k! R0 w( Q-0.0000
" V6 A8 O7 G0 }. e9 \0.0000 + o0 V) [) p; x1 K8 r' Y# W# \9 T
1.0000
4 c; }: @+ w3 s  c. T
9 \$ r  ~" J+ {" L四.方程组的非负最小二乘解 7 P5 Z: x$ @) {2 S, N9 q: Y0 @
在某些条件下,所求的线性方程组的解出现负数是没有意义的。虽然方程组可以得到精确解,但却不能取负值解。在这种情况下,其非负最小二乘解比方程的精确解更有意义。在MATLAB中,求非负最小二乘解常用函数nnls,其调用格式为: " C8 U( U- O7 o1 `: \$ w" K
(1)X=nnls(A,b)返回方程Ax=b的最小二乘解,方程的求解过程被限制在x 的条件下;
% h# ?0 @( A" |1 i(2)X=nnls(A,b,TOL)指定误差TOL来求解,TOL的默认值为TOL=max(size(A))*norm(A,1)*eps,矩阵的-1范数越大,求解的误差越大;
) o- d7 J6 G! s* l9 a& Y: P# Q(3)[X,W]=nnls(A,b) 当x(i)=0时,w(i)<0;当下x(i)>0时,w(i)0,同时返回一个双向量w。 6 A: ^  h4 W! ?' b$ G
【例9】求方程组的非负最小二乘解 0 r2 L( U+ B( |8 A
A=[3.4336 -0.5238 0.6710
$ r4 Y& ?1 X. s3 O( P4 n0 p-0.5238 3.2833 -0.7302 5 p# s* s+ N' z" e7 ?7 N; T5 r% z
0.6710 -0.7302 4.0261];
6 ^3 T4 {& C; u* k  v5 f! hb=[-1.000 1.5000 2.5000]; ) q0 c( a* ]! L! K- u7 \, R
[X,W]=nnls(A,b) 8 F% z# w5 r8 A9 P  A0 F) j
X= 3 J1 \, }# J2 q# h: y! A
0 5 s4 S5 k4 [" c* [4 B; N% x3 d
0.6563
" Q  D; p1 }. }' J" }, k- L0.6998 6 I+ \- A7 r. k  y" q; L
W= . p6 Q- k2 F  a. ]; }
-3.6820
0 K9 y+ |+ G5 `& r9 x0 `  i-0.0000   c: ?2 O0 y* G+ b/ ^
-0.0000 ! R8 E" L# F, F0 c6 l' B8 `  p; A3 K
x1=A\b
+ Q0 K0 }  t4 c. x; d8 k7 B& }9 _x1=
* _! |9 \7 ~, {0 |* F-0.3569 7 u9 Y. t9 x& f- X( D, [, \
0.5744
+ `  S, _* {7 M% o: G8 J' h0.7846 4 {" a2 N* S8 v8 U% ?
A*X-b
; f  E0 Y0 v8 k/ v1 z" l4 g2 Zans= ! l: M! \( m8 o! }. m/ l
1.1258 7 `4 x( a  J! I( y0 B
0.1437 7 g+ V6 }( Y2 q+ Y6 K' X
-0.1616 . e$ T* S2 @; K# i
A*x1-b
* E2 X% ^# B4 S4 uans=
3 h5 u! ~; ^% J' l, V# u1.0e-0.15   B; k1 h1 t4 O9 ]+ d5 C
-0.2220 5 A; D* b+ J5 H% M: p: c% N
0.4441 $ @0 i- C1 _" {4 R5 n
0
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

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

EDA365公众号

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

GMT+8, 2025-7-11 00:15 , Processed in 0.109375 second(s), 23 queries , Gzip On.

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

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

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