用matlab求解非线性方程组的几种方法之程序

第二章非线性方程（组）的数值解法的MATLAB程序9．高等教育出版社教育电子音像出版社作者：任玉杰本章主要介绍方程根的有关概念，求方程根的步骤，确定根的初始近似值的方法（作图法，逐步搜索法等），求根的方法（二分法，迭代法，牛顿法，割线法，米勒（Müller）法和迭代法的加速等）及其MATLAB程序，求解非线性方程组的方法及其MATLAB程序.2.12.12.12.1方程方程方程方程（（（（组组组组））））的根及其的根及其的根及其的根及其MATLABMATLABMATLABMATLAB命令命令命令命令2.1.22.1.22.1.22.1.2求解方程求解方程求解方程求解方程（（（（组组组组））））的的的的solve命令命令命令命令求方程f(x)=q(x)的根可以用MATLAB命令:x=solve('方程f(x)=q(x)',’待求符号变量x’)求方程组fi(x1,…,xn)=qi(x1,…,xn)(i=1,2,…,n)的根可以用MATLAB命令:E1=sym('方程f1(x1,…,xn)=q1(x1,…,xn)');…………………………………………………….En=sym('方程fn(x1,…,xn)=qn(x1,…,xn)');[x1,x2,…,xn]=solve(E1,E2,…,En,x1,…,xn)2.1.32.1.32.1.32.1.3求解求解求解求解多项式多项式多项式多项式方程方程方程方程（（（（组组组组））））的的的的roots命令命令命令命令如果)(xf为多项式，则可分别用如下命令求方程0)(=xf的根，或求导数)('xf（见表2-1）.表2-1求解多项式方程（组）的roots命令命令功能xk=roots(fa)输入多项式)(xf的系数fa（按降幂排列），运行后输出xk为0)(=xf的全部根.dfa=polyder(fa)输入多项式)(xf的系数fa（按降幂排列），运行后输出dfa为多项式)(xf的导数)('xf的系数.dfx=poly2sym(dfa)输入多项式)(xf的导数)('xf的系数dfa（按降幂排列），运行后输出dfx为多项式)(xf的导数)('xf.2.1.42.1.42.1.42.1.4求解方程求解方程求解方程求解方程（（（（组组组组））））的的的的fsolve命令命令命令命令如果非线性方程（组）是多项式形式，求这样方程（组）的数值解可以直接调用上面已经介绍过的roots命令.如果非线性方程（组）是含有超越函数，则无法使用roots命令，需要调用MATLAB系统中提供的另一个程序fsolve来求解.当然，程序fsolve也可以用于多项式方程（组），但是它的计算量明显比roots命令的大.fsolve命令使用最小二乘法（leastsquaresmethod）解非线性方程（组）(FX=)0的数值解，其中X和F（X）可以是向量或矩阵.此种方法需要尝试着输入解X的初始值（向量或矩阵）X0，即使程序中的迭代序列收敛，也不一定收敛到(FX=)0的根（见例2.1.8）.fsolve的调用格式的调用格式的调用格式的调用格式：X=fsolve(F,X0)第二章第二章第二章第二章非线性方程非线性方程非线性方程非线性方程（（（（组组组组））））的数值解法的数值解法的数值解法的数值解法第二章非线性方程（组）的数值解法的MATLAB程序10．高等教育出版社教育电子音像出版社作者：任玉杰输入函数)(xF的M文件名和解X的初始值（向量或矩阵）X0，尝试着解方程（组）(FX=)0，运行后输出(FX=)0解的估计值（向量或矩阵）X.要了解更多的调用格式和功能请输入:helpfsolve,查看说明.2.22.22.22.2搜索根的方法及其搜索根的方法及其搜索根的方法及其搜索根的方法及其MATLABMATLABMATLABMATLAB程序程序程序程序求解非线性方程根的近似值时，首先需要判断方程有没有根？如果有根，有几个根？如果有根，需要搜索根所在的区间或确定根的初始近似值（简称初始值）.搜索根的近似位置的常用方法有三种:作图法、逐步搜索法和二分法等，使用这些方法的前提是高等数学中的零点定理.2.2.12.2.12.2.12.2.1作图法及其作图法及其作图法及其作图法及其MATLABMATLABMATLABMATLAB程序程序程序程序作函数的图形的方法很多，如用计算机软件的图形功能画图，或用高等数学中应用导数作图，或用初等数学的函数叠加法作图等.下面介绍两种作图程序.作函数作函数作函数作函数)(xfy====在区间在区间在区间在区间[[[[a,ba,ba,ba,b]]]]的图形的的图形的的图形的的图形的MATLABMATLABMATLABMATLAB程序一程序一程序一程序一x=a:h:b;%h是步长y=f(x);plot(x,y)grid,gtext('y=f(x)')说明：⑴此程序在MATLAB的工作区输入，运行后即可出现函数)(xfy=的图形.此图形与x轴交点的横坐标即为所要求的根的近似值.⑵区间[a,b]的两个端点的距离b-a和步长h的绝对值越小，图形越精确.作函数作函数作函数作函数)(xfy====在区间在区间在区间在区间[[[[a,ba,ba,ba,b]]]]上的图形的上的图形的上的图形的上的图形的MATLABMATLABMATLABMATLAB程序二程序二程序二程序二将)(xfy=化为)()(xgxh=，其中)()(xgxh和是两个相等的简单函数x=a:h:b;y1=h(x);y2=g(x);plot(x,y1,x,y2)grid,gtext('y1=h(x),y2=g(x)')说明：此程序在MATLAB的工作区输入，运行后即可出现函数)()(21xgyxhy==和的图形.两图形交点的横坐标即为所要求的根的近似值.下面举例说明如何用计算机软件MATLAB的图形功能作图.2.2.22.2.22.2.22.2.2逐步搜索法逐步搜索法逐步搜索法逐步搜索法及其及其及其及其MATLABMATLABMATLABMATLAB程序程序程序程序逐步搜索法也称试算法.它是求方程0)(=xf根的近似值位置的一种常用的方法.逐步搜索法依赖于寻找连续函数)(xf满足)(af与)(bf异号的区间],[ba.一旦找到区间，无论区间多大，通过某种方法总会找到一个根.MATLAB的库函数中没有逐步搜索法的程序，现提供根据逐步搜索法的计算步骤和它的收敛判定准则编写其主程序，命名为zhubuss.m.逐步搜索法的逐步搜索法的逐步搜索法的逐步搜索法的MATLABMATLABMATLABMATLAB主程序主程序主程序主程序输入区间端点a和b的值，步长h和精度tol，运行后输出迭代次数k=(b-a)/h+1,方程0)(=xf根的近似值r.function[k,r]=zhubuss(a,b,h,tol)%输入的量---a和b是闭区间[a,b]的左、右端点；%---h是步长；%---tol是预先给定的精度.%运行后输出的量---k是搜索点的个数；%---r是方程在[a,b]上的实根的近似值，其精度是tol；X=a:h:b;Y=funs(X);n=(b-a)/h+1;m=0;X(n+1)=X(n);Y(n+1)=Y(n);第二章非线性方程（组）的数值解法的MATLAB程序11．高等教育出版社教育电子音像出版社作者：任玉杰fork=2:nX(k)=a+k*h;Y(k)=funs(X(k));%程序中调用的funs.m为函数sk=Y(k)*Y(k-1);ifsk=0,m=m+1;r(m)=X(k);endxielv=(Y(k+1)-Y(k))*(Y(k)-Y(k-1));if(abs(Y(k))tol)&(xielv=0)m=m+1;r(m)=X(k);endend例例例例2.2.42.2.42.2.42.2.4用逐步搜索法的MATLAB程序分别求方程0332223=--+xxx和0)2sin(cos3=x在区间]2,2[-上的根的近似值，要求精度是0.0001.解解解解在MATLAB工作窗口输入如下程序[k,r]=zhubuss(-2,2,0.001,0.0001)运行后输出的结果k=4001r=-1.2240-1.0000-1.0000-0.99901.2250即搜索点的个数为k=4001，其中有5个是方程⑴的近似根，即r=-1.2240，-1.0000，-1.0000，-0.9990，1.2250，其精度为0.0001.在程序中将y=2.*x.^3+2.*x.^2-3.*x-3用y=sin(cos(2.*x.^3))代替，可得到方程⑵在区间]2,2[-上的根的近似值如下r=-1.9190-1.7640-1.5770-1.3300-0.92200.92301.33101.57801.76501.9200如果读者分别将方程⑴的结果与例2.2.3比较，方程⑵的结果与例2.1.2比较，将会发现逐步搜索法的MATLAB程序的优点.如果精度要求比较高，用这种逐步搜索法是不合算的.2.32.32.32.3二分法及其二分法及其二分法及其二分法及其MATLABMATLABMATLABMATLAB程序程序程序程序2.3.12.3.12.3.12.3.1二分法二分法二分法二分法二分法也称逐次分半法逐次分半法逐次分半法逐次分半法.它的基本思想是：先确定方程0)(=xf含根的区间（a,b），再把区间逐次二等分.我们可以根据式（2.3b）、区间[a,b]和误差ε，编写二分法求方程根的迭代次数的MATLAB命令.已知闭区间[a,b]和误差ε.用二分法求方程误差不大于ε的根的迭代次数k的MATLAB命令k=-1+ceil((log(b-a)-log(abtol))/log(2))%ceil是向+∞方向取整，abtol是误差ε.2.3.22.3.22.3.22.3.2二分法的二分法的二分法的二分法的MATLABMATLABMATLABMATLAB程序程序程序程序二分法需自行编制程序，现提供用二分法求方程f(x)=0的根*x的近似值kx的步骤和式（2.3a）编写一个名为erfen.m的二分法的MATLAB主程序如下.二分法的二分法的二分法的二分法的MATLABMATLABMATLABMATLAB主程序主程序主程序主程序求解方程0)(=xf在开区间（a,b）内的一个根的前提条件是)(xf在闭区间[a,b]上连续，且0)()(⋅bfaf.输入的量：a和b是闭区间[a,b]的左、右端点，abtol是预先给定的精度.运行后输出的量：k是使用二分法的次数.x是方程在（a,b）内的实根x*的近似值，第二章非线性方程（组）的数值解法的MATLAB程序12．高等教育出版社教育电子音像出版社作者：任玉杰其精度是abtol.wuca=|bk-ak|/2是使用k次二分法所得到的小区间的长度的一半，即实根x*的近似值x的绝对精度限,满足wuca≤abtol.yx=f(xk)，即方程f(x)=0在实根x*的近似值x处的函数值.function[k,x,wuca,yx]=erfen(a,b,abtol)a(1)=a;b(1)=b;ya=fun(a(1));yb=fun(b(1));%程序中调用的fun.m为函数ifya*yb0,disp('注意：ya*yb0,请重新调整区间端点a和b.'),returnendmax1=-1+ceil((log(b-a)-log(abtol))/log(2));%ceil是向∞+方向取整fork=1:max1+1a;ya=fun(a);b;yb=fun(b);x=(a+b)/2;yx=fun(x);wuca=abs(b-a)/2;k=k-1;[k,a,b,x,wuca,ya,yb,yx]ifyx==0a=x;b=x;elseifyb*yx0b=x;yb=yx;elsea=x;ya=yx;endifb-aabtol,return,endendk=max1;x;wuca;yx=fun(x);例例例例2.3.2.3.2.3.2.3.3333确定方程x3-x+4=0的实根