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 。
- 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。