matlab求解代数方程组解析

1、第三讲Matlab求解代数方程组理论介绍：直接法+迭代法，简单介绍相关知识和应用条件及注意事项软件求解：各种求解程序讨论如下表示含有n个未知数、由n个方程构成的线性方程组:aiiXi - ai2X2ainXn 二匕a21X1a22X2a2nxn = b2(D、直接法1.高斯消元法:aniXi . an2X2, annXn =bn高斯消元法的基本原理:在(1)中设an =0,将第一行乘以-a1,加到第k(k =2,3,，n),得: ana% +ai(?X2 +ai(?Xn =bia22)Xii 十a(1n)Xn=b22)(2)aS)X2- - an2) Xn = bn2)其中ag =源,“(1)

2、 =b再设a22) #0,将(2)式的第二行乘以号(k = 31，n)加 a22到第k行，如此进行下去最终得到:a(1)Xiai(2)X2露 Xn =b(1)a22)X,或 Xn=b22)(3)(n劣(n(nan,njXnan J,n Xn - bn 4c(n).(n)ann Xn - bnA，并如此进行从(3)式最后一个方程解出代入它上面的一个方程解出下去，即可依次将Xn，Xn4,，X2，x全部解出，这样在ak?¥0(k=1,2,，n)的假设下，由上而下的消元由下而上的回代，构成了方程组的高斯消元法 高斯消元法的矩阵表示: 若记 A = (q)n*,x=(Xi，Xn)T,b=(bi

