3-7变量非线性方程组的逆Broyden解法matlab程序
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解(非)线性、微分方程(组)
实验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实验一：非线性方程求解-牛顿法
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程序
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方法解非线性方程组
用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。
非线性方程组求解及matlab实现讲解
x0
X
例:牛顿法计算x^2-25=0的解
f(x)=x2-25，则f’(x)=2x 可构造迭代公式如下：
xi2 25 xi 1 xi 2 xi
取x0=2代入上式，得x1=7.25，继续递推， 依次得5.35、5.0114、5.000001、5.0000000001 …
牛顿法注意事项
逐步扫描法计算示例-方程x2-2=0的正数解
计算方程 x 2 2 0 的正数解
二分法
若函数f(x)在区间[a,b]内单调连续，且f(a)f(b)<0， 则在闭区间[a,b]内必然存在方程f(x)=0的根x*
二分法的图形解释 二分法的MATLAB程序
k=0; while abs(b-a)>eps x=(a+b)/2; if sign(f(x))==sign(f(b)) b=x; else a=x; end k=k+1; end
f '( x) 0, f "( x) 连续且不变号，则只 在有根区间[a,b]上， 要选取的初始近似根x0满足 f ( x0 ) f "( x0 ) 0 ，切线法 必定收敛。 在单根附近，牛顿公式恒收敛，而且收敛速度很快。 但是需要注意如果初始值不在根的附近，牛顿公式 不一定收敛 在实际使用中，牛顿法最好与逐步扫描法结合起来， 先通过逐步扫描法求出根的近似值，然后用牛顿公 式求其精确值，以发挥牛顿法收敛速度快的优点
c x
不动点迭代法
从给定的初值x0，按上式可以得到一个数列： { x0, x1, x2, …, xk, … }
如果这个数列有极限，则迭代格式是收敛的。 * x xk 就是方程的根 这时数列{xk}的极限 lim k 上述求非线性代数方程式数值解的方法称为直 接迭代法（或称为不动点迭代法）。这个方法 虽然简单，但根本问题在于当k->∞时，xk是否 收敛于x*，也就是必须找出收敛的充分条件
用matlab求解非线性方程组的几种方法之程序
教育电子音像出版社 作者：任玉杰 第二章 非线性方程（组）的数值解法的 MATLAB 程序
第二章
非线性方程（ 非线性方程（组）的数值解法
本章主要介绍方程根的有关概念，求方程根的步骤，确定根的初始近似值的方法（作图 法，逐步搜索法等） ，求根的方法（二分法，迭代法，牛顿法，割线法，米勒（Müller）法 和迭代法的加速等）及其 MATLAB 程序，求解非线性方程组的方法及其 MATLAB 程序.
3 2
2.3 二分法及其 MATLAB 程序
2.3.1 二分法
二分法也称逐次分半法.它的基本思想是：先确定方程 f ( x ) = 0 含根的区间（a,b） ，再 把区间逐次二等分. 编写二分法求方程根的迭代次数的 MATLAB 我们可以根据式 （2.3b） 、 区间[a,b]和误差 ε ， 命令. 已 知闭 区间[a,b] 和 误差 ε .用二分法求方程 误差 不大于 ε 的根的迭代 次 数 k 的 MATLAB 命令 k=-1+ceil((log(b-a)- log(abtol))/ log(2)) % ceil 是向+ ∞ 方向取整，abtol 是误差 ε.
10．
高等教育出版社
教育电子音像出版社 作者：任玉杰 第二章 非线性方程（组）的数值解法的 MATLAB 程序
for k=2:n X(k)=a+k*h;Y(k)=funs(X(k)); %程序中调用的funs.m为函数 sk=Y(k)*Y(k-1); if sk<=0, m=m+1;r(m)=X(k); end xielv=(Y(k+1)-Y(k))*(Y(k)-Y(k-1)); if (abs(Y(k))<tol)&( xielv<=0) m=m+1;r(m)=X(k); end end 例2.2.4 用逐步搜索法的MATLAB程序分别求方程 2 x + 2 x − 3 x − 3 = 0 和 sin(cos 2x 3 ) = 0 在区间 [ −2,2] 上的根的近似值，要求精度是0.000 1. 解 在MATLAB工作窗口输入如下程序 >> [k,r]=zhubuss(-2,2,0.001,0.0001) 运行后输出的结果 k =4001 r = -1.2240 -1.0000 -1.0000 -0.9990 1.2250 即搜索点的个数为k =4 001，其中有5个是方程⑴的近似根，即r = -1.224 0，-1.000 0，-1.000 0，-0.999 0，1.225 0，其精度为0.000 1. 在程序中将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.9220 0.9230 1.3310 1.5780 1.7650 1.9200 如果读者分别将方程⑴的结果与例2.2.3比较，方程⑵的结果与例2.1.2比较，将会发现 逐步搜索法的MATLAB程序的优点.如果精度要求比较高，用这种逐步搜索法是不合算的.
MATLAB应用 求解非线性方程
第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 。
matlab非线性方程的解法(含牛拉解法)
非线性方程的解法(含牛拉解法)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方法和逆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种)
数值分析中求解非线性方程的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求解非线性方程组
１ 引言
［ｘ， ｆｖａｌ， ｅｘｉｔｆｌａｇ， ｏｕｔｐｕｔ］＝ｆｓｏｌｖｅ（…）返 回 一 个 包 含
非 线 性 方 程 组 解 的 几 何 意 义 与 线 性 方 程 最优化信息的输出结构 ｏｕｔｐｕｔ。
组 类 似， 方 程 组 中 每 个 方 程 定 义 了 一 个“曲 ”超 ［ｘ， ｆｖａｌ， ｅｘｉｔｆｌａｇ， ｏｕｔｐｕｔ， ｊａｃｏｂｉａｎ］＝ｆｓｏｌｖｅ（… ）返 回
４ 迭代方法程序
数可以建立符号变量、表达式和矩阵。Ｍａｔｌａｂ 的
一个多世纪以来， 迭代法一直被人们研
符 号 处 理 功 能 可 以 对 符 号 对 象 进 行 因 式 分 解 、 究、使用和发展。近些年来， 求解非线性方程组
替 换 、化 简 等 处 理 以 及 进 行 微 积 分 、求 极 限 、线 的迭代法越来越受到人们的重视， 并为许多计
ｆ（ｘ）及 ｘ＝ｆ（ｔ）、ｆ（ｘ）＝０ 构 成 的 参 数 曲 线 ， ｅｚｐｏｌａｒ 可 迭代初值的选取方法， 三是证明迭代方法的保
以 绘 制 ｒ＝ｆ（θ）的 极 坐 标 曲 线 ， ｅｚｐｌｏｔ３ 可 绘 制 ｙ＝ｆ 正性， 还有一些经典迭代法和外推迭代法的最
（ｔ）、ｘ＝ｆ（ｔ）、ｚ＝ｆ（ｔ）构 成 的 参 数 曲 线 ， 还 有 ｅｚｓｕｒｆ、 佳参数问题、在实际使用迭代法时如何建立可
平面， 非线性方程组的解为所有超平面的交点， 一个基于解的雅可比行列式 ｆｕｎ。
但是这些曲面可能相交， 也可能不相交， 情况比 求 解 方 程 之 前 ， 需 要 建 立 一 个 ｍ 程 序 定 义
平面复杂。通常对一个二维或三维没有解析解 “ｆｕｎ”， 即所求的非线性方程组，程序如下：
非线性方程组的逆Broyden秩1拟Newton方法及其在MATLAB中的实现
Ak +1 = Ak + ΔA k .
( 5)
其中 ΔA k 为增量矩阵且 rank (ΔA k) ≥1 . 式 (3) ～ (5) 便构成 Bro yden 拟 Newto n 法公式.
这里仅考虑 rank (ΔAk) = 1 时的方法 ,即 B ro yden 秩 1 拟 Newto n 方法.
)
T
Bk.
从而可得到与式 (13) 相应的 Broy den 秩 1 公式 :
x ( k +1) = x ( k) - Bk F ( x( k) ) ,
p ( k) = x ( k+ 1) - x ( k) ,
q ( k) = F ( x ( k +1) ) - F ( x ( k) ) ,
[ A k + u ( k) ( v ( k) ) T ] - 1 =
( Ak + ΔA k) - 1 =
A
k
1 +1.
( 14) ( 15)
由式(9 ) , (10 ) 得
A
-1 k
u
(
k)
=
A
k
1
q ( k) - A kp ( k) ( v ( k) ) T p ( k)
=
A
k
1
q
(
k)
-
Matlab实例源码教程：如何用MATLAB求解非线性微分方程
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]即一堆数，所以需要我们进行后处理比如画图之类的。
解线性与非线性方程(组)Matlab代码
线性方程组数值方法(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('方程为欠定方程。
MATLAB解方程组(线性与非线性方程组)
例7-9 求下列非线性方程组在(0.5,0.5) 附近的数值解。 (1) 建立函数文件myfun.m。 function q=myfun(p) x=p(1); y=p(2); q(1)=x-0.6*sin(x)-0.3*cos(y); q(2)=y-0.6*cos(x)+0.3*sin(y); (2) 在给定的初值x0=0.5,y0=0.5下，调用fsolve函数求方程的根。 x=fsolve('myfun',[0.5,0.5]',optimset('Display','off')) x= 0.6354 0.3734
2．Gauss-Serdel迭代法 在Jacobi迭代过程中，计算时，已经得到，不必再用，即原来的迭代
公式Dx(k+1)=(L+U)x(k)+b可以改进为Dx(k+1)=Lx(k+1)+Ux(k)+b， 于是得到：
x(k+1)=(D-L)-1Ux(k)+(D-L)-1b 该式即为Gauss-Serdel迭代公式。和Jacobi迭代相比，Gauss-Serdel
7.1.2 迭代解法 迭代解法非常适合求解大型系数矩阵的方程组。在数值分析中，迭代
解法主要包括 Jacobi迭代法、Gauss-Serdel迭代法、超松弛迭代法 和两步迭代法。
1．Jacobi迭代法 对于线性方程组Ax=b，如果A为非奇异方阵，即aii≠0(i=1,2,…,n)，则
可将A分解为A=D-L-U，其中D为对角阵，其元素为A的对角元素， L与U为A的下三角阵和上三角阵，于是Ax=b化为： x=D-1(L+U)x+D-1b 与之对应的迭代公式为：
(2) QR分解
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
3-7变量非线性方程组的逆Broyden解法
matlab程序
function n_broyden
clear all %清内存
clc %清屏
format long
i=input('请输入非线性方程组个数(3-7)i=');
switch i
case 3
syms x y z
x0=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;
end
x1
k
case 4
syms a b c d
x0=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;
end
x1
k
case 5
syms a b c d e
x0=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;
end
x1
k
case 6
syms a b c d e f
x0=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;
end
k
case 7
syms a b c d e f g
format long
x0=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;
end
x1
k
end
end。