数值分析算法 MATLAB 实践 非线性方程(组)求解

4.1 二分法

function [c,err,yc]=bisectmethod(fun,a,b,delta)
% 二分法function [c,err,yc]= bisectmethod (fun,a,b,delta)
%a为左区间，b为右区间，delta为区间误差限 
%c求解的值，err误差，yc:c在f(x)处的值
    ya=feval('fun',a);
    yb=feval('fun',b);
    if yb==0
        c=b;
        return;
    end
    if(ya*yb)>0
        disp('(a,b)不是有根区间');
        return;
    end
    max1=1+round((log(b-a)-log(delta))/log(2));
    for k=1:max1
        c=(a+b)/2;
        yc=feval('fun',c);
        if yc==0
            a=c;
            b=c;
            return;
        elseif(yb*yc)>0
                b=c;
                yb=yc;
        else
                a=c;
                ya=yc;
        end
        if(b-a)<delta
            break;
        end
    end
    c=(a+b)/2;
    err=abs(b-a);
    yc=feval('fun',c);
end

% 二分法function [c,err,yc]=bisectmethod(fun,a,b,delta)
%a为左区间，b为右区间，delta为区间误差限 
%c求解的值，err误差，yc:c在f(x)处的值
eps = 1e-4;
[x0,err,y_x0]=bisectmethod(@fun,-2,-1,eps);
disp(['方程组的解：x0= ']);
disp(x0)

function [xc,k_cnt] = bisect(f,a,b,tol)
% 二分法function [xc,k_cnt] = bisect(f,a,b,tol)
%a为左区间，b为右区间，tol为区间误差限 
%xc求解的值
    if sign(f(a))*sign(f(b))>=0
        error('f(a)f(b)eps
    x0=r;
    Fx = subs(F,transpose(symvar(F)),x0);
    dFx = subs(dF,transpose(symvar(F)),x0);
    r=vpa(x0-inv(dFx)*Fx');
    tol=norm(r-x0)
    if(n>N)
disp('迭代步数太多，可能不收敛！');break;
    end
    allx(:,n)=x0;
    ally(:,n)=Fx;
    n=n+1;
end

function f = Multifun(x)


% 只需要将fun.m文件中的 k更改为所需维数，文件中的方程替换为你自己所需的方程
% x(1)->x x(2)->y x(3)->z
%{
    f(1)=3*x(1)-cos(x(1)*x(2))-1/2;
    f(2)=x(1)^2-81*(x(2)+0.1)^2+sin(x(3))+1.06;
    f(3)= exp(-1*x(1)*x(2))+20*x(3)+(10*pi-3)/3;
%}
     k=3;
    for i=1:k
      x(i)=sym (['x',num2str(i)]);
    end 

  f(1)=3*x(1)-cos(x(1)*x(2))-1/2;
  f(2)=x(1)^2-81*(x(2)+0.1)^2+sin(x(3))+1.06;
  f(3)= exp(-1*x(1)*x(2))+20*x(3)+(10*pi-3)/3;
end
allx0 = [1,1,1];
[allx,ally,r,k_cnt]=MultiNewton(Multifun,allx0);
disp(['迭代次数：k_cnt= ']);
disp(k_cnt)
disp(['方程组的解：allx= ']);
disp(allx)