%=======================================================================
% ch11_false_position.m
%
% Find one root of f(x) = 0 by the false position (regula falsi) method.
%
% Origin:
%   Based on the false-position routine in Root5.m by Alexander Godunov
%   (February 2022). This standalone version uses the example and
%   status convention of ch11_false_position.f90.
%
% Example:
%   f(x) = sqrt(x) - cos(x) on [0,1].
%   Reference root = 0.64171437087288... .
%
% The program reports the iteration count, root estimate, and f(root).
% The requested tolerance applies to the residual |f(root)|. Unlike
% bisection, ordinary false position may leave one endpoint fixed for
% many iterations, so the bracket width is not the stopping test.
%=======================================================================

clear; close all; clc;

% Initial bracket and requested residual tolerance
x1  = 0.0;
x2  = 1.0;
eps = 1.0e-6;

fprintf('A single root of f(x)\n');
fprintf('Method - False position\n');

[root, flag] = false_position(@f, x1, x2, eps);

% Report a missing bracket or invalid input before printing a root
if flag == 0
    fprintf('invalid interval or no initial sign change\n');
    return
end

fprintf('iterations         = %d\n', abs(flag));
fprintf('residual tolerance = %.5e\n', eps);
if flag < 0
    fprintf('iteration limit reached before the requested residual\n');
end

fprintf('root               = %.5e\n', root);
% Check the residual of the reported estimate
fprintf('f(root)            = %.5e\n', f(root));


function [root, flag] = false_position(fun, x1, x2, eps)
%=======================================================================
% Locate one root of fun(x) on [x1,x2] by ordinary false position.
%
% Input:
%   fun   - function handle, continuous on the initial bracket
%   x1    - left endpoint of the initial bracket
%   x2    - right endpoint of the initial bracket
%   eps   - requested maximum residual |fun(root)| (positive)
%
% Output:
%   root  - most recent false-position estimate when flag is nonzero
%   flag  - convergence indicator and iteration count
%           flag > 0: |fun(root)| <= eps; flag is the iteration count
%           flag = 0: invalid interval/tolerance or no sign change
%           flag < 0: iteration limit reached; abs(flag) is the count
%
% Notes:
%   A sign change guarantees a root only when fun is continuous.
%   The maximum number of iterations is 200. One bracket endpoint
%   can remain nearly fixed, so a small interval width is not assumed.
%   An endpoint root is returned immediately with a positive flag.
%=======================================================================

iter = 200;
flag = 0;
root = NaN;  % No estimate is returned when the initial bracket is invalid.

% Check that the inputs describe a usable initial bracket
if x1 >= x2 || eps <= 0
    return
end

a = x1;
b = x2;
fa = fun(a);
fb = fun(b);
if fa*fb > 0
    return
end

% Handle an exact zero at an endpoint before forming a secant line.
% Flag 1 follows the existing convention that positive means success.
if fa == 0
    root = a;
    flag = 1;
    return
end
if fb == 0
    root = b;
    flag = 1;
    return
end

% Interpolate between the current endpoints and retain the sign change.
for i = 1:iter
    c = b - fb*(b-a)/(fb-fa);
    fc = fun(c);

    % False position uses the residual, not the bracket width, to stop.
    if abs(fc) <= eps
        root = c;
        flag = i;
        return
    end

    if fa*fc <= 0
        b = c;
        fb = fc;
    else
        a = c;
        fa = fc;
    end
end

root = c;
flag = -iter;
end


function y = f(x)
%-----------------------------------------------------------------------
% Example function. The chosen bracket must stay in the domain x >= 0.
% To try the alternative from the original file, use [-1,0].
%-----------------------------------------------------------------------
% y = x + cos(x);
y = sqrt(x) - cos(x);
end