计算机应用基础-3-线性与非线性方程(组)求解

1、第三章第三章 代数方程（组）的代数方程（组）的Matlab求解求解3.1 线性方程（组）求解线性方程（组）求解一方程组解个数的判定一方程组解个数的判定二高斯消元法求解二高斯消元法求解三三LU矩阵分解求解矩阵分解求解四四Jacobi迭代求解迭代求解五五 调用内部函数求解调用内部函数求解3.2 非线性方程（组）求解非线性方程（组）求解3.1 线性方程（组）求解线性方程（组）求解一、方程组解个数及求解一、方程组解个数及求解Ax=Bmn2m1mn22221n11211aaaaaaaaaAmp2m1mp22221p11211bbbbbbbbbB对于矩阵对于矩阵3.1 线性方程（组）求解线性方程（组）求解

2、解的判定矩阵解的判定矩阵,又叫增广矩阵又叫增广矩阵mp2m1mmn2m1mp22221n22221p11211n11211bbbaaabbbaaabbbaaac线性方程组解的判定定理线性方程组解的判定定理: (1) 当当 m=n 且且 rank(A)=rank(C)=n 时时 方程组有唯一解方程组有唯一解 x=A-1B x=inv(A)*B (2) 当当rank(A)=rank(C)=rn 时时方程组有无穷多解方程组有无穷多解:为为任任意意数数其其中中rn2 , 1i,axaxaxaxirnrn2211Z=null(A) % 求解求解A矩阵的化零矩阵矩阵的化零矩阵Z=null(A, r) %

3、求解求解A矩阵的化零矩阵的规范形式矩阵的化零矩阵的规范形式若若rank(A)norm(A*x-B)=7.47e-015 1 2 3 4 5 1 4 3 2 1 X= 4 2 1 3 2 4 3 3 4 1 3 2 2 4 A=1 2 3 4; 4 3 2 1; 1 3 2 4; 4 1 3 2; B=5 1; 4 2; 3 3; 2 4; X=inv(A)*BX1 = -9/5, 12/5 28/15, -19/15 58/15, -49/15 -32/15, 41/153.1 线性方程（组）求解线性方程（组）求解精确解：精确解： x1=inv(sym(A)*B检验：检验： norm(doub

4、le(A*x1-B)=0 【例例3-23-2】 A=1 2 3 4; 2 2 1 1; 2 4 6 8; 4 4 2 2 B=1; 3; 2; 6; C=A B; rank(A), rank(C)rank(A)=rank(C)=2Z=null(A,r)Z = 2.0000 3.0000 -2.5000 -3.5000 1.0000 0 0 1.0000X0= 0.9542; 0.7328; -0.0763; -0.2977X0=pinv(A)*B3.1 线性方程（组）求解线性方程（组）求解解出规范化零空间：解出规范化零空间：得出一个特解：得出一个特解：3.1 线性方程（组）求解线性方程（组）求

5、解解出全部解：解出全部解： 2 3 0.9542X = a1 -2.5 + a2 -3.5 + 0.7328 1 0 -0.0763 0 1 -0.2977验证得出的解：验证得出的解： a1=randn(1); a2=rand(1); x=a1*Z(:, 1)+a2*Z(:, 2)+x0 norm(A*x-B)【例例3-33-3】3.1 线性方程（组）求解线性方程（组）求解 1 2 3 4 1 2 2 1 1 X = 2 2 4 6 8 3 4 4 2 2 4 B=1:4; C=A B; rank(A), rank(C)最小二乘解：最小二乘解： x=pinv(A)*B检验：检验： norm(A

6、*x-B)=7.47e-015 二二 高斯消去法求解线性方程组高斯消去法求解线性方程组 高斯消去法它不需要方程组解的初值，也不需要重复迭代计高斯消去法它不需要方程组解的初值，也不需要重复迭代计算，它通过消去和回代两个过程就可以直接求出方程组的算，它通过消去和回代两个过程就可以直接求出方程组的解。和前面的迭代法一样，在计算过程中，也有可能发散解。和前面的迭代法一样，在计算过程中，也有可能发散而得不到方程组的解，但可以通过一些其它方法解决。而得不到方程组的解，但可以通过一些其它方法解决。设一个设一个n n元方程组，以矩阵形式表示为元方程组，以矩阵形式表示为：n21n213n2n1nn22221n1

7、1211bbb xxx aaaaaaaaa3.1 线性方程（组）求解线性方程（组）求解1 消元过程消元过程 aaaaaaaaaA3n12n11n1n212212111n11211111 n121111bbbb（1）令）令 ，进行矩阵的代数运算，可消去第一，进行矩阵的代数运算，可消去第一列的元素（列的元素（2 n) 11111 i1 iaam(2) 同理，可消去第二列的元素（同理，可消去第二列的元素（3 n)（3） 可得任意步消元公式，假设第可得任意步消元公式，假设第k步，有：步，有：3.1 线性方程（组）求解线性方程（组）求解 kkikki1kikkjikkij1kijkkkkikikbmbb

