使用正割法的 MATLAB 程序有时似乎输出无意义的值

MATLAB program using Secant Method sometimes seems to output nonsense values

为了学习 MATLAB,我正在开发一个使用正割法求零根的程序。在我的代码下方:

function [zero, x_vector] = secant_method(fx, x1, x2, nmax, tol)
%function that tries to find root of given function using the secant method. 
%Mandatory inputs: FX is the function to be evaluated, X1 is first input
%value of the function, X2 is second input value of the function.
%Optional inputs: NMAX is maximum number of iterations (default = 5000), tol
%is maximum tolerance (default = 1e-4).
%Outputs: ZERO is the root of the function within the tolerance or after
%the maximum iterations has been reached, X_VECTOR is a vector of all
%consecutive evaluted inputs of the function. 

%code by Steven

if ~exist('tol')
    tol = 1e-4;
end

if ~exist('nmax')
    nmax = 5000;
end

y_vector = [];
x_vector = [];
x_vector(1)=x1;
x_vector(2)=x2;

diff=tol+1;         %this ensures at least one iteration will take place
i=1;
while diff > tol && i < (nmax)
    i = i+1;

    %Creating new straight line from two previous coordinates, using the conventional
    % y=ax+b for paramter notation.
    y_vector(i-1) = feval(fx,x_vector(i-1));
    y_vector(i) = feval(fx,x_vector(i));
    a = (y_vector(i)-y_vector(i-1))/(x_vector(i)-x_vector(i-1));
    b = y_vector(i)-a*x_vector(i);

    x_vector(i+1) = -b/a; 

    diff = abs(x_vector(i+1)-x_vector(i));
end

if i == (nmax) && diff > tol
    fprintf(['Function stopped without converging to the desired tolerance\n',... 
        'because the number of maximum iterations was reached.']);
end

zero = x_vector(i+1);

return 

现在对于大多数输入,代码都可以正常工作。例如

secant_method('sin', 1, 5, 50, 1e-10) = 3.141592653589793 足够接近 pi.


但是对于一些输入,我不明白为什么代码 returns 它做了什么。例如

secant_method('sin', -8, 12, 500, 1e-10) = 2.002170932439574e+30

but sin(2.002170932439574e+30)= 0.019481942752609 which is not really close to 0.

非常感谢任何关于为什么会发生这种情况或我哪里出错的想法。

正割法规则

if fx(x1) < 0 ---> choose x2 such as fx(x2) > 0
if fx(x1) > 0 ---> choose x2 such as fx(x2) < 0

示例 1:

x1 = 1; sin(x1) =  0.8415;
x2 = 5; sin(x2) = -0.9589;

fx(x1) > 0; fx(x2) < 0 No problem

示例 2:

x1 = -8; sin(x1) =   -0.9894;
x2 = 12; sin(x2) = -0.5366;

fx(x1) < 0; fx(x2) < 0 Change the boundary values