3-7变量非线性方程组的逆Broyden解法matlab程序

B r o y d e n方法求解非线性方程组的M a t l a b实现-CAL-FENGHAI.-(YICAI)-Company One1Broyden方法求解非线性方程组的Matlab实现注：matlab代码来自网络，仅供学习参考。

1.把以下代码复制在一个.m文件上function [sol, it_hist, ierr] = brsola(x,f,tol, parms)% Broyden's Method solver, globally convergent% solver for f(x) = 0, Armijo rule, one vector storage%% This code comes with no guarantee or warranty of any kind.%% function [sol, it_hist, ierr] = brsola(x,f,tol,parms)%% inputs:% initial iterate = x% function = f% tol = [atol, rtol] relative/absolute% error tolerances for the nonlinear iteration% parms = [maxit, maxdim]% maxit = maxmium number of nonlinear iterations% default = 40% maxdim = maximum number of Broyden iterations% before restart, so maxdim-1 vectors are% stored% default = 40%% output:% sol = solution% it_hist(maxit,3) = scaled l2 norms of nonlinear residuals% for the iteration, number function evaluations,% and number of steplength reductions% ierr = 0 upon successful termination% ierr = 1 if after maxit iterations% the termination criterion is not satsified.% ierr = 2 failure in the line search. The iteration% is terminated if too many steplength reductions% are taken.%%% internal parameter:% debug = turns on/off iteration statistics display as% the iteration progresses%% alpha = , parameter to measure sufficient decrease%% maxarm = 10, maximum number of steplength reductions before % failure is reported%% set the debug parameter, 1 turns display on, otherwise off%debug=1;%% initialize it_hist, ierr, and set the iteration parameters%ierr = 0; maxit=40; maxdim=39;it_histx=zeros(maxit,3);maxarm=10;%if nargin == 4maxit=parms(1); maxdim=parms(2)-1;endrtol=tol(2); atol=tol(1); n = length(x); fnrm=1; itc=0; nbroy=0;%% evaluate f at the initial iterate% compute the stop tolerance%f0=feval(f,x);fc=f0;fnrm=norm(f0)/sqrt(n);it_hist(itc+1)=fnrm;it_histx(itc+1,1)=fnrm; it_histx(itc+1,2)=0; it_histx(itc+1,3)=0;fnrmo=1;stop_tol=atol + rtol*fnrm;outstat(itc+1, :) = [itc fnrm 0 0];%% terminate on entry%if fnrm < stop_tolsol=x;returnend%% initialize the iteration history storage matrices%stp=zeros(n,maxdim);stp_nrm=zeros(maxdim,1);lam_rec=ones(maxdim,1);%% Set the initial step to -F, compute the step norm%lambda=1;stp(:,1) = -fc;stp_nrm(1)=stp(:,1)'*stp(:,1);%% main iteration loop%while(itc < maxit)%nbroy=nbroy+1;%% keep track of successive residual norms and% the iteration counter (itc)%fnrmo=fnrm; itc=itc+1;%% compute the new point, test for termination before% adding to iteration history%xold=x; lambda=1; iarm=0; lrat=.5; alpha=;x = x + stp(:,nbroy);fc=feval(f,x);fnrm=norm(fc)/sqrt(n);ff0=fnrmo*fnrmo; ffc=fnrm*fnrm; lamc=lambda;%%% Line search, we assume that the Broyden direction is an% ineact Newton direction. If the line search fails to% find sufficient decrease after maxarm steplength reductions % brsola returns with failure.%% Three-point parabolic line search%while fnrm >= (1 - lambda*alpha)*fnrmo && iarm < maxarm % lambda=lambda*lrat;if iarm==0lambda=lambda*lrat;elselambda=parab3p(lamc, lamm, ff0, ffc, ffm);endlamm=lamc; ffm=ffc; lamc=lambda;x = xold + lambda*stp(:,nbroy);fc=feval(f,x);fnrm=norm(fc)/sqrt(n);ffc=fnrm*fnrm;iarm=iarm+1;end%% set error flag and return on failure of the line search%if iarm == maxarmdisp('Line search failure in brsola ')ierr=2;it_hist=it_histx(1:itc+1,:);sol=xold;return;end%% How many function evaluations did this iteration require %it_histx(itc+1,1)=fnrm;it_histx(itc+1,2)=it_histx(itc,2)+iarm+1;if(itc == 1) it_histx(itc+1,2) = it_histx(itc+1,2)+1; end;it_histx(itc+1,3)=iarm;%% terminate%if fnrm < stop_tolsol=x;rat=fnrm/fnrmo;outstat(itc+1, :) = [itc fnrm iarm rat];it_hist=it_histx(1:itc+1,:);% it_hist(itc+1)=fnrm;if debug==1disp(outstat(itc+1,:))endreturnend%%% modify the step and step norm if needed to reflect the line % search%lam_rec(nbroy)=lambda;if lambda ~= 1stp(:,nbroy)=lambda*stp(:,nbroy);stp_nrm(nbroy)=lambda*lambda*stp_nrm(nbroy);end%%% it_hist(itc+1)=fnrm;rat=fnrm/fnrmo;outstat(itc+1, :) = [itc fnrm iarm rat];if debug==1disp(outstat(itc+1,:))end%%% if there's room, compute the next search direction and step norm and % add to the iteration history%if nbroy < maxdim+1z=-fc;if nbroy > 1for kbr = 1:nbroy-1ztmp=stp(:,kbr+1)/lam_rec(kbr+1);ztmp=ztmp+(1 - 1/lam_rec(kbr))*stp(:,kbr);ztmp=ztmp*lam_rec(kbr);z=z+ztmp*((stp(:,kbr)'*z)/stp_nrm(kbr));endend%% store the new search direction and its norm%a2=-lam_rec(nbroy)/stp_nrm(nbroy);a1=1 - lam_rec(nbroy);zz=stp(:,nbroy)'*z;a3=a1*zz/stp_nrm(nbroy);a4=1+a2*zz;stp(:,nbroy+1)=(z-a3*stp(:,nbroy))/a4;stp_nrm(nbroy+1)=stp(:,nbroy+1)'*stp(:,nbroy+1);%%%else%% out of room, time to restart%stp(:,1)=-fc;stp_nrm(1)=stp(:,1)'*stp(:,1);nbroy=0;%%%end%% end whileend%% We're not supposed to be here, we've taken the maximum% number of iterations and not terminated.%sol=x;it_hist=it_histx(1:itc+1,:);ierr=1;if debug==1disp(' outstat')endfunction lambdap = parab3p(lambdac, lambdam, ff0, ffc, ffm)% Apply three-point safeguarded parabolic model for a line search. %% This code comes with no guarantee or warranty of any kind.%% function lambdap = parab3p(lambdac, lambdam, ff0, ffc, ffm)%% input:% lambdac = current steplength% lambdam = previous steplength% ff0 = value of \| F(x_c) \|^2% ffc = value of \| F(x_c + \lambdac d) \|^2% ffm = value of \| F(x_c + \lambdam d) \|^2%% output:% lambdap = new value of lambda given parabolic model%% internal parameters:% sigma0 = .1, sigma1=.5, safeguarding bounds for the linesearch %%% set internal parameters%sigma0=.1; sigma1=.5;%% compute coefficients of interpolation polynomial%% p(lambda) = ff0 + (c1 lambda + c2 lambda^2)/d1%% d1 = (lambdac - lambdam)*lambdac*lambdam < 0% so if c2 > 0 we have negative curvature and default to% lambdap = sigam1 * lambda%c2 = lambdam*(ffc-ff0)-lambdac*(ffm-ff0);if c2 >= 0lambdap = sigma1*lambdac; returnendc1=lambdac*lambdac*(ffm-ff0)-lambdam*lambdam*(ffc-ff0);lambdap=-c1*.5/c2;if (lambdap < sigma0*lambdac) lambdap=sigma0*lambdac; endif (lambdap > sigma1*lambdac) lambdap=sigma1*lambdac; end2.应用举例把以下代码复制在command 窗口中x=[1 2 3]’;f=@(x)[3*x(1)-cos(x(2)*x(3))-1/2;x(1)^2-81*(x(2)+^2+sin(x(3))+;exp(-x(1)*x(2))+20*x(3)+(10*pi-3)/3;];tol=[3,-5];[sol, it_hist, ierr] = brsola(x,f,tol)说明：以上应用举例只是给出了上文中代码的一个应用实例，具体能否得到方程的满意数值解还需要进一步调节初始给的x和tol的值。

