数值分析学习(一)之交叉法(Bracketing Methods)求根

秒速五厘米 2022-08-06 08:27 108阅读 0赞

在使用计算机求根的时候，由于计算机的数字离散性，通常需要通过试一试的方法确定一个初始值，然后根据这个初始值重复迭代过程，使其最后的函数值趋近于0。但是更多的时候，我们需要让这个过程在计算机中自动进行。目前使用普遍的初始值的猜想方法主要有两个：

1、交叉法(Bracketing Methods)：基于两个猜想的初始值坐落在两边，假想根在这两个初始值的中间；

2、开型法(Open Methods)：可以涉及到一个或多个初始值的猜想，但是没必要将初始值放在实际的根的两边；

比较两种方法，Bracketing Methods几乎在所有情况下都能使用，但是收敛速度较慢，需要更多的迭代次数；而Open Methods在有些情况下不能使用，但是收敛速度较快。

下面首先介绍几种Bracketing Methods：

1、增量搜索(Intcremental Search)

假如两个上下初始值有式子 ***f******(******x******l******)******f******(******x******u******) <*** ***0，*xl**表示较小的初始值猜想，**xu**表示较大的初始值猜想。则至少有一个实根在**xl**和**xu**之间。增量搜索有一个搜索精度问题：通常情况下，值得确定出现在***f******(******x******l******)******f******(******x******u******)***的值改变方向的时候，但是如果每次搜索的步长太小，就会导致迭代次数过多，收敛时间很长；如果搜索步长较大就可能导致空间上最接近，最精确的那个跟被丢掉。如下图：  
![20150108092438880][]

上图中最右边的根有可能因为搜索步长太大而被丢失。

下面是Matlab的Incremental Search函数：

function xb = incsearch(func,xmin,xmax,ns)
    % incsearch: incremental search root locator
    % xb = incsearch(func,xmin,xmax,ns):
    % finds brackets of x that contain sign changes
    % of a function on an interval
    % input:
    % func = name of function
    % xmin, xmax = endpoints of interval
    % ns = number of subintervals (default = 50)
    % output:
    % xb(k,1) is the lower bound of the kth sign change
    % xb(k,2) is the upper bound of the kth sign change
    % If no brackets found, xb = [].
    
    if nargin < 3
        error('At least 3 arguments required');
    end
    
    if nargin < 4
        ns = 50; 
    end
    
    x = linspace(xmin,xmax,ns);
    f = func(x);
    nb = 0;
    xb = [];
    
    for k = 1 : length(x) - 1
        if sign(f(k)) ~= sign(f(k + 1))
            nb = nb + 1;
            xb(nb,1) = x(k);
            xb(nb,2) = x(k + 1);
        end
    end
    
    if isempty(xb)
        disp('No brackets found')
        disp('Check interval or increase ns')
    else
        disp('Number of brackets: ')
        disp(nb)
    end

举个例子：

***f**(**x**)* =  sin*(*10*x**)* \+  cos*(*3*x**)***

*输入:*

**incsearch(@(x) sin(10\*x)+cos(3\*x),3,6)** 得到如下答案:

Number of brackets:  
5  
ans =  
3.2449 3.3061  
3.3061 3.3673  
3.7347 3.7959  
4.6531 4.7143  
5.6327 5.6939  
这是根的位置图:

![20150108094042946][]

2、对分法(Bisection Method)

对分法如下图所示：首先猜想出两个初始值，找到两个初始值的中间位置，两个初始值分别与中间位置坐标的函数值乘积，如果函数值乘积小于0，则选择函数值乘积小于0的那一边进行迭代，以下的步奏以此类推。

![20150108094505775][]

function [root, fx, ea, iter] = bisect(func,xl,xu,es,maxit,varargin)
    % bisect: root location zeroes
    % [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):
    % uses bisection method to find the root of func
    % input:
    % func = name of function
    % xl, xu = lower and upper guesses
    % es = desired relative error (default = 0.0001%)
    % maxit = maximum allowable iterations (default = 50)
    % p1,p2,... = additional parameters used by func
    % output:
    % root = real root
    % fx = function value at root
    % ea = approximate relative error (%)
    % iter = number of iterations
    
    if nargin < 3
        error('At least 3 input arguments required');
    end
    
    test = func(xl,varargin{:}) * func(xu,varargin{:});
    
    if test > 0
        error('No sign change');
    end
    
    if nargin < 4 | isempty(es)
        es = 0.0001;
    end
    
    if nargin < 5 | isempty(maxit) 
        maxit = 50;
    end
    
    iter = 0;
    xr = xl;
    ea = 100;
    
    while(1)
        xrold = xr;
        xr = (xl + xu)/2;
        iter = iter + 1;
        if xr ~= 0
            ea = abs((xr - xrold)/xr) * 100;
        end
    
        test = func(xl,varargin{:}) * func(xr,varargin{:});
    
        if test < 0
            xu = xr;
        elseif test > 0
             xl = xr;
        else
             ea = 0;
        end
    
        if ea <= es | iter >= maxit
            break;
        end
    end
    
    root = xr;
    fx = func(xr,varargin{:});

这里同时定义了相对误差：

![20150108101031380][]

x\_new是当前迭代值，x\_old是前一次的迭代值。迭代的时候需要满足一定的的误差条件，只要满足一定的误差范围，可以中断迭代，得到最终的值。

3、试位法(False Position Method)，也可以叫线性插值法(Linear Interpolation Method)

![20150108102130015][]

通过函数值上面的两个点画出直线，与**x**轴的交点找到**xr**。根据**xr**得到一个函数值，并且用**xr**替换与**xl**和**xu**中函数值符号相同的端值。以此进行迭代。**xr**的计算方程如下：

![20150108102707566][]

[20150108092438880]: /images/20220806/bada5241aa9649ae9dadc216179e0e1f.png
[20150108094042946]: /images/20220806/1541f66b62e549eda67277deceaee9b3.png
[20150108094505775]: /images/20220806/1b5b490ef9fb4e37b296dd87f64ac799.png
[20150108101031380]: /images/20220806/8b4e4151b946496492603bd6010c9846.png
[20150108102130015]: /images/20220806/0d9d0e79abb842cb9c0f038fa1d43fbf.png
[20150108102707566]: /images/20220806/a049f030b8734ef199f427ab03f4903d.png