数值分析算法 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)<0 not satisfied!') % 停止运行
    end
    fa = f(a);
    fb = f(b);
    k_cnt=0;
    while (b-a)/2 > tol
        c = (a+b)/2;
        fc = f(c);
        if fc == 0
            break
        end
  if sign(fc)*sign(fa) < 0 
% [a,c]是新的区间
           b = c;
          fb = fc;
        else     % [b,c]是新的区间
            a = c;
            fa = fc;
        end
        k_cnt=k_cnt+1;
    end
xc = (a+b)/2; 
 % 新的中点就是最优估计
end

% 二分法function [xc,k_cnt] = bisect(f,a,b,tol)
%a为左区间，b为右区间，tol为区间误差限 
%xc求解的值
[x1,k1_cnt] = bisect(@fun,-2,-1,eps);
disp(['迭代次数：k1_cnt= ']);
disp(k1_cnt)
disp(['方程组的解：x1= ']);
disp(x1)

function y=fun(x)
    y=x.^3-3*x-1;
end

4.2 牛顿迭代法

function[p1,err,k,y]=NewtonFunmethod(fun,dfun,p0,delta,max1)
%fun是给定的非线性函数  %p0为初始值  %delta为给定误差界
%max1迭代次数的上限  %p1为所求得的方程的近似解  %err为p1-p0的绝对值
%k为所需要的迭代次数  %y=fun(p1)
for k=1:max1
p1 = p0-fun(p0)/dfun(p0);
%一定要有fun,dfun .m文件
%p1=p0feval('fun',p0)/feval('dfun',p0);
 err=abs(p1-p0);
  p0=p1;
  %y=feval('fun',p1);
    y= fun(p1);
   if(err<delta)|(y==0)
       break;
     end
end
end

%function [p1,err,k,y]=NewtonFunmethod(fun,dfun,p0,delta,max1)
%NewtonFunmethod对于初值比较敏感，选择不好，出现不收敛情况
delta = 1e-7;
p0 = 1.2;
it_max = 100;
[x2,err,k2_cnt,y]=NewtonFunmethod(@fun1,@dfun1,p0,delta,it_max);
disp(['迭代次数：k2_cnt= ']);
disp(k2_cnt)
disp(['方程组的解：x2= ']);
disp(x2)

function y=fun1(x)
    y=x.^3-3*x+2;
end
function y=dfun1(x)
    y=3*x.^2-3;
end

4.3 割线法

function[p1,err,k,y]=secantFunmethod(fun,p0,p1,delta,max1)
%割线法fun是给定的非线性函数  %p0，p1为初始值  %delta为给定误差界
%max1迭代次数的上限  %p1为所求得的方程的近似解  %err为p1-p0的绝对值
%k为所需要的迭代次数  %y=fun(p1)
for  k=1:max1
    p2 = p1-fun(p1)*(p1-p0)/(fun(p1)-fun(p0));
    err=abs(p2-p1);
    p0=p1;
    p1=p2;
    y = fun(p1);
    if(err<delta)||(abs(y)<0.00001)
          break;
    end
end

delta = 1e-6;
p0 = 1.0;
p1 = 1.5;
it_max = 100;
[x3,err,k3_cnt,y]=secantFunmethod(@fun2,p0,p1,delta,it_max);
disp(['迭代次数：k3_cnt= ']);
disp(k3_cnt)
disp(['方程组的解：x3= ']);
disp(x3)

function y=fun2(x)
    y=x.^3+4*x.^2-10;
end

4.4 非线性方程组-牛顿法

%牛顿迭代法 计算非线性方程组
%输入 x0 为迭代初值
%tol 为误差容限 如果缺省 默认为 10 的-10 次方
%data 用来存放计算的中间数据便于计算收敛情况分析
function [x,n,data]=new_ton(x0,tol)
    if nargin==1
         tol=1e-10;
    end
    x1=x0-f1(x0)/df1(x0);
    n=1;
    %迭代过程
    while (norm(x1-x0)>tol)
        x0=x1;
        x1=x0-f1(x0)/df1(x0);
        n=n+1;
        %data 用来存放中间数据
        data(:,n)=x1;
    end
    x=x1;
end
%牛顿迭代法的 方程函数
function f=f1(x0)
x=x0(1);
y=x0(2);
f1=x^2-2*x-y+0.5;
f2=x^2+4*y^2-4;
%最后方程函数 以行向量输出
f=[f1 f2];
end
% 求解非线性方程表达式的雅可比矩阵
function f=df1(x0)
    x=x0(1);
    y=x0(2);
    f=[2*x-2 -1
2*x 8*y];
end

x1=[1 1];
[x_1,k1_cnt,data]=new_ton(x1);
disp(['迭代次数：k1_cnt= ']);
disp(k1_cnt)
disp(['方程组的解：x_1= ']);
disp(x_1)
%抽取 data1 中第一个变量数据 画出曲线
subplot(2,1,1)
plot(data(1,:)),title('x 在迭代中的变化')
%抽取 data 中的第二个变量数据 画出其变化曲线
subplot(2,1,2)
plot(data(2,:)),title('y 在迭代中的变化')

function [allx,ally,r,n]=MultiNewton(F,x0,eps)
% allx是用来记录每一步迭代的点的矩阵，ally是用来记录迭代点对应计算出来的函数值的，
% r是满足精度要求的最终迭代点 n迭代次数
  if nargin==2
    eps=1.0e-4;
  end
  x0 = transpose(x0);
  Fx = subs(F,transpose(symvar(F)),x0);
  var = transpose(symvar(F));
  dF = jacobian(F,var);
  dFx = subs(dF,transpose(symvar(F)),x0);
  n=dFx;
  r=x0-inv(dFx)*Fx';
  n=1;
  tol=1;
  N=100;
  symx=length(x0);
  ally=zeros(symx,N);
  allx=zeros(symx,N);
  while tol>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)