实验2 利用matlab 解（非）线性、微分方程（组）一、实验目的1、线性方程组的解法：直接求解法和迭代法；2、非线性方程以及非线性方程组的求解；3、微分方程的数值解。

二、实验内容1、对于下列线性方程组：⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡66136221143092x （1） 请用直接法求解；>> A=[2 9 0;3 4 11;2 2 6];>> b=[13 6 6]';>> x=A\bx =7.4000-0.2000-1.4000（2） 请用LU 分解方法求解；>> [L,U]=lu(A);>> x1=U\(L\b)x1 =7.4000-0.2000-1.4000>> [L,U,P]=lu(A);>> x2=U\(L\P*b)x2 =7.4000-0.2000-1.4000（3） 请用QR 分解方法求解；>> [Q,R]=qr(A);>> x1=R\(Q\b)x1 =7.4000-0.2000-1.4000>> [Q,R,E]=qr(A);>> x2=E*(R\(Q\b))x2 =7.4000-0.2000-1.4000（4） 请用Cholesky 分解方法求解。

>> R=chol(A)Error using cholMatrix must be positive definite.因此系数矩阵A 不是正定的，故不能用Cholesky 分解法2、设迭代精度为10-6，分别用Jacobi 迭代法、Gauss-Serdel 迭代法求解下列线性方程组，并比较此两种迭代法的收敛速度。