8、n, 2k, 1kj , iamaaaam（4） 经过经过n-1步消元后，系数矩阵称为上三角矩阵步消元后，系数矩阵称为上三角矩阵 a000aaaaa0aaaaAnnnn1n1nn222n , 1n1n2n222221n111n1121111n3.1 线性方程（组）求解线性方程（组）求解 nn2211nbbbb回归过程：回归过程：（1）首先解最后一个方程）首先解最后一个方程（2）依次回代，可得到任意公式）依次回代，可得到任意公式: 1, 3n, 2nk,a/xabxkkkn1kjjkkjkkk3.1 线性方程（组）求解线性方程（组）求解高斯消元法通用程序高斯消元法通用程序function x=U

9、pTriangleFun(A,b)% 求上三角系数矩阵的线性方程组求上三角系数矩阵的线性方程组Ax=b的解的解% A为线性方程组的系数矩阵为线性方程组的系数矩阵% b为线性方程组中的常数向量为线性方程组中的常数向量% x为线性方程组中的解为线性方程组中的解n=size(A);N=n(1);for k=N:-1:1 if(k1) s=A(k,1:(k-1)*x(1:(k-1),1); else s=0; end x(k,1)=(b(k)-s)/A(k,k);end3.1 线性方程（组）求解线性方程（组）求解高斯全主元消去法高斯全主元消去法对于求解精度较高，或方程数量较多，易采用全主元法对于求解精

10、度较高，或方程数量较多，易采用全主元法n21n213n2n1nn22221n11211bbb xxx aaaaaaaaa其增广矩阵为其增广矩阵为(A(1)|b(1),寻找寻找A(1)中绝对值最大为主元素中绝对值最大为主元素 0amaxa1j , inj1ni111j ,1i将此元素进行行列交换，换为第一个元素，再进行消将此元素进行行列交换，换为第一个元素，再进行消元过程元过程3.1 线性方程（组）求解线性方程（组）求解对于任意的第对于任意的第k-1次，其最大的主元素为：次，其最大的主元素为： 0amaxakj , injknikkkj ,ki进过进过n-1次主元消去之后，对后得到上三角矩阵。次

11、主元消去之后，对后得到上三角矩阵。通用算法通用算法3.1 线性方程（组）求解线性方程（组）求解function x,xA= AllGaussFun(A,b) % AllGaussFun.m为用高斯全主元消去法求解线性方程组的解为用高斯全主元消去法求解线性方程组的解% A为线性方程组的系数矩阵为线性方程组的系数矩阵% b为线性方程组的常数向量为线性方程组的常数向量% x为线性方程组的解为线性方程组的解% xA为消元后的系数矩阵为消元后的系数矩阵N = size(A);n = N(1);indexl = 0;indexr = 0;order = 1:n; %记录未知数顺序的向量记录未知数顺序的向量

12、for i=1:(n-1) me = max(max(abs(A(i:n,i:n); %选取全主元选取全主元 for k=i:n for r=i:n if(abs(A(k,r)=me) indexl = k; indexr = r; k=n; break; end end endtemp = A(i,1:n); A(i,1:n) = A(indexl,1:n); A(indexl,1:n) = temp; bb = b(indexl); b(indexl)=b(i); b(i) = bb; %交换主行交换主行 temp = A(1:n,i); A(1:n,i) = A(1:n,indexr);

13、 A(1:n,indexr) = temp; %交换主列交换主列 pos = order(i); order(i) = order(indexr); order(indexr) = pos; %主列的交换会造成未知数顺序的变化主列的交换会造成未知数顺序的变化 for j=(i+1):n if(A(i,i)=0) disp(对角元素为对角元素为0！); return; end l = A(j,i); m = A(i,i); A(j,1:n)=A(j,1:n)-l*A(i,1:n)/m; b(j)=b(j)-l*b(i)/m; endendx=UpTriangleFun(A,b);y=zeros(

14、n,1);for i=1:n for j=1:n if(order(j)=i) y(i)=x(j); end endend x=y;xA =A;例例【3-4】求多元线性方程组的解求多元线性方程组的解16x410 x3x21x3xx56xx2xx34434314321求解程序求解程序“li3_4.m”x =（2.40, -19.20, -1.00, 4.00)3.1 线性方程（组）求解线性方程（组）求解三三 LU矩阵分解求解矩阵分解求解线性方程组变为：线性方程组变为：LUx=b令：令：Y=Ux则：则：LY=b3.1 线性方程（组）求解线性方程（组）求解解两个三角行方程组，解两个三角行方程组，可得

15、解。可得解。A=LU1ll1l1L2n1n21nnn222n11211uuuuuuU其中其中function x=LUfun(A,b)% 实现实现LU分解求线性方程组分解求线性方程组Ax=b的解的解% A为线性方程组的系数矩阵为线性方程组的系数矩阵% b为线性方程组中的常数向量为线性方程组中的常数向量% x为线性方程组中的解为线性方程组中的解r=rank(A);% flag=isexist(A,b) %判断方程组解的情况判断方程组解的情况if flag=0 disp(该方程组无解！该方程组无解！); x=; return;else m,n=size(A); L,U,P=lu(A); %lu函数

16、为函数为MATLAB提供的函数实现矩阵提供的函数实现矩阵LU分分解解 b=P*b; % 解解Ly=b y(1)=b(1); if m1 for k=2:m y(k)=b(k)-L(k,1:k-1)*y(1:k-1); end end y=y; % 解解Ux=y得原方程组的一个特解得原方程组的一个特解 x0(r)=y(r)/U(r,r); if r1 for k=r-1:-1:1; x0(k)=(y(k)-U(k,k+1:r)*x0(k+1:r)/U(k,k); end end x0=x0; % 如果方程组有唯一解如果方程组有唯一解 if flag=1 x=x0; return; else %如

17、果方程组有无穷多解如果方程组有无穷多解 format rat; z=null(A,r); %求出对应齐次方程组的基础解系求出对应齐次方程组的基础解系 mz,nz=size(z); x0(r+1:n)=0; for k=1:nz t=sym(char(107 48+k); K(k)=t; %取取K=k1,k2,.; end x=x0; for k=1:nz x=x+K(k)*z(:,k); %将方程组的通解表示为特解加对应齐次通解形式将方程组的通解表示为特解加对应齐次通解形式 end endend 例例 利用利用LU分解求解如下线性方程组的解分解求解如下线性方程组的解4362x-3522532x

18、1753x43214343214321xxxxxxxxxxx3.1 线性方程（组）求解线性方程（组）求解四四 Jacobi迭代法求解迭代法求解迭代法是对任意给定的初始向量，按照某种方法逐步生迭代法是对任意给定的初始向量，按照某种方法逐步生成近似解序列，是解序列的极限为方程组的解。成近似解序列，是解序列的极限为方程组的解。设线性方程组：设线性方程组： Ax=b若若A非奇异，非奇异，b不等于不等于0，有唯一解，有唯一解x*，上式变换为：，上式变换为：x=Tx+C取一初始向量取一初始向量x(0)为上述方程的近似解为上述方程的近似解x（k+1）=Tx(k)+C满足条件：满足条件： *kkxxlim3.

19、1 线性方程（组）求解线性方程（组）求解此迭代法收敛，此迭代法收敛，x*为原方程组的解。为原方程组的解。Jacobi迭代法迭代法 工程常用此方法解决稀疏矩阵的迭代问题。工程常用此方法解决稀疏矩阵的迭代问题。n21n21nn2n1nn22221n11211bbb xxx aaaaaaaaa系数矩阵为系数矩阵为 aaaaaaaaaAnn2n1nn22221n112113.1 线性方程（组）求解线性方程（组）求解若若A非奇异，方程组非奇异，方程组Ax=b有唯一解有唯一解x=A-1b.方程组通过变换为：方程组通过变换为：nnn1nnn1n,n1nn1nn222n22n2122212111n11n121

20、1121atxaaxaax atxaaxaaxatxaaxaax令令iiiiiiijijatgaab3.1 线性方程（组）求解线性方程（组）求解nnnn21nn2nn222121nn12121gxbxbx gxbxbxgxbxbx aaaDnn2211令令有：有：bDg,ADIADDB1113.1 线性方程（组）求解线性方程（组）求解gBxx算法：算法： （1）设定一个初始向量）设定一个初始向量x(0), 将初始值带入矩阵方程，得到新将初始值带入矩阵方程，得到新向量向量x(1）。）。 （2）将新向量）将新向量x（1）带入矩阵方程，得到另一个向量带入矩阵方程，得到另一个向量x(2)。 (3)对于

21、任意对于任意k步迭代，迭代格式为：步迭代，迭代格式为： gBxxk1k如果序列如果序列x(k)收敛于收敛于x*，则说明，则说明x*为方程组的解为方程组的解Jacobi迭代程序：迭代程序：3.1 线性方程（组）求解线性方程（组）求解function x,n=Jacobifun(A,b,x0,eps,var)% Jacobifun为编写的雅可比迭代函数为编写的雅可比迭代函数% A为线性方程组的系数矩阵为线性方程组的系数矩阵% b为线性方程组的常数向量为线性方程组的常数向量% x0为迭代初始向量为迭代初始向量% eps为解的精度为解的精度% var为迭代步数控制为迭代步数控制% x为线性方程组的解为

22、线性方程组的解% n为求出所需精度的解实际的迭代步数为求出所需精度的解实际的迭代步数if nargin=3 eps= 1.0e-6; M = 200;elseif nargin=eps x0=x; x=B*x0+f; n=n+1; if(n=M) disp(Warning: 迭代次数太多，不收敛！迭代次数太多，不收敛！); return; endend3.1 线性方程（组）求解线性方程（组）求解例例【3-5】用用Jacobi迭代法计算方程组迭代法计算方程组8x5x2x2xx3x212xxx4321321321设初始值为设初始值为x0=0, 0, 0计算程序计算程序“li3_5.m”X=2.39

23、, 2.18, -0.253.1 线性方程（组）求解线性方程（组）求解五五 调用调用Matlab内部函数内部函数5.1 共轭梯度法，通过最优化求解，适用于系数矩阵为共轭梯度法，通过最优化求解，适用于系数矩阵为方阵方阵m, 调用格式如下：调用格式如下：x, flag, relres, iter=bicg(A, b, tol)方方程程组组的的解解方方程程组组成成功功的的标标志志解解的的相相对对误误差差迭迭代代的的次次数数求求解解器器系系数数矩矩阵阵误误差差限限常常数数向向量量3.1 线性方程（组）求解线性方程（组）求解8x5x2x2xx3x212xxx4321321321 A=4, 1, -1;

24、2, -3, 1; 1, 2, -5; b=12, -2, 8; tol=1e-7; x, flag, relres, iter=bicg(A, b, tol)x =2.3929，2.1786 -0.2500flag = 0 relres = 1.1384e-014 iter = 33.1 线性方程（组）求解线性方程（组）求解例例3-6对于方程：对于方程：Ax+xAT=-C可采用可采用Lyapunov方法求解，其调用格式：方法求解，其调用格式： x=lyap(A,C)5.2 lyap函数函数5.3 minres 函数函数 利用残差法求解线性方程组，其调用格式为利用残差法求解线性方程组，其调用格

25、式为x, flag, relres, iter=minres(A, b, tol)3.1 线性方程（组）求解线性方程（组）求解3.2 非线性方程（组）求解非线性方程（组）求解一、对分法一、对分法二、直接迭代法二、直接迭代法三、牛顿迭代法三、牛顿迭代法四、割线法四、割线法五、牛顿迭代法解方程组五、牛顿迭代法解方程组六、符号求解（内部函数）六、符号求解（内部函数）七、非线性方程组的应用七、非线性方程组的应用一一 对分法对分法 对分法或称二分法是求方程近似解的一对分法或称二分法是求方程近似解的一种简单直观的方法。设函数种简单直观的方法。设函数f(x)在在a,b上连续，上连续，且且f(a)f(b)0，

26、则，则f(x)在在a,b上至少有一零点，上至少有一零点，这是微积分中的介值定理，也是使用对分法的这是微积分中的介值定理，也是使用对分法的前提条件。计算中通过对分区间，逐步缩小区前提条件。计算中通过对分区间，逐步缩小区间范围的步骤搜索零点的位置。间范围的步骤搜索零点的位置。 如果我们所要求解的方程从物理意义上来如果我们所要求解的方程从物理意义上来讲确实存在实根，但又不满足讲确实存在实根，但又不满足f(a)f(b)0，这，这时，我们必须通过改变时，我们必须通过改变a和和b的值来满足二分法的值来满足二分法的应用条件。的应用条件。 3.2 非线性方程（组）求解非线性方程（组）求解计算计算f(x)=0的

27、一般计算步骤如下：的一般计算步骤如下： 1、输入求根区间、输入求根区间a,b和误差控制量和误差控制量，定义函数，定义函数f(x)。 2、判断、判断: 如果如果f(a)f(b)0则转下，否则，重新输入则转下，否则，重新输入a和和b的值。的值。 3、计算中点、计算中点 x=(a+b)/2以及以及f(x)的值的值 分情况处理分情况处理（1）|f(x)|：停止计算：停止计算x*=x，转向步骤，转向步骤4（2）f(a)f(x)0：修正区间：修正区间a,xa,b，重复，重复3（3）f(x)f(b)0) t=(a+b)/2;r=compute_bisect(f,t,b,tol); else if(fa*fa

28、b=0) r=(a+b)/2; else if(abs(b-a)clear all;tol=1e-6;g=bisectfun(x3-3*x+1, -1, 1,tol)g = 0.34733.2 非线性方程（组）求解非线性方程（组）求解二二 直接迭代法直接迭代法 对给定的方程对给定的方程f(x)=0，将它转换成等价形式：，将它转换成等价形式： 。给定初。给定初值值x0，由此来构造迭代序列，由此来构造迭代序列 k=1,2,，如果迭代收敛，即，如果迭代收敛，即 有有 ,则就是方程则就是方程f(x)=0的根。在计算中当的根。在计算中当 小于给定的小于给定的精度控制量时，取精度控制量时，取 为方程的根。

29、为方程的根。 对于方程对于方程 构造的多种迭代格式构造的多种迭代格式 ， 怎样判断构怎样判断构 造的迭代格式是否收敛？收敛是否与迭代的初值有关？根据数造的迭代格式是否收敛？收敛是否与迭代的初值有关？根据数学知识，我们可以直接利用以下收敛条件：学知识，我们可以直接利用以下收敛条件： )(xx)x(xk1kbxxkkkk)(limlim1)b(bkkxx11kxb0)(xf3.2 非线性方程（组）求解非线性方程（组）求解 2 在在a,b上可导，并且存在正数上可导，并且存在正数L1，使任意，使任意的的 ，有，有 则在则在a,b上有唯一的点上有唯一的点 满满足足 ， 称称 为为 的不动点。而且迭代格式

30、对任的不动点。而且迭代格式对任意初值均收敛于的不动点，并有下面误差估计式：意初值均收敛于的不动点，并有下面误差估计式：,baxa)(xb)(x, baxLx | )(|1、 当当 有有)(*xx*x)(x|1|01*xxLLxxkk要构造满足收敛条件的等价形式一般比较困难。事实上，要构造满足收敛条件的等价形式一般比较困难。事实上，如果如果 为为f(x)的零点，若能构造等价形式的零点，若能构造等价形式 ，而而 ，由，由 的边疆性，一定存在的邻域的边疆性，一定存在的邻域 ，其，其上有上有 ,这时若初值这时若初值 迭代也就收敛了。由此构造迭代也就收敛了。由此构造收敛迭代格式有两个要素，其一，等价形式

31、收敛迭代格式有两个要素，其一，等价形式 应满足；应满足； 其二，初值必须取自其二，初值必须取自 的充分小邻域，其大小决定的充分小邻域，其大小决定于函数于函数f(x)，及做出的等价形式，及做出的等价形式)(xx*x1| )(|*x)(x,*xx,*0 xxx1| )(| Lx)(xx1| )(|*x*x)(xx3.2 非线性方程（组）求解非线性方程（组）求解上述算法的通用程序为上述算法的通用程序为“fixptfun.m”迭代法的算法：迭代法的算法： 1. 取初始点取初始点x0，最大迭代次数，最大迭代次数N和精度要求和精度要求 ， 设设k=0; 2. 计算计算xk+1= (xk); 3. 若若 ,

32、则停止计算则停止计算 4. 若若k=N, 则停止计算，否则则停止计算，否则k=k+1,转步骤转步骤2k1kxx3.2 非线性方程（组）求解非线性方程（组）求解function x,k,e=fixptfun(f,x0,tol,max1)X(1)=x0;for k=2:max1 X(k)=feval(f,X(k-1); k, e=abs(X(k)-X(k-1) x=X(k); if(etol) break; end if k=max1 disp(超过最大迭代次数！超过最大迭代次数！); endendx=X;3.2 非线性方程（组）求解非线性方程（组）求解例例【3-8】利用迭代法求非线性方程组在利用

33、迭代法求非线性方程组在x0=0.4的近似的近似解，误差为解，误差为1e-6，最大迭代次数为，最大迭代次数为30。0 xxsin采用程序采用程序”ME_31_1.m”,计算计算X=0.8767,迭代次数为迭代次数为12.3.2 非线性方程（组）求解非线性方程（组）求解三三 牛顿迭代法牛顿迭代法对方程对方程f(x)=0可构造多种迭代格式可构造多种迭代格式 ，牛顿迭代法是，牛顿迭代法是借助于对函数借助于对函数f(x)=0的泰勒展开而得到的一种迭代格式。的泰勒展开而得到的一种迭代格式。 将将f(x)=0在初始值在初始值x0做泰勒展开得：做泰勒展开得：)(1kkxx)(! 2)()()()(00 000

34、 xxxfxxxfxfxf取展开式的线性部分作为的近似值，则有：取展开式的线性部分作为的近似值，则有： 0 xf)xx)(x(f)x(f000设设f(x0)K0,)()(0001xfxfxx类似地，再将类似地，再将f(x)=0在在x1作泰勒展开并取其线性部分得到作泰勒展开并取其线性部分得到：)()(1112xfxfxx3.2 非线性方程（组）求解非线性方程（组）求解得到牛顿法的迭代格式得到牛顿法的迭代格式, 2 , 1,)()(1kxfxfxxkkkk牛顿法的算法：牛顿法的算法： 1. 取初始点取初始点x0，最大迭代次数，最大迭代次数N和精度要求和精度要求 ，设，设k=0; 2. 如果如果f(

35、x0)=0，则停止计算，否则计算：，则停止计算，否则计算： 3. 若若 ,则停止计算则停止计算 4. 若若k=N, 则停止计算，否则则停止计算，否则k=k+1,转步骤转步骤2k1kxx)x(f)x(fxxkkk1k3.2 非线性方程（组）求解非线性方程（组）求解function x,e,k,y=newtonfun(f,df,x0,delta,max1)% newtonfun为编写利用牛顿法求解非线性方程组的根的函数为编写利用牛顿法求解非线性方程组的根的函数% f为非线性函数；为非线性函数； % df为为f的微商；的微商；% x0为初始值为初始值% delta为给定允许的误差；为给定允许的误差；

36、 % max1为迭代的最大次数为迭代的最大次数% x为用牛顿法求得的非线性方程的近似解为用牛顿法求得的非线性方程的近似解% e为为x0的误差估计；的误差估计； % k为迭代次数；为迭代次数；% y=f(x0)x0,feval(f,x0)for k=1:max1 x1=x0-feval(f,x0)/feval(df,x0); e=abs(x1-x0); x0=x1; x1,e,k,y=feval(f,x1) if(etol) n=n+1; x1=r; fx=subs(sym(f),x1); s=fx*fa; if(s0) r=b-(x1-b)*fb/(fx-fb); else r=a-(x1-a

37、)*fa/(fx-fa); end err=abs(r-x1); endend% secantfun.m为编写的割线法为编写的割线法计算非线性方程的函数计算非线性方程的函数% f为非线性方程为非线性方程% a为非线性方程的左区间为非线性方程的左区间% b为非线性方程的右区间为非线性方程的右区间% tol为误差容限为误差容限,当误差容限缺当误差容限缺省时省时 默认为默认为10的的-5次方次方% r为求得的非线性方程的近似为求得的非线性方程的近似根根 % n为迭代的次数为迭代的次数3.2 非线性方程（组）求解非线性方程（组）求解例例【3-10】利用割线法求解非线性方程在利用割线法求解非线性方程在【

38、-1， 2】上的根，精度为上的根，精度为1e-6。06x8e2x2计算程序计算程序“ME_31_3.m”，结果：，结果： r=1.6731 n=183.2 非线性方程（组）求解非线性方程（组）求解非线性方程组解法大多来自于非线性方程，只不过计算非线性方程组解法大多来自于非线性方程，只不过计算换成了矩阵，其中最重要的就是雅克比矩阵。非线性方换成了矩阵，其中最重要的就是雅克比矩阵。非线性方程组一般形式如下：程组一般形式如下：0)x,x,x(f0)x,x,x(f0)x,x,x(fn21nn212n211五五 迭代求解非线性方程组迭代求解非线性方程组为了方便，记为了方便，记Tn21Tn21)f,f,f

39、(F)x,x,x(x则，非线性方程组为则，非线性方程组为 F(x)=03.2 非线性方程（组）求解非线性方程（组）求解 非线性方程组可构造等价形式非线性方程组可构造等价形式)x,x,x(x)x,x,x(x)x,x,x(xn21nnn2122n2111设上述函数及其一阶导数在设上述函数及其一阶导数在x*区域连续，且满足区域连续，且满足（1）初始值）初始值x(0)足够接近足够接近x*;(2) 对于任意在对于任意在x*领域的领域的x点，点， 为为 F(x)的的Jacobi矩阵，即矩阵，即 1xF xF3.2 非线性方程（组）求解非线性方程（组）求解 nn1nn111xfxfxfxfxF则迭代序列则迭

40、代序列xk均收敛到均收敛到x*5.1 牛顿迭代法牛顿迭代法3.2 非线性方程（组）求解非线性方程（组）求解假定假定x*为方程组的近似解，牛顿迭代法是借助于多元函数的泰为方程组的近似解，牛顿迭代法是借助于多元函数的泰勒展开而得到的一种迭代格式。勒展开而得到的一种迭代格式。 假设它的假设它的k次近似解为次近似解为x(k)，其，其泰勒展开得：泰勒展开得： 0 xxRlim0 xxRxxxFxFxF0 xk*k*kk*其中其中当当x(k)接近性接近性x*，可略去，可略去 从而得到线性方程从而得到线性方程 k*xxR kkxFxxF3.2 非线性方程（组）求解非线性方程（组）求解上述方程解上述方程解 x

41、近似误差近似误差 k*xx 1 , 0kxFxFxxxxxxxk1kkk1kk1k有有记记实际迭代时，采用如下的方法实际迭代时，采用如下的方法 kkkk1kxFxxFxxx3.2 非线性方程（组）求解非线性方程（组）求解牛顿迭代法的算法：牛顿迭代法的算法： 1. 取初始点取初始点x0，最大迭代次数，最大迭代次数N和精度要求和精度要求 ，设，设k=0; 2. 计算计算F(x(k)和和F(x(k）)； 3. 若若Jacobi矩阵矩阵F(x(k）)奇异，则停止计算奇异，则停止计算 ,否则计算否则计算 求得求得 x（k);令令 4若若 ，停止计算，停止计算 4. 若若k=N, 则停止计算，否则则停止计

42、算，否则k=k+1,转步骤转步骤2k1kxx kkxFxxF k1kxxx3.2 非线性方程（组）求解非线性方程（组）求解function x,n,err=NewtonGroudfun(f,jf,x0,tol,max1)x=x0;y=feval(f,x0);for k=1:max1 j=feval(jf,x0); q=x0-(jy); y=feval(f,q); err=norm(q-x0); x0=q; x=x,x0; n=k; if(erreps fx = subs(F,findsym(F),x0); J = zeros(m,m); for i=1:m x1 = x0; x1(i) = x

43、1(i)+h; J(:,i) = (subs(F,findsym(F),x1)-fx)/h; end % FastDownGroudfun.m为编写利用最速下降为编写利用最速下降算法求非线性方程组的解的函数算法求非线性方程组的解的函数% F为非线性方程组为非线性方程组% x0为给定的初始值为给定的初始值% h数值微分增量步数值微分增量步% eps为解的精度为解的精度% r为求得非线性方程组的一组解为求得非线性方程组的一组解% n为迭代步数为迭代步数3.2 非线性方程（组）求解非线性方程（组）求解lamda = fx/sum(diag(transpose(J)*J); r=x0-J*lamda;

44、 fr = subs(F,findsym(F),r); tol=dot(fr,fr); x0 = r; n=n+1; if(n100000) disp(迭代步数太多，可能不收敛！迭代步数太多，可能不收敛！); return; endendformat short;3.2 非线性方程（组）求解非线性方程（组）求解六六 Matlab内部函数（符号函数）内部函数（符号函数）1. 内部符号函数内部符号函数fzero()可以求解非线性方程，其调用可以求解非线性方程，其调用格式如下：格式如下：2, 1, 0,zeroppoptionxfunfx 例例【3-11】求解非线性方程在求解非线性方程在x0=1 时解时解求解程序求解程序”ME_31_9.m” X=-3.043.2 非线性方程（组）求解非线性方程（组）求解3sinexx2. 内部符号函数内部符号函数fsolve()可以求解非线性方程组，采可以求解非线性方程组，采用的是最小二乘法，其调用格式如下：用的是最小二乘法，其调用格式如下：2p, 1p,option, 0 x,funfslovex 例例【3-11】求解非线性方程组在求解非线性方程组在x0=1 1时解时解0 xx3xx20 xxx4x2122221211求解程序求解程序”ME_31_6.m” X=0.25 03.2 非线性方程（组）求解非线性方程（组）求解3.2 非线性方