二分法 | matlab 实现

function bisec(file_name, a, c, xmin, xmax, n_points)
% file_name 是要计算的函数
% a, c, b 是含根区间上下限和中间值
% xmin，xmax 是图形横坐标的最小和最大值
% n_points 是自变量X的采样数
% Y_a, Y_c 是当前端点的y值
% it 是迭代的次数


clf, hold off
clear Y_a, clear Y_c

wid_x = xmax - xmin;
dx = (xmax - xmin) / n_points;
xp = xmin: dx: xmax;
yp = feval(file_name, xp);
plot(xp, yp, 'r'); % y=f(x)
xlabel('x'); ylabel('f(x)'); title('Bisection Method');
hold on

ymin = min(yp);
ymax = max(yp);
wid_y = ymax - ymin;
yp = 0.0 * xp;
plot(xp, yp); % x坐标轴

fprintf('Bisection Scheme \n\n');
fprintf('it  a  b  c   fa = f(a)');
fprintf('  fc = f(c)  abs(fc - fa)  abs(c - a)/2\n');
tolerance = 0.000001; % 绝对误差限
it_limit = 30; % 迭代的最高次数
it = 0;
Y_a = feval(file_name, a);
Y_c = feval(file_name, c);
plot([a, a], [Y_a, 0], 'black');
text(a, -0.1 * wid_y, 'x=a'); 
plot([c, c], [Y_c, 0], 'black');
text(c, -0.3 * wid_y, 'x=c');
if(Y_a * Y_c > 0)
    fprintf('f(a)f(c)>0 \n');
else
while 1
    it = it + 1;
    b = (a + c) / 2;
    Y_b = feval(file_name, b);
    plot([b, b], [Y_b, 0], ":")
    plot(b, 0, 'o')
    if it<4
        text(b, wid_y/20, [num2str(it)])
    end
        fprintf('%3.0f %10.6f %10.6f', it, a, b);
        fprintf('%10.6f %10.6f %10.6f', c, Y_a, Y_c);
        fprintf(' %12.3e %12.3e\m', abs(Y_c-Y_a), abs(c-a)/2);
    if(abs(c-a)/2 <= tolerance)
        fprintf('Tolerance is satisfied. \n');
        break
    end
    if(it>it_limit)
        fprintf('Iteration limit exceeded. \n');
        break
    end
    if(Y_a * Y_b <= 0)
        c = b;
        Y_c = Y_b;
    else
        a = b;
        Y_a = Y_b;
    end
end
fprintf('Final result ; Root = %12.6f \n', b)
end
x = b;
plot([x,x], [0.05 * wid_y, 0.2 * wid_y], 'g');
text(x, 0.25 * wid_y, 'Final solution');
plot([x, (x-wid_x * 0.004)], [0.05 * wid_y, 0.09 * wid_y], 'r');
plot([x, (x+wid_x * 0.004)], [0.05 * wid_y, 0.09 * wid_y], 'r');
    

function y=ex_5x2(x)
y=exp(x) - 5.* x .^2