⎪⎩⎪⎨⎧=+-=-+-=-510272109103232121x x x x x x x 解，Ax=b,新建函数如下>> A=[10 -1 0;-1 10 -2;0 -2 10];>> b=[9 7 5]';>> eps=10e-6;>> [x,n]=jacobi(A,b,[0,0,0]',eps)x =0.99370.93680.6874n =9>> [x,n]=gauseidel(A,b,[0,0,0]',eps)x =0.99370.93680.6874n =6故本例中Gauss-Serdel 迭代法优于Jacobi 迭代法3、求解非线性方程010=-+x xe x 在2附近的根。

matlab实验一：非线性方程求解-牛顿法实验一:非线性方程求解程序1：二分法：syms f x;f=input（'请输入f（x)='）;A=input（’请输入根的估计范围[a,b］=’）;e=input（'请输入根的误差限e=’);while (A（2)—A（1)）>ec=(A（1）+A(2)）/2；x=A（1）；f1=eval（f)；x=c；f2=eval（f）；if (f1＊f2）〉0A(1)=c；elseA(2)=c;endendc=（A(1)+A（2）)/2；fprintf(’c=％。

6f\na=%.6f\nb=%.6f\n'，c，A)用二分法计算方程：1．请输入f（x）=sin（x)—x^2/2请输入根的估计范围[a，b］=[1,2］请输入根的误差限e=0。

5e-005c=1.404413a=1。

404411b=1.4044152．请输入f(x)=x^3—x-1请输入根的估计范围［a，b］=[1，1.5]请输入根的误差限e=0.5e-005c=1。

324717a=1。

324715b=1。

324718程序2：newton法：syms f x；f=input（’请输入f（x)=');df=diff(f);x0=input('请输入迭代初值x0=’)；e1=input('请输入奇异判断e1=’）;e2=input（'请输入根的误差限e2='）；N=input（'请输入迭代次数限N=’)；k=1；while （k〈N）x=x0；if abs（eval（f）)〈e1 fprintf(’奇异！\nx=％。

6f\n迭代次数为:％d\n’,x0，k）breakelsex1=x0—eval(f)/eval(df）;if abs（x1-x0)<e2fprintf（'x=%。

6f\n迭代次数为:％d\n’，x1,k）breakelsex0=x1；k=k+1;endendendif k〉=Nfprintf（'失败\n'）end用newton法计算方程:1．请输入f(x）=x*exp（x)—1请输入迭代初值x0=0。

3-7变量非线性方程组的逆Broyden解法matlab程序function n_broydenclear all %清内存clc %清屏format longi=input('请输入非线性方程组个数(3-7)i=');switch icase 3syms x y zx0=input('请输入初值（三维列向量[x;y;z])=');eps=input('请输入允许的误差精度eps=');f1=input('请输入f1(x,y,z)=');f2=input('请输入f2(x,y,z)=');f3=input('请输入f3(x,y,z)=');F=[f1;f2;f3];J=jacobian(F,[x,y,z]); %求jacobi矩阵Fx0=subs(F,{x,y,z},x0);Jx0=subs(J,{x,y,z},x0);H=inv(Jx0);x1=x0-H*Fx0;k=0;while norm(x1-x0)>eps %用2范数来做循环条件p=x1-x0;q=subs(F,{x,y,z},x1)-subs(F,{x,y,z},x0);H=H-((H*q-p)*p'*H)*(p'*H*q)^-1;x0=x1;Fx0=subs(F,{x,y,z},x0);x1=x1-H*Fx0;k=k+1;endx1kcase 4syms a b c dx0=input('请输入初值（列向量[a;b;c;d])=');eps=input('请输入允许的误差精度eps=');f1=input('请输入f1(a,b,c,d)=');f2=input('请输入f2(a,b,c,d)=');f3=input('请输入f3(a,b,c,d)=');f4=input('请输入f4(a,b,c,d)=');F=[f1;f2;f3;f4];J=jacobian(F,[a,b,c,d]); %求jacobi矩阵Fx0=subs(F,{a,b,c,d},x0);Jx0=subs(J,{a,b,c,d},x0);H=inv(Jx0);x1=x0-H*Fx0;k=0;while norm(x1-x0)>eps %用2范数来做循环条件 p=x1-x0;q=subs(F,{a,b,c,d},x1)-subs(F,{a,b,c,d},x0);H=H-((H*q-p)*p'*H)*(p'*H*q)^-1;x0=x1;Fx0=subs(F,{a,b,c,d},x0);x1=x1-H*Fx0;k=k+1;endx1kcase 5syms a b c d ex0=input('请输入初值（列向量[a;b;c;d;e])='); eps=input('请输入允许的误差精度eps=');f1=input('请输入f1(a,b,c,d,e)=');f2=input('请输入f2(a,b,c,d,e)=');f3=input('请输入f3(a,b,c,d,e)=');f4=input('请输入f4(a,b,c,d,e)=');f5=input('请输入f5(a,b,c,d,e)=');F=[f1;f2;f3;f4;f5];J=jacobian(F,[a,b,c,d,e]); %求jacobi矩阵Fx0=subs(F,{a,b,c,d,e},x0);Jx0=subs(J,{a,b,c,d,e},x0);H=inv(Jx0);x1=x0-H*Fx0;while norm(x1-x0)>eps %用2范数来做循环条件 p=x1-x0;q=subs(F,{a,b,c,d,e},x1)-subs(F,{a,b,c,d,e},x0); H=H-((H*q-p)*p'*H)*(p'*H*q)^-1;x0=x1;Fx0=subs(F,{a,b,c,d,e},x0);x1=x1-H*Fx0;k=k+1;endx1kcase 6syms a b c d e fx0=input('请输入初值（列向量[a;b;c;d;e;f])='); eps=input('请输入允许的误差精度eps=');f1=input('请输入f1(a,b,c,d,e,f)=');f2=input('请输入f2(a,b,c,d,e,f)=');f3=input('请输入f3(a,b,c,d,e,f)=');f4=input('请输入f4(a,b,c,d,e,f)=');f5=input('请输入f5(a,b,c,d,e,f)=');f6=input('请输入f6(a,b,c,d,e,f)=');F=[f1;f2;f3;f4;f5;f6];J=jacobian(F,[a,b,c,d,e,f]); %求jacobi矩阵Fx0=subs(F,{a,b,c,d,e,f},x0);Jx0=subs(J,{a,b,c,d,e,f},x0);H=inv(Jx0);x1=x0-H*Fx0;k=0;while norm(x1-x0)>eps %用2范数来做循环条件 p=x1-x0;q=subs(F,{a,b,c,d,e,f},x1)-subs(F,{a,b,c,d,e,f},x0); H=H-((H*q-p)*p'*H)*(p'*H*q)^-1;x0=x1;Fx0=subs(F,{a,b,c,d,e,f},x0);x1=x1-H*Fx0;k=k+1;endkcase 7syms a b c d e f gformat longx0=input('请输入初值（列向量[a;b;c;d;e;f;g])=');eps=input('请输入允许的误差精度eps=');f1=input('请输入f1(a,b,c,d,e,f,g)=');f2=input('请输入f2(a,b,c,d,e,f,g)=');f3=input('请输入f3(a,b,c,d,e,f,g)=');f4=input('请输入f4(a,b,c,d,e,f,g)=');f5=input('请输入f5(a,b,c,d,e,f,g)=');f6=input('请输入f6(a,b,c,d,e,f,g)=');f7=input('请输入f7(a,b,c,d,e,f,g)=');F=[f1;f2;f3;f4;f5;f6;f7];J=jacobian(F,[a,b,c,d,e,f,g]); %求jacobi矩阵Fx0=subs(F,{a,b,c,d,e,f,g},x0);Jx0=subs(J,{a,b,c,d,e,f,g},x0);H=inv(Jx0);x1=x0-H*Fx0;k=0;while norm(x1-x0)>eps %用2范数来做循环条件p=x1-x0;q=subs(F,{a,b,c,d,e,f,g},x1)-subs(F,{a,b,c,d,e,f,g},x0); H=H-((H*q-p)*p'*H)*(p'*H*q)^-1;x0=x1;Fx0=subs(F,{a,b,c,d,e,f,g},x0);x1=x1-H*Fx0;k=k+1;endx1kendend。

用Broyden 方法求解非线性方程沈欢北京大学工学院，北京1008712011年10月26日摘要用Broyden 方法求解给定的非线性方程组。

1问题描述用Broyden 方法计算非线性方程组(x 1+3)(x 32−7)+18=0(1)sin (x 2e x 1−1)=0(2)的解，取初始猜测为−→x 0=(−0.5,1.4)T 。

2问题分析通过观察不难看出(0,1)点是方程的一个精确解。

取初始猜测−→x 0=(−0.5,1.4)T 在精确解附近，可以用Broyden 方法求解。

设：−−−−−−→F (x 1,x 2)=F 1(x 1,x 2)−→i +F 2(x 1,x 2)−→j 其中F 1(x 1,x 2)=(x 1+3)(x 32−7)+18;F 2(x 1,x 2)=sin (x 2ex 1−1)。

3Broyden 算法Broyden 算法类似于牛顿方法，它通过在已知的雅可比矩阵上加上一个修正项得到新的雅可比矩阵，从而避免了牛顿方法计算雅可比矩阵的复杂度。

Broyden 算法如图一所示。

Broyden 算法说明：a )初始迭代点取−→x 0=(−0.5,1.4)T ;初始的雅可比矩阵的近似矩阵取为真实的雅可比矩阵，即：A =∇F (−→x )= x 32−73(x 1+3)x 22cos (x 2e x 1−1)∗x 2∗e x 1cos (x 2e x 1−1)∗e x 1(3)1图1:Broyden方法的算法简图A0=∇F(−→x)=−4.256014.70000.83950.5996(4)取精度为10−8。

