零点的存在性定理
早在高中阶段，我们就学习过函数的零点存在性定理。简单地说,对于区间[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小时内删除。
作 者 | 郭志龙
编 辑 | 郭志龙
校 对 | 郭志龙