MATLAB二分法求方程的根(实例加程序)

零点的存在性定理
早在高中阶段，我们就学习过函数的零点存在性定理。简单地说,对于区间[a,b]上的连续函数f(x),如果满足f(a)f(b)<=0，那么函数在[a,b]上至少存在一个零点。

根据函数与方程的关系我们可以得到，对于相应的方程f(x)=0。如果方程的左侧在a,b处不同号，那么，方程在[a,b]上存在零点。

二分法的思想
在得到根的存在性之后，我们就希望找到或者逼近方程的根。这种情况下比较显然的一种方式就是二分法。二分法的基本步骤如下：

实例1

程序

clc;
clear all;
close all;
Re = 1e4;%赋值Re的值
C = 0.57;%%赋值C的值
%第二问程序
f= @(beta) 0.5959+0.0312.*beta.^(2.1)-0.184*beta.^8+(91.71.*beta.^(2.5)./(Re^0.75))-C;%设置目标函数
a =0.9;%赋值a
b = 1;%赋值b
eps = 1e-6;%赋值eps
T = bisect(f,a,b,eps);%调用函数
data = [T(:,6) T(:,end)+C T(:,end)]; %输出求解的beta  C error
save('data93552.txt','data','-ascii');

function T = bisect(f,a,b,eps)
%%
%输入
%f代表输入的函数 a，b代表区间范围[a,b],eps是输入的误差
%T代表输出的参数 
%包括迭代次数  左区间a  a点函数值  右区间b b点函数值  区间a和b的中点值xk   xk点函数值
%%
k=1%设置初始值
x=(a+b)/2;%设置初始区间中点
fprintf('        k        a        f(a)       b        f(b)       xk        f(xk)\n ');%输出变量的名字
T=[k,a,f(a),b,f(b),x,f(x)];%对T赋值
while abs(T(k,4)-T(k,6))>eps/2   %判断和误差的大小
    k=k+1;%循环计数
    if  f(x)*f(a)==0   %判断当函数值为0的时候
        a=a;%左区间重新赋值
        b=x;%右区间重新赋值
        x=(a+b)/2;%区间中点重新赋值
        T=[T;k,a,f(a),b,f(b),x,f(x)];%对T赋值
        break
    elseif  f(x)*f(a)>0  %判断当f(x)和f(a)同号的情况
        a=x;%左区间重新赋值
        b=b;%右区间重新赋值
        x=(a+b)/2;%区间中点重新赋值
        T=[T;k,a,f(a),b,f(b),x,f(x)];%对T赋值
    elseif  f(x)*f(a)<0   %判断当f(x)和f(a)异号的情况
        a=a;%左区间重新赋值
        b=x;%右区间重新赋值
        x=(a+b)/2;%区间中点重新赋值
        T=[T;k,a,f(a),b,f(b),x,f(x)];%对T赋值
    end
 
end
disp(T);%输出变量T
fprintf('经过%d次迭代，函数方程根的近似解为：x=%.8f\n',k-1,T(k-1,6))%输出迭代过程
error = T(:,7);%误差
figure;%新建一个窗口
plot(1:k,error,'r');%画图
xlabel('k');%设置横轴坐标
ylabel('error value');%设置纵轴坐标
end

运行结果

结果：

k =

     1

        k        a        f(a)       b        f(b)       xk        f(xk)
     1.0000    0.9000    0.0422    1.0000   -0.0352    0.9500    0.0125
    2.0000    0.9500    0.0125    1.0000   -0.0352    0.9750   -0.0087
    3.0000    0.9500    0.0125    0.9750   -0.0087    0.9625    0.0025
    4.0000    0.9625    0.0025    0.9750   -0.0087    0.9688   -0.0029
    5.0000    0.9625    0.0025    0.9688   -0.0029    0.9656   -0.0002
    6.0000    0.9625    0.0025    0.9656   -0.0002    0.9641    0.0012
    7.0000    0.9641    0.0012    0.9656   -0.0002    0.9648    0.0005
    8.0000    0.9648    0.0005    0.9656   -0.0002    0.9652    0.0002
    9.0000    0.9652    0.0002    0.9656   -0.0002    0.9654    0.0000
   10.0000    0.9654    0.0000    0.9656   -0.0002    0.9655   -0.0001
   11.0000    0.9654    0.0000    0.9655   -0.0001    0.9655   -0.0000
   12.0000    0.9654    0.0000    0.9655   -0.0000    0.9655   -0.0000
   13.0000    0.9654    0.0000    0.9655   -0.0000    0.9654   -0.0000
   14.0000    0.9654    0.0000    0.9654   -0.0000    0.9654   -0.0000
   15.0000    0.9654    0.0000    0.9654   -0.0000    0.9654    0.0000
   16.0000    0.9654    0.0000    0.9654   -0.0000    0.9654    0.0000
   17.0000    0.9654    0.0000    0.9654   -0.0000    0.9654   -0.0000
   18.0000    0.9654    0.0000    0.9654   -0.0000    0.9654   -0.0000

经过17次迭代，函数方程根的近似解为：x=0.96543503

实例2

程序

clc;
clear all;
close all;
syms U L;    %将区间上下限定为变量
f=@(x) exp(x)-x^2+3*x-2;    %求给定的函数
U=1;    
L=0;
while U-L>1e-10   %设定精度
    root=(U+L)/2;    %当根的区间大于所给精度时，利用二分法重新规划求根区间
    if f(root)==0    
        break;    %r恰好为所求根，直接跳出循环
    end
    if f(root)*f(U)<0    %用零点存在定理判断根所在的区域
        L=root;
    else
        U=root;
    end
end
root  
%结果 root =0.2575

运行结果

%结果 root =0.2575

实例3

程序

clc;
clear all;
close all;

% -------------- inputs -------------------
f = @(x) 3*x^2-x-3;
a = 0;
b = 2;
% tolerance / max iter
TOL = 1e-4; NI = 50;
% -------------------------------------------------------
% STEP 1: initialization
i = 1;
fa = f(a);
converge = false; % convergence flag
% STEP 2: iteration
while i<=NI
% STEP 3: compute p at the i's step
p = a+(b-a)/2;
fp = f(p);
% STEP 4: check if meets the stopping criteria
if (abs(fp)<eps || (b-a)/2 < TOL) % eps is Matlab-machine zero
converge = true; % bisection method converged!
break; % exit out of while loop
else
% STEP 5
i = i+1;
% STEP 6
if fa*fp > 0
a = p; fa = fp;
else
b = p;
end
end
end
b
f(b)

运行结果

b =

    1.1805


ans =

   4.9619e-04

实例4

clc;
clear all;
close all;
a = 1;
b = 1.5;
tol = 1e-8;
x = half(a, b, tol)
function x = half(a, b, tol)% tol 是 tolerance 的缩写，表示绝对误差
c = (a + b) / 2; k = 1;
m = 1 + round((log(b - a) - log(2 * tol)) / log(2)); % <1>
while k <= m
    if f(c) == 0
        c = c;
    return;
    elseif f(a) * f(c) < 0
        b = (a + b) / 2;
    else
        a = (a + b) / 2;
    end
    c = (a + b) / 2; k = k + 1;
end
x = c; % 这里加分号是为了不再命令行中输出
k % 不加分号就会在控制台输出
c
end
function y = f(x)
y = x^3 - x -1;
end

运行结果

k =

    27


c =

    1.3247


x =

    1.3247

>> 

本文内容来源于网络，仅供参考学习，如内容、图片有任何版权问题，请联系处理，24小时内删除。

作 者 | 郭志龙
编 辑 | 郭志龙
校 对 | 郭志龙