Todd's Matlab讲义第5讲:二分法和找根

二分法和if ... else ... end 语句

先回顾一下二分法。要求方程(f(x)=0)的根。假设(c = f(a) < 0)(d = f(b) > 0),如果(f(x))是连续函数,那么方程的根(x^*)一定位于(a)(b)之间。然后,我们看一下(a)(b)中点(x=(a+b)/2),计算函数值(y=f(x)),如果函数不为0,比较(c)(d)(y)的符号,确定新的二分区间。具体来说,如果(c)(y)同号,新的二分区间就是([x,b]),如果(c)(y)异号,新的二分区间就是([a,x])。如图5.1所示。


图5.1 二分法

在不同的情况下做不同的事情,在程序里,这叫做流程控制。达到这个目的,最通常的方法是通过if ... else ... end语句来实现,这个语句是我们前面讲过的if ... end语句的扩展。

误差界

二分法的一个好处是,我们一直知道方程的真实的根(x^*)一定位于二分区间内,这样我们就可以知道最大误差可以是多少,一定小于区间大小的一半,即

egin{equation*} mid x-x^* mid lt frac{b-a}{2} end{equation*}

其中(x=(a+b)/2),为区间的中点。

二分法程序如下:

function [x e] = mybisect (f,a,b,n)
% function [x e] = mybisect (f,a,b,n)
% Does n iterations of the bisection method for a function f
% Inputs : f -- an inline function
% a,b -- left and right edges of the interval
% n -- the number of bisections to do.
% Outputs : x -- the estimated solution of f(x) = 0
% e -- an upper bound on the error
format long
% evaluate at the ends and make sure there is a sign change
c = f(a); d = f(b);
if c*d > 0.0
error (’Function has same sign at both endpoints .’)
end
disp (’ x y’)
for i = 1:n
% find the middle and evaluate there
x = (a + b )/2;
y = f(x);
disp ([ x y])
if y == 0.0 % solved the equation exactly
e = 0;
break % jumps out of the for loop
end
% decide which half to keep , so that the signs at the ends differ
if c*y < 0
b=x;
else
a=x;
end
end
% set the best estimate for x and the error bound
x = (a + b )/2;
e = (b-a )/2;

此程序不仅给出了方程的近似根,还给出了最大可能误差。

二分法总是可以进行下去,而牛顿法,如果初值(x_0)不足够接近(x^*)的话,程序可能不会收敛。而二分法,每一步的区间都可以缩一半,最终想缩多小就能缩多小。

找根

二分法和牛顿法都可以用来找根的任意近似值,但都需要先给定初值。二分法需要两个初值(a)(b),并要求真实的根位于二者之间,牛顿法只需要一个初值(x_0),要求不能离真实的根太远。如何确定合适的初值?看情况而定。如果你一次只解一个方程,最好先把函数画出来,通过图上的函数零点,可以很容易定出合适的初值。

有些情况下你不是一次就解一个方程,而是多次求解同一个方程,只是系数不同。当你开发某特定用途的软件的时候,常常会遇到这种情况。在这种情况下,首先你要利用问题的自然条件,比如选定的区间要符号实际情况。然后通过区间内某些点的符号,很容易就可以逐渐接近方程的真实的根了。方程的根一定位于符号相异的两点之间。下面的程序就是在一给定的区间([a_0,b_0])之间找方程的根:

function [a,b] = myrootfind (f,a0 ,b0)
% function [a,b] = myrootfind (f,a0 ,b0)
% Looks for subintervals where the function changes sign
% Inputs : f -- an inline function
% a0 -- the left edge of the domain
% b0 -- the right edge of the domain
% Outputs : a -- an array , giving the left edges of subintervals
% on which f changes sign
% b -- an array , giving the right edges of the subintervals
n = 1001; % number of test points to use
a = []; % start empty array
b = [];
% split the interval into n -1 intervals and evaluate at the break points
x = linspace (a0 ,b0 ,n);
y = f(x);
% loop through the intervals
for i = 1:(n -1)
if y(i)*y(i +1) < 0 % The sign changed , record it
a = [a x(i)];
b = [b x(i +1)];
end
end
if size (a ,1) == 0
warning (’no roots were found ’)
end

如果没有任何已知信息来找方程的根,事情就比较难办,好在这种情况在工程问题里不常见。

一旦确定出根所在的区间([a, b])(a)(b)就可以作为二分法和割线法(见下一讲)的初值。对于牛顿法,初值(x_0)一个明显的选法是(x_0 = (a + b)/2)。一个更好的选法是利用割线法选(x_0)

练习

1 修改程序mybisect,在给定公差下解方程。应用你的程序,求解方程(f(x) = 2x^3 + 3x − 1=0)的根,初始区间是([0,1]),公差为(10^{-8})。程序需要运行多少步能达到公差的要求?最终残量有多大?

2 用纸和计算器对于函数(f(x) = x^3 − 4)和初始区间([1,3])进行3步二分法迭代,计算每步结果的误差和相对误差,并与第3讲练习3的结果比较。

原文地址:https://www.cnblogs.com/joyfulphysics/p/5733464.html