图

多区间二分法

预备知识 二分法

   这里介绍一种多区间二分法, 可以求出连续函数在某区间内几乎全部的根. 方法就是把这个区间等分为若干个相等的小区间, 然后分别判断这些小区间两端函数值的符号, 对所有两端异号的区间使用二分法即可. 显然, 小区间的个数越多, 越有可能找到所有的根. 例程如下.

  bisectionN.m

1
2
3
4
5
6
7
8
9
10
11
12
13
function roots = bisectionN(f, int, N)
x = linspace(int(1), int(2), N); % 划分区间
y = arrayfun(f, x); % 求所有 f(x(ii))
figure; plot(x, y, '.-') % 画图
title('f(x)')
Sign = sign(y);
ind = find(Sign(1:end-1) .* Sign(2:end) <= 0); % 找符合条件的区间序号
Nroot = numel(ind);
roots = zeros(1, Nroot); % 预赋值
for ii = 1:Nroot
    roots(ii) = fzero(f, [x(ind(ii)),x(ind(ii)+1)]);  
end
end

   函数的前两个输入变量分别是需要求根的函数句柄和求根区间(二元行矢量或列矢量), 第三个变量 $N$ 是子区间端点的个数(即子区间的个数加一). 函数中先求出所有的端点 x, 以及对应的函数值 y, 然后画图. 第 6-7 行寻找所有两端异号或有一端为 0 的区间的序号, 然后在第 10 行的循环中对这些区间逐个使用二分法. 为了提高运算效率, 这里并没有使用“二分法” 中的例程, 而是使用了 Matlab 自带的 fzero 函数.

   bisectionN 的画图功能是为了让用户判断是否有可能出现漏根, 以下举两个例子说明.

1
2
3
>> f = @(x)exp(-0.2*x)*sin(x);
>> roots = bisectionN(f, [0, 15], 50)
roots = 0    3.1416    6.2832    9.4248   12.5664

图
图1:运行结果

   运行结果如图 1 , 由于画出的曲线较为光滑, 可判断漏根的可能性很小. 再看另一个例子

1
2
3
>> f = @(x)sin(1/x);
>> roots = bisectionN(f, [0, 0.3], 50)
roots = 0.0245  0.0398  0.0455  0.0531  0.0637  0.0796  0.1061  0.1592

图
图2:运行结果

   我们已经知道函数 $\sinRound{1/x}$ 在该区间上有无数个根, 且越接近 $x = 0$, 相邻根之间的距离越小. 运行结果如图 2 , 可见在区间 $[0, 0.1]$ 内, 子区间端点的函数值非常不平滑, 极有可能出现漏根. 为了求得更多的根, 我们可以增加子区间的个数.

致读者: 小时物理百科一直以来坚持所有内容免费且不做广告,这导致我们处于日渐严重的亏损状态。长此以往很可能会最终导致我们不得不选择商业化,例如大量广告,内容付费,会员制,甚至被收购。因此,我们鼓起勇气在此请求广大读者热心捐款,使网站得以健康发展。如果看到这条信息的每位读者能慷慨捐助 10 元,我们几天内就能脱离亏损状态,并保证网站能在接下来的一整年里向所有读者继续免费提供优质内容。感谢您的支持。
—— 小时(项目创始人)

编辑词条 返回目录 返回主页 捐助项目 © 小时物理百科 保留一切权利