在 matlab 中，如何创建一个循环来重复我的 Newton Raphson 算法，直到找到所有根？

Question

我写了下面的代码，它使用 Newton Raphson 方法求根。它可以找到 1 个根，但随后会中断。应该如何修改代码以重复该算法，直到在所需范围内找到所有根？

我知道我应该使用 for 或 while 循环我只是想不通。

'''

function [xn nrfail] = newraphson(fun,xg,xl,xr,tol)
% fun is the function of interest = (0.01-1) + x*cot(x), xg is the initial root guess, xl and xr are the bounds of the range I am using
% Initializing
i=0;
nrfail=0;
check=1;
h=1e-4;

% Loop
    while tol<check % Tolerence check
        i=i+1;
        fp=derivative(fun,xg,h); %derivative is another function I wrote that finds the derivative of the function of interest
        y=feval(fun,xg);
        xn=xg-y/fp; % NR method

        if xn<xl || xn>xr || i>10 || fp==0 % Method check
            nrfail=1;
            break
        end

        if abs(xn) <= 1
            check=abs(xg-xn); % Absolute error
        else
            check=abs(1-xg/xn); % Relative error
        end

        xg=xn; % New guess
        end
    end 

'''

Answer 1

首先，无法确定 newton-raphson 方法会收敛初始输入。一旦函数的导数达到abs (f') <= errThreshold，算法就会收敛。这意味着算法失败

1. 非常“扁平”的功能或
2. 函数域中导数变平坦的区域，即局部平坦的函数。

   第二个问题是，对于具有多个最大值或最小值的函数，函数会在多个点变平。所以结果实际上取决于起点。 我给出了这些初始警告，因为除非您完全 100% 确定该函数定义得很好，否则您很可能会陷入无限循环。所以最好的方法是，而不是设置无限循环，允许maximumNumberOfIterations，之后算法“崩溃”并返回“找不到一些根' ...

   也就是说，根据您的计算，请考虑改用 bissection method。它的收敛速度往往较慢，但话又说回来，收敛对导数的依赖性较小，因此它往往会在所有情况下收敛。

   您正在寻找的算法如下（在 Python 上）：

   from numpy.random import rand
   
   #
   # The maximum number of iterations before the 
   # algorithm breaks.
   ##
   maximumNumberOfIterations : int
   
   #
   # The function that we are evaluating.
   ##
   function : Callable
   
   #
   # The lower bound.
   ##
   lowerBound : float;
   
   #
   # The upper bound.
   ##
   upperBound : float;
   
   #
   # The tolerance.
   ##
   tolerance : float;
   
   #
   # An array to store the outcomes 
   # for each iteration
   ##
   roots: list[float] = [] 
   
   #
   # For each iteration DO:
   ##
   for iteration in range (maximumNumberOfIterations):
   
       #
       # Generate a 'random' initial guess
       # specified within the allowed boundary
       ##
       initialGuess : float = (upperBound - lowerBound) * rand()
   
       #
       # Actually try to run the algorithm and append
       # the root to the list.
       ##
       try:
           roots.append (
               newraphson(function, guess, lowerBound,
               upperBound, tolerance)
           )
   
       #
       # It is likely that the algorithm will fail and
       # return a RuntimeFault for 'some' cases. We
       # must handle these... Not a big deal, just do
       # nothing and continue.
       ##
       except RuntimeError:
           pass
   
   #
   
   # Now we just need to handle the repeated roots. Being
   # lazy I can just use a set and display the results...
   ##
   roots = set (roots)
   
   for root in roots:
       print (roots)
       
   #
   # If you were expecting a certain number of roots ...
   ##
   
   expectedNumberOfRoots : int
   
   if len (roots) < expectedNumberOfRoots:
       print ({n} roots could not be found"
       .format (n=expectedNumberOfRoots-len (roots)))