3、,，bn)T ,则(1)式可表为 Ax = b.于是高斯消元法的过程可用矩阵表示为:Mn1 M2M1Ax =MnM2M1b.其中:M1a a21 尹 alln1an高斯消元法的Matlab程序：%顺序 gauss消去法,gauss函数 functionA,u=gauss(a,n)for k=1:n-1%肖去过程for i=k+1:nfor j=k+1:n+1%如果a(k,k)=0,则不能削去if abs(a(k,k)>1e-6%计算第k步的增广矩阵a(i,j)=a(i,j)-a(i,k)/a(k,k)*a(k,j);else% a(k,k)=0,顺序gauss消去失败disp('

4、顺序gauss消去失败)； pause;exit;endendendend%回代过程x(n)=a(n,n+1)/a(n,n);for i=n-1:-1:1s=0;for j=i+1:ns=s+a(i,j)*x(j);endx(i)=(a(i,n+1)-s)/a(i,i);end%返回gauss消去后的增广矩阵A=triu(a);%返回方程组的解u=x;练习和分析与思考：用高斯消元法解方程组：x1 +2x2 +4x3 +x4 +7x5 =152x1 +3x2 +x4 +8x5 =144x1 x2 7x3 6x4 x5 = 19x1 +x2 +2x4 +x5 =5x1 3x2 x4 x5 = 62

5、.列主元素消元法在高斯消元法中进行到第k步时，不论aikk)是否为0,都按列选择| aikk) | (i =k,，n)中最大的一个，称为列主元，将列主元所在行与第k行交换再按高斯消元法进行下去称为列主元素消元法。列主元素消元法的matlab程序痛!J主元guass消去函数functionA,u=gauss(a,n)%肖去过程for k=1:n-1%选主元c=0;for q=k:nif abs(a(q,k)>cc=a(q,k);l=q；endend呦口果主元为0,则矩阵A不可逆if abs(c)<1e-10disp( 'error ');pause;exitend冽果

6、l不等于k,则交换第l行和第k行if l=kfor q=k:n+1temp=a(k,q);a(k,q)=a(l,q);a(l,q)=temp;endend%计算第k步的元素值for i=k+1:nfor j=k+1:na(i,j)=a(i,j)-a(i,k)/a(k,k)*a(k,j);endendend%回代过程x(n)=a(n,n+1)/a(n,n);s=0;for j=i+1:ns=a+a(i,j)*x(j);endx(i)=(a(i,n+1)-s)/a(i,i);end%返回列主元gauss消去后的增广矩阵A=triu(a);%返回方程组的解u=x;练习和分析与思考：用列主元消去法重新

7、求解gauss消元法求解的上述问题。二、迭代法1 .迭代法的总体思想(1)迭代公式的构造：对线性方程组Ax = b ,可以构造迭代公式X "刊=BX+ f ,给出X(0),由迭代公式的X(k),如果X(k)收敛于X*,则X*就是原方程组的解。(2)矩阵的分解：设线性方程组Ax=b,其中A非奇异，则可以把A矩阵分解：A=D-L-U ,D =diaga11,a22,，ann,Z0a210L = -a31%20999pman 2an30a1200；al3a1na23a2n .mm0an 口 ,n0 于是Ax = b化为x = D(L +U )x + D "b,与之对应的迭代公式为

8、:x(k 1) =d'(L U)x(k) D-1b.2 .雅可比(Jacobia。迭代法公式:B =D(L U), f =D,b用A矩阵的元素表示为:nbi -%aijXj)x(k 1),(i =1,2, ,n)_ j :,jiaJacobian迭代法的Matlab程序function y,n=jacobi(A,b,x0,eps)%K差if nargin=3eps=1.0e-6;elseif nargin<3error returnend%求A的对角矩阵，下三角阵，上三角阵D=diag(A);diag(diag(A)?L=-tril(A,-1);U=-triu(A,1);B=D(

9、L+U);f=Db;y=B*x0+f;%迭代的次数n=1;%当误差没有满足要求时继续迭代while norm(y-x0)>=epsx0=y;y=B*x0+f;n=n+1;end练习和分析与思考：利用Jacobian迭代法求解方程组：x1 +3x3 =12x2 + x3 +2x4 = 63x1 x2 15x3 = 52x2 4x4 =83 .高斯-塞德尔(Gauss-Seide)迭代法公式：将 Jacobi迭代公式 Dx(k*)=(L +U )x(k) +b 改进为 Dx(k也=Lx(k川 +Ux(k) +b ,于是得到 x(k 1) =(D-L)Ux(k) (D-L)b.B =(D -L

10、),U, f =(D -L)b用A矩阵的元素表示为:i -1nbiaijxjk) %ajk),(i =12 ,n)_ j -j - iaHGauss-Seidel迭代法的Matlab程序function y,n=gauseidel(A,b,x0,eps)%I差if nargin=3 eps=1.0e-6;elseif nargin<3 error returnend%求A的对角矩阵，下三角阵，上三角阵D=diag(A);diag(diag(A)?L=-tril(A,-1);U=-triu(A,1);G=(D-L)U;f=(D-L)b;y=G*x0+f;名迭代的次数n=1;%当误差没有满足

11、要求时继续迭代while norm(y-x0)>=epsx0=y;y=G*x0+f;n=n+1;end4 .超松弛(Successive Over Relaxation method 迭代法(SOR)三、求解非线性方程的方法1 .根的隔离二分法利用函数的性质确定根的大致范围(a,b，取(a,b)的中点x0=g;辿，若2f (x0) = 0 ,则%即是根，否则如果f (a) - f (x0) <0,令a1 = a, h = %,如果f ( b) f (ox ) 金a=x0, b =b,则在(a, b)内至少有一根，且(a,b)二(a ,bi),再取区匕)的中点与=亘也，如此进行下去，

12、包含根的区间(anh)(n=1,2，)的长 2度每次缩小一半，n足够大时即可获得满意精度。2 .切线法对于方程f(x)=0将f(x)在“作Taylor展式保留线性项，即_ _ _ 'f (x) = f (xk) f (xk)(x -xk),设f'(xk)¥0,然后用xk4代替右端的x就得到迭代公式：xr=xk-f这种 f (xk)方法称为牛顿(Newton)切线法。3 .割线法在牛顿切线法中用差商f(xk) f(x")代替f(xk),则xk - xk 4xkxk - f(xk)(xk -xk”)(割线法)f (xk) - f (xy)练习和分析与思考：求方程

13、f (x) = x2+x 14的一个正根。四、求解非线性方程组的牛顿法方程组记作F(X)= 0,其中X =(为；凶T F X十f( X( f)2X),fn X '(诩X(k) =(Xi(k)，Xnk)T,且方程组F(x)=0的第k步近似解，与牛顿切线法类似，在。°作Taylor展式，线性化后用x(k41)代替x可得:2 ,n)开1 ：Xn 寸2 ：Xnfi(x(k 1)="产)(研 1)”的) ；X1，f fexCX2身2身2(F (x) =CX1CX2-则有 F(x(k*) =F(x(k) +F'(x(k)(x(k.)-x(k),若 F'(Xk )

14、可逆，我们得到解非线性方程组的牛顿迭代公式：x(k1) =x(k) _F,(x(k)/F(x(k).注：实际上在计算过程的第k步，往往先计算F(x(k)和F'(x(k),再解方程 组F'(x(k)Vx(k) =-F(x(k)得到Vx后，令x(k+)=x(k) +Vx(k)即可.迭代过程需要分 析其收敛性、分支与混沌.练习和分析与思考：用牛顿迭代法解非线性方程组：X2 - 10x1 , X； x3 7=0 224x1 x2 + x3 - 2x3 = 02,2 人 ，2-x1 +x2 -3x2 +x3 = 0五、Matlab求解(非)线性方程组的内置工具箱1 .线性方程组求解(1)

15、直接求解：线性方程组 Ax = b,解法x=A b,(注意：此处 不是/符号) 若A为方阵，x = inv(A)*b,若A不是方阵，x=pinv(A)*b.练习和分析与思考：求方程组的解.5x1 +6x2 = 1 x1 +5x2 +6x3 = 0 « x2 +5x3 +6x4 = 0 x3 + 5x4 + 6x5 = 0 x4 5x5 =1(2)利用矩阵分解求解线性方程组矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的 乘积.常见的矩阵分解有LU分解、QR分解、Cholesky分解以及Schur分解、 HessenberS解、奇异分解等.LU分解：矩阵的LU分解就是将

16、一个矩阵表示为一个交换下三角矩阵和一个 上三角矩阵的乘积形式.只要方阵A非奇异，矩阵的LU分解总是可以进行的.Matlab提供的lu函数用于对矩阵进行LU分解，其调用格式为：L,U=lu(A):产生一个上三角阵U和一个变换形式的下三角矩阵 L(行交换) 使之满足A=LU ,注意：这里的矩阵A必须是方阵.L,U,P=lu(A):产生一个上三角阵U和一个变换形式的下三角矩阵 L以及 一个置换矩阵P使之满足PA=LU,注意：这里的矩阵A必须是方阵.实现LU分解后，线性方程组 Ax = b的解为：x = U (L b)或x=U (L Pb).QR分解：矩阵的QR分解就是将一个矩阵A分解成一个正交矩阵Q

17、和一个 上三角矩阵R的乘积形式.QR分解只能对方阵进行.Matlab的函数qr用于对矩 阵进行QR分解，其调用格式为：Q,R=qr(A):产生一个正交矩阵Q和一个上三角矩阵R使之满足A=QR.Q,R,E=qr(A):产生一个正交矩阵Q和一个上三角矩阵R以及一个置换矩 阵E使之满足AE=QR方阵.实现QR分解后，线性方程组人*4的解为：x=R (Q b)或x = E(R (Q b).Cholesky分解：如果矩阵A是对称正定的，贝U Cholesky分解将矩阵A分解成 一个下三角矩阵R和上三角矩阵的乘积 R'(R的转置)，即A= R' R. Matlab的函数 chol(A)用于

18、对矩阵进行Cholesky分解，其调用格式为：R=chol(A):产生一个上三角矩阵 R使之满足A=R 'R.若矩阵A不是对称正定 的，则输出一个出错信息.R,p=chol(A):这个命令格式不输出出错信息.当A是对称正定的，则p=0, R 与上述格式得到的结果相同，否则 p为一个正整数，如果 A为满秩矩阵，则R 为一个阶数为q=p-1的上三角阵，且满足R'R=A(1:q,1:q).实现 Cholesky 分解后，Ax = b 变为 R'Rx = b ,所以 x = R (R' b).2 .单变量非线性方程求解Matlab的函数fzero可用于对矩阵求单变量非线

19、性方程的根，其调用格式为： z=fzero( 'fname' ,x0,tol,trace)其中fname是待求根的函数文件名，x0为搜索起点，一个函数可能有多个 根，但fzero函数只给出离x0最近的那个根.tol控制结果的相对精度，缺省时取 tol=eps,trace指定迭代信息是否在运算中显示，为 1时显示，为0时不显示，缺 省时 trace=0.练习和分析与思考：求 f(x)=x-10x +2在x0=0.5 附近的根.3 .非线性方程组的求解对于非线Tt方程组F (x) = 0用fsolve函数求其数值解.fsolve函数的调用格式为：x=fsolve( Un'

20、,xption)其中x为返回的解，fun用于定义需求解的非线性方程组的函数文件名，x0是求根过程的初值，option为最优化工具箱的选项设定.最优化工具箱提供了 20 多个选项，用户可以使用optimset命令将它们显示出来.如果想改变其中某个选 项，可以调用 optimset()函数来完成.例如 optimset( 'display ' ., ' off ') 六、实例赏析a=1 -1 4 -2; 1 -1 -1 2;3 1 7 -2;1 -3 -126; rref(a)b=2 3 1;1 -2 4;3 8 -2;4 -1 9b1=4 -5 13 -6'b,b1rref(b,b1)x,y=solve('5*yA2-6*x*y+5*xA2-16','yA2-x*y+2*xA2-y-x-4') a*xA2+b*x+c=0求方程 xA4+7xA3 +9x-20=0的全部根a=1,1.5,2,9,7 ; 0,3.6,0.5,-4,4 ; 7,10,-3,22,33; 3,7,8.5,21,6; 3,8,0,90,-20;b=3;-