bisect_vis.m

Go to the documentation of this file.

%%
clear, clc, close all
 
%% define function and intervals
syms x
% f(x) = sin(x) - 0.5 * x;  a = 0.01; b = 2;
% f(x) = sqrt(sin(x)) - 0.5 * x;  a = 0.01; b = 2;
f(x) = sqrt(sin(x)) - 0.5 * x^6;  a = 0.01; b = 2;
% f(x) = -atan(100*(x-1.23456789)); a = 0.1; b = 2;
% f(x) = (1.234567-x)*(x+0.987654); a = 1; b = 2;
 
%% define interval
a = 0.01;
b = 2;
 
%% graphics
f0=figure();
fplot(f,[a,b]); axis ([a,b,-1,2])
hold on
plot([a,b],[0,0],'-r')
 
% '$f(a,b,c)=',latex(f(a,b,c)),'$'
title(['f(x) ',char(f)]); xlabel("x")
legend("f(x)")
 
%%
fnum=matlabFunction(f);
% Define the tolerance for the bisection method:  Co-pilot
eps = 1e-6;
%% Call the bisection method to find the root of the function f:  Co-pilot
%rootApprox = Bisect(fnum, a, b, eps, false);
rootApprox = Bisect(fnum, a, b, eps, 0);
fprintf('Equation to solve: %s = 0. \n',char(f));
disp(['mid point: Approximate root: ', num2str(rootApprox)]);
fprintf('mid point:            root: %f  after %i iterations. \n',rootApprox(end),length(rootApprox));
%% Linear improved bisection
%rootLinear = Bisect(fnum, a, b, eps, true);
rootLinear = Bisect(fnum, a, b, eps, 1);
fprintf('Equation to solve: %s = 0. \n',char(f));
disp(['linear: Approximate root: ', num2str(rootLinear)]);
fprintf('linear:            root: %f  after %i iterations. \n',rootLinear(end),length(rootLinear));
 
%
rootQuad = Bisect(fnum, a, b, eps, 2);
disp(['quad: Approximate root: ', num2str(rootQuad)]);
fprintf('quad:            root: %f  after %i iterations. \n',rootQuad(end),length(rootQuad));
 
%%
f3=figure();
% Evaluate the function at the approximate root and plot the result
% semilogy(rootApprox, f_val, '*b', 'MarkerSize', 10);
semilogy(abs(rootApprox-rootApprox(end)), '-*m', 'MarkerSize', 10);
hold on
semilogy(max(eps/10,abs(rootLinear-rootLinear(end))), '-*b', 'MarkerSize', 10);
semilogy(max(eps/10,abs(rootQuad-rootQuad(end))), '-or', 'MarkerSize', 10);
 
title("Iteration history"); xlabel("iteration k"); ylabel("|x^\ast-x_k|")
legend("midpoint", "linear", "quadratic");
grid on;
 
 
%%
function xvec = Bisect(func,a,b,eps,method)
% f(a) > = 0 >= f(b)
% Initialize the bisection method:  Co-pilot
xvec = [];
fc = 2*eps;   % guarantee entry into the loop
while abs(fc) >= eps && (b - a) / 2 > eps
    if 2==method
        Visualize_Interval_subdivision(func,a,b)
    end
    fa=func(a); fb=func(b);
    method_loc = method;
    TFAC=2;
    if (1==method && (-fa/fb>TFAC || -fb/fa>TFAC) && length(xvec)<4)  % improve linear
        method_loc = 0;
        fprintf('switch method to bisect at [%f, %f]\n',a,b);
    end
    switch method_loc
        case 0                        % midpoint
            c = (a + b) / 2;
        case 1                        % linear
            c = a-(b-a)*fa/(fb-fa);
        case 2                        % quadratic
            c = (a + b) / 2;
            dd =[1 a a^2; 1 b b^2; 1 c c^2]\[fa;fb;double(func(c))];
            cr = getRoot(dd,a,b);
            if ~isempty(cr)
                c = cr(1);
            else
                disp("no root found")
            end
        otherwise
            disp("Unknown method!  Choose from [0,2]")
    end
    fc = func(c);
    if abs(fc) < eps
        % % we found the root
        % xvec = [xvec, c]; % Exact root found
        % break;
    elseif fa * fc < 0
        b = c; % Root is in the left half
    else
        a = c; % Root is in the right half
    end
    xvec = [xvec, c]; % Store the midpoint
end
end
 
%%
function r = getRoot(dd,a,b)
r = roots(dd(end:-1:1));
r = r(a<=r & r<=b);
if length(r)==2
    disp("!!  2 roots  !!")
end
end
%%
function Visualize_Interval_subdivision(func,a,b)
%
fa=func(a); fb=func(b);
 
c0 = (a + b) / 2;              % midpoint
c1 = a-(b-a)*fa/(fb-fa);       % linear
 
% c2 = c0;              % quadratic via midpoint
c2 = c1;              % quadratic via linear
dd =[1 a a^2; 1 b b^2; 1 c2 c2^2]\[fa;fb;double(func(c2))];
cr = getRoot(dd,a,b);
if ~isempty(cr)
    c2 = cr(1);
else
    disp("no root found")
end
 
figure()
plot([a,b],[fa,fb],'sg')
hold on
plot(c0,func(c0),'*m')
plot(c1,func(c1),'db')
plot(c2,func(c2),'or')
%
x = linspace(a,b,101);
plot(x,func(x),'-g')
plot([a,b],[fa,fb],'-b')
y = polyval(dd(end:-1:1),x);
plot(x,y,'-r')
plot([a,b],[0,0],'-k')
title("Visualize 3 strategies in bisection step")
legend("given data","midpoint","linear","quadratic")
 
end
 
%%
function xvec = Bisect_org(func,a,b,eps,linear)
% Initialize the bisection method:  Co-pilot
xvec = [];
while (b - a) / 2 > eps
    if linear
        midpoint = double(a-(b-a)*func(a)/(func(b)-func(a)));
    else
        midpoint = (a + b) / 2;
    end
    % if subs(func, x, midpoint) == 0
    % if func(midpoint) == 0
    if abs(func(midpoint)) < eps
        xvec = [xvec, midpoint]; % Exact root found
        break;
        % elseif subs(func, x, a) * subs(func, x, midpoint) < 0
    elseif func(a) * func(midpoint) < 0
        b = midpoint; % Root is in the left half
    else
        a = midpoint; % Root is in the right half
    end
    xvec = [xvec, midpoint]; % Store the midpoint
end
end