并计算得到−→x1=[−0.0553,1.0281]T。

b)在每步迭代中用−→x k和−−→x k−1来估计新的雅可比近似矩阵。

A k=A k−1+g k−A k−1.y kyk.y k∗y T k(5)其中：g k=F(−−→x k−1)−F(−→x k),y k=−−→x k−1−−→x k。

第7章 求解非线性方程7.1 多项式运算在MATLAB 中的实现一、多项式的表达n 次多项式表达为：n a +⋯⋯++=x a x a x a p(x)1-n 1-n 1n 0，是n+1项之和 在MATLAB 中，n 次多项式可以用n 次多项式系数构成的长度为n+1的行向量表示[a0, a1,……an-1,an]二、多项式的加减运算 设有两个多项式na +⋯⋯++=x a x a x a p1(x)1-n 1-n 1n 0和m b +⋯⋯++=x b x b x b p2(x)1-m 1-m 1m 0。

它们的加减运算实际上就是它们的对应系数的加减运算。

当它们的次数相同时，可以直接对多项式的系数向量进行加减运算。

当它们的次数不同时，应该把次数低的多项式无高次项部分用0系数表示。

例2 计算()()1635223-+++-x x x xa=[1, -2, 5, 3]; b=[0, 0, 6, -1]; c=a+b例 3 设()6572532345++-+-=x x x x x x f ，()3532-+=x x x g ，求f(x)+g(x)f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; g1=[0, 0, 0, g];%为了和f 的次数找齐f+g1, f-g1三、多项式的乘法运算conv(p1,p2)例4 在上例中，求f(x)*g(x)f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; conv(f, g)四、多项式的除法运算[Q, r]=deconv(p1, p2)表示p1除以p2，给出商式Q(x)，余式r(x)。

Q,和r 仍为多项式系数向量 例4 在上例中，求f(x)/g(x)f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; [Q, r]=deconv(f, g) 五、多项式的导函数p=polyder(P)：求多项式P 的导函数 p=polyder(P,Q)：求P ·Q 的导函数[p,q]=polyder(P,Q)：求P/Q 的导函数，导函数的分子存入p ，分母存入q 。

非线性方程的解法(含牛拉解法)1引 言数学物理中的许多问题归结为解函数方程的问题，即，0)(=x f (1.1) 这里，)(x f 可以是代数多项式，也可以是超越函数。

若有数*x 为方程0)(=x f 的根，或称函数)(x f 的零点。

设函数)(x f 在],[b a 内连续，且0)()(<b f a f 。

根据连续函数的性质知道，方程0)(=x f 在区间],[b a 内至少有一个实根；我们又知道，方程0)(=x f 的根，除了极少简单方程的根可以用解析式表达外，一般方程的根很难用一个式子表达。

即使能表示成解析式的，往往也很复杂，不便计算。

所以，具体求根时，一般先寻求根的某一个初始近似值，然后再将初始近似值逐步加工成满足精度要求为止。

如何寻求根的初始值呢？简单述之，为了明确起见，不妨设)(x f 在区间],[b a 内有一个实的单根，且0)(,0)(><b f a f 。

我们从左端出点a x =0出发，按某一预定的步长h 一步一步地向右跨，每跨一步进行一次根的“搜索”，即检查每一步的起点k x 和1+k x (即，h x k +)的函数值是否同号。

若有：0)(*)(≤+h x f x f k k (1.2) 那么所求的根必在),(h x x k k +内，这时可取k x 或h x k +作为根的初始近似值。

这种方法通常称为“定步长搜索法”。

另外，还是图解法、近似方程法和解析法。

2 迭代法2.1 迭代法的一般概念迭代法是数值计算中一类典型方法，不仅用于方程求根，而且用于方程组求解，矩阵求特征值等方面。

迭代法的基本思想是一种逐次逼近的方法。

首先取一个精糙的近似值，然后用同一个递推公式，反复校正这个初值，直到满足预先给定的精度要求为止。

对于迭代法，一般需要讨论的基本问题是：迭代法的构造、迭代序列的收敛性天收敛速度以及误差估计。

这里，主要看看解方程迭代式的构造。

对方程(1.1)，在区间],[b a 内，可改写成为：)(x x ϕ= (2.1) 取],[0b a x ∈，用递推公式：)(1k k x x ϕ=+, ,2,1,0=k (2.2) 可得到序列：∞==0210}{,,,,k k k x x x x x (2.3)当∞→k 时，序列∞=0}{k k x 有极限x ~，且)(x ϕ在x ~附近连续，则在式(2.2)两边极限，得， )~(~x x ϕ= 即，x ~为方程(2.1)的根。

一种改进的逆Broyden算法景书杰;赵建卫;韩学锋【摘要】对解非线性方程组Broyden方法和逆Broyden方法进行了改进,构造了求解非线性方程组F(x)=0的一个迭代公式,并讨论了其收敛性,说明该算法是有效的.【期刊名称】《吉林师范大学学报（自然科学版）》【年(卷),期】2015(036)001【总页数】4页(P66-68,71)【关键词】非线性方程组;Broyden方法;逆Broyden方法;收敛性【作者】景书杰;赵建卫;韩学锋【作者单位】河南理工大学数学与信息科学学院,河南焦作454003;河南理工大学数学与信息科学学院,河南焦作454003;河南理工大学数学与信息科学学院,河南焦作454003【正文语种】中文【中图分类】O2210 引言设映像F：D⊂Rn→Rn，求解非线性方程组F(x)=0(1)的Broyden法是对拟Newton法中的Ak进行修正的计算方案，它每步迭代只需要计算n个分量函数值及O(n2)次算术运算，大大降低了Newton法每步的计算量.它的具体形式如下：(2)其中，sk=xk+1-xk,yk=F(xk+1)-F(xk).王斌[1]对逆Broyden秩1拟Newton方法给出了新的推导思路；司智勇[2]、陈新一[3]、刘慧[4]等对Newton型迭代从不同的方面进行了修正，并得到了理想的结果；本文正是受[2]的启发对逆Broyden算法进行了新的修正，借鉴陈兰平[5-6]、赵云彬[7]、谢铁军[8]等对收敛性证明的思路，参考相关资料[10-12]知识，得到了本篇文章.先假定已求得的方程(1)的k次近似解xk∈D，记xk=xk,0，再取另外t个辅助点xk,1，…，xk,t，我们用xk,t来作为xk+1的值，则Broyden算法变为如下新算法：(3)其中，sk=xk+1-xk,yk=F(xk+1)-F(xk)引理0：当矩阵A∈L(Rn)非奇异，u,ν∈Rn，则A+uνT非奇异，且当σ=1+νTA-1u≠0时，有这个引理的证明是显然的.由引理0可知：存在，令则可由下式给出其中(sk)TBkyk≠0.那么(3)就可改写成逆Broyden新算法：(4)其中xk,0=xk.1 算法步1 给出初始近似x0及精度要求ε1及ε2；步2 计算初始矩阵B0(一般B0=[F′(x0)]-1或者B0=I，I为单位矩阵)；步3 k=0；计算F(x0)；步4 xk+1=xk,t，计算sk=-BkF(xk)，及xk+1=xk+sk；步5 计算F(xk+1)；检验‖F(xk+1)‖≤ε2或者‖sk‖≤ε1，若成立则转步8，否则做步6；步6 计算yk=F(xk+1)-F(xk)，并由(4)计算Bk+1；步7 k+1→k，F(xk+1)→F(xk)，Bk+1→Bk，xk+1→xk转步4；步8 x*=xk+1，打印x*，‖F(xk+1)‖，‖sk‖后结束.2 收敛性分析下面讨论新算法的收敛性，首先讨论算法(3)的收敛性，这里首先给出如下引理：引理1 设F：D⊂Rn→Rn在开凸集D0⊂D上连续可导，且存在x*⊂D0，使得F(x*)=0，F′(x*)非奇异，{Ak}⊂L(Rn)为非奇异矩阵列，若对某个x0∈D0，由拟Newton法生成的序列xk⊂D0收敛于x*，则此序列超线性收敛的充要条件是(*)引理2 设F：D⊂Rn→Rn满足引理1的条件，又在x*处存在常数L>0及η>0使得‖F′(x)-F′(x*)‖≤L‖x-x*‖，∀x∈D0，(5)‖Ak-F′(x*)‖≤η‖xk-x*‖，k=0，1，…，则由拟Newton法生成的序列{xk}至少平方收敛于x*.引理3 若x*∈D0是F(x)=0的解.F：D⊂Rn→Rn在开凸集D0⊂D上连续可导，且F′(x*)非奇异，条件(5)成立.U(x,A)是定义在D×DM上F的修正函数，若存在α1>0及α2>0使(6)其中则对每个γ∈(0，1)，∃ε(γ)>0，及δ(γ)>0，使当x0∈D0，A0∈DM满足‖x0-x*‖≤ε(γ),‖A0-F′(x*)‖≤δ(γ)时，逆Newton法生成的序列{xk}是适定的且收敛于x*，并有‖xk+1-x*‖2≤γ‖xk-x*‖2,k=0，1，…序列{‖Ak‖}及一致有界.引理1，引理2，引理3的证明可以参考[9]．定理设F：D⊂Rn→Rn在x*凸领域D0⊂D内连续可导，F(x*)=0，且F′(x*)非奇异，F′(x)满足条件(5)，则改进的算法(3)产生的序列{xk}超线性收敛到x*.证明：先证明(2)生成的序列{xk}1收敛于x*，这只需要验证条件(6)成立即可，为此我们取DM=Qyt，并注意上式两端取模，得(7)由中值定理有又于是由(7)得这表明条件(6)当α1=0及α2=L时成立，于是引理3条件全部满足，故(3)生成的序列{xk}1收敛于x*.值得注意的是{xk}和{xk}1的区别是：前者每项中的每项都是由后者中每个对应项的t个辅助点来近似的.因而由序列{xk}1收敛于x*，那么必有{xk}收敛于x*.下面进一步证明{xk}是超线性收敛的，根据引理1，只要证明(*)成立即可.为此利用直接计算得∀E∈L(Rn),现在利用不等式(α2-β2)1/2≤α-(2α)-1β2，(α≥β≥0)即可推得若令E=Ak-F′(x*)而ηk=‖E‖F，将这些结果用到不等式(7)则得σk=max{‖xk+1-x*‖,‖xk-x*‖}．由于ηk+1≤ηk+Lσk且{xk}线性收敛，故{ηk}有界，若记ηk≤η(对一切k)，则有于是又为单调增序列，故收敛，于是φk=0，即由引理1可得序列{xk}超线性收敛.从前面的推导可以看出，在定理0的条件下(3)和(4)相通的，我们的定理已经证明了算法(3)是超线性收敛的，同理可证算法(4)也是超线性收敛的.参考文献【相关文献】[1]王斌.非线性方程组的逆Broyden秩1拟Newton方法及其在MATLAB中的实现[J].云南大学学报(自然科学版)，2008,30(S2):144～148.[2]司智勇，阿不都热西提·阿不都外力，张知难.Newton型方法的推广与改进[J].数值计算与计算机应用,2008,29(3):171～175.[3]陈新一.Newton迭代法的一个改进[J].数学的实践与认识，2006,36(2):291～294.[4]刘慧，张立杰.解非线性方程组的一种新的Newton型迭代法[J].新疆工学院学报，1997,18(3):171～175.[5]陈兰平，焦宝聪.一类非拟Newton算法及其收敛性[J].应用数学与计算数学学报，1997,11(2):9～17．[6]陈兰平，王丽伟.广义拟牛顿算法对一般目标函数的收敛性[J].应用数学,2002,15(3)：69～75．[7]赵云彬，段虞荣.伪Newton-δ族算法对一般目标函数的收敛性[J].数值计算与计算机应用,1996,17(1):36～47.[8]谢铁军，陈明文，刘任平.无记忆拟牛顿方法的收敛性[J].运筹与管理,2000，9(4)：57～61．[9]李庆扬，莫孜中，祁力群.非线性方程组的数值解法[M].北京:科学出版社，1987.[10]J.M.Ortega，W.C.Rheinboldt.多元非线性方程组的迭代解法[M].北京:科学出版社，1987.[11]赵东方.数学模型与计算[M].北京：科学出版社，2007.[12]杨晓光.一类拟牛顿算法的收敛性[C].清华大学应用数学论文集,1991，113～117．。

数值分析中求解非线性方程的MATLAB求解程序（6种）1.求解不动点function [k,p,err,P]=fixpt(g,p0,tol,max1)%求解方程x=g(x) 的近似值，初始值为p0%迭代式为Pn+1=g(Pn)%迭代条件为：在迭代范围内满足|k|<1(根及附近且包含初值)k为斜率P(1)=p0;for k=2:max1P(k)=feval(g,P(k-1));err=abs(P(k)-P(k-1));relerr=err/(abs(P(k))+eps);p=P(k);if (err<tol)|(relerr<tol)break;endendif k==max1disp('超过了最长的迭代次数')endP=P';2.二分法function [c,err,yc]=bisect(f,a,b,delta)%二分法求解非线性方程ya=feval(f,a);yb=feval(f,b);if ya*yb>0break;endmax1=1+round((log(b-a)-log(delta))/log(2));for k=1:max1c=(a+b)/2;yc=feval(f,c);if yc==0a=c;b=c;elseif yb*yc>0b=c;yb=yc;elsea=c;ya=yc;endif b-a<deltabreak;endendc=(a+b)/2;err=abs(b-a);yc=feval(f,c);3.试值法function [c,err,yc]=regula(f,a,b,delta,epsilon,max1) %试值法求解非线性方程%f(a)和飞(b)异号ya=feval(f,a);yb=feval(f,b);if ya*yb>0disp('Note:f(a)*f(b)>0');endfor k=1:max1dx=yb*(b-a)/(yb-ya);c=b-dx;ac=c-a;yc=feval(f,c);if yc==0break;elseif yb*yc>0b=c;yb=yc;elsea=c;ya=yc;enddx=min(abs(dx),ac);if abs(dx)<delta|abs(yc)<epsilonbreak;endendc;err=abs(b-a)/2;yc=feval(f,c);4.求解非线性方程根的近似位置function R=approot(X,epsilon)%求解根近似位置%为了粗估算方程f(x)=0在区间[a,b]的根的位置，%使用等间隔采样点(xk,f(xk))和如下的评定准则：%f(xk-1)与f(xk)符号相反，%或者|f(xk)|足够小且曲线y=f(x)的斜率在%(xk,f(xk))附近改变符号。

Matlab实例源码教程：如何用MATLAB求解非线性微分方程Matlab实例源码教程:如何用MATLAB求解非线性微分方程做一个最基本的假设:你们都看过高数。

一。

老湿发话了:童鞋们，求解一下这个方程，判断她是否稳定。

要是稳定，那么她是否存在极限环:一看明白了，这不就是传说中的范德普方程。

地球人都知道她稳定并有极限环。

现在我们就看看如何用MATLAB求解她的轨迹。

二。

一般的计算机求解方程的方法无外乎是这样:首先把该方程改写成一个规范的形式，一般使用状态空间表示法;而后调用已有的算法进行求解;最后对得出的结果进行处理，比如画图之类的。

接下来就对这三大步分别作出解释。

三。

输入待求解的方程。

首先我们知道，状态空间的标准形式(自由系统)是:这里X是列向量，F是作用于列向量的函数，可以是线性也可以是非线性。

范德普方程可以改写成这样的标准形式:MATLAB中关于输入输入待解方程的语句特别简单。

需要先定义一个普通函数，函数名任意，姑且叫做myFcn，格式如下 function xdot = myFcn (t, x) 需要注意的是，函数必须含有t, x两个参数，名称可以自己任意定。

xdot是这个函数本身的返回值，只出现在这个函数内部，因此也可以任意定。

定义中的x被视为是一个列向量，x(i)表示列向量中的第i个分量。

那么F函数的每一个分量即简单地用表达时给出即可。

其中的自变量可以引用x(i)。

以范德普方程为例:xdot = [x(2) ; u(1-x(1)^2)*x(2)-x(1)]于是，这两句话便构成了待解函数。

四。

调用MATLAB函数进行求解通常人工求解微分方程需要知道初始值，计算机求解也不例外。

另外，由于非线性方程一般只有数值解，故计算精度也可以调整。

这些都是可以自己调整的参数。

调用MATLAB计算求解常微分方程的模式很简单，格式为:[t, x] = ode45(@myFcn, [开始时间结束时间], [初始值列向量], options) 注意到求解出来的结果是[t, x]即一堆数，所以需要我们进行后处理比如画图之类的。

线性方程组数值方法(1) % jacobi 迭代法计算线性方程组function [x,k]=Fjacobi(A,b,x0,tol)% tol 为输入误差容限，x0为迭代初值max1= 300; %默认最多迭代300，超过要300 次给出警告D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=D\(L+U);f=D\b;x=B*x0+f;k=1; %迭代次数while norm(x-x0)>=tolx0=x;x=B*x0+f;k=k+1;if(k>=max1)disp('迭代超过300 次，方程组可能不收敛');return;end%[k x'] %如果想显示每一步迭代的结果可以把这行命令保留end(2) %高斯-塞德尔迭代法计算线性方程组function [x,k]=Fgseid(A,b,x0,tol)% tol 为误差容限max1= 300; %默认最高迭代300 次D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);G=(D-L)\U;f=(D-L)\b;x=G*x0+f;k=1;while norm(x-x0)>=tolx0=x;x=G*x0+f;k=k+1;if(k>=max1)disp('迭代次数太多，可能不收敛');return;end% [k,x'] %如果要显示每一步迭代结果，可以这行命令。

end(3) %超松弛迭代法计算线性function [x,k]=Fsor(A,b,x0,w,tol)%08 年/4 月/14 日%tol 需要的计算精度max = 300; %迭代最大次数%对松弛因子做出限制if(w<=0 || w>=2)error;return;endD=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵B=inv(D-L*w)*((1-w)*D+w*U);f=w*inv((D-L*w))*b;x=B*x0+f;k=1; %迭代次数while norm(x-x0)>=tolx0=x;x =B*x0+f;k=k+1;if(k>=max)disp('迭代次数太多，SOR方法可能不收敛');return;end% [k,x'] %如果要显示中间结果，可以保留这行语句end(4) % Gauss 消去法计算线性方程组function X=Fgauss(A,b)zengguang=[A b]; n=length(b);ra=rank(A);rz=rank(zengguang);temp1=rz-ra;if temp1>0,disp('无一般意义下的解，系数矩阵与增广矩阵的秩不同') returnendif ra==rzif ra==nX=zeros(n,1);C=zeros(1,n+1);for p= 1:n-1for k=p+1:nm= zengguang(k,p)/ zengguang(p,p);zengguang(k,p:n+1)= zengguang(k,p:n+1)-m* zengguang(p,p:n+1); endendb=zengguang(1:n,n+1);A=zengguang(1:n,1:n);X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('方程为欠定方程。