Skip to content

Joyful Physics

格软物质 推硬物理

菜单1 菜单1
  • 首页
  • 关于
  • 连载
  • 友链
  • 导航
Expand Search Form

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

瞿立建 2016-08-03

二分法和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^*$一定位于二分区间内,这样我们就可以知道最大误差可以是多少,一定小于区间大小的一半,即

\begin{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$。

练习

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

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

分类目录 Todd讲Matlab 标签 if else语句, 二分法, 根, 误差, 迭代, 非线性方程
Previous: Todd’s Matlab讲义第4讲:控制误差和条件语句
Next: Todd’s Matlab讲义第6讲:割线法

功能

  • 登录
  • 项目feed
  • 评论feed
  • WordPress.org

近期文章

  • 美国物理学会流体视频大赛获胜视频 2017-10-03
  • 薄膜干涉 2017-04-10
  • 最速降线gif图片 2017-02-19
  • DNA 解链 2016-12-26
  • 生物化学机器DNA 2016-12-26

分类目录

  • Matlab (1)
  • Todd讲Matlab (6)
  • 博客教程 (3)
  • 备课讲义 (27)
    • 光学近代物理备课讲义 (7)
    • 电磁学备课讲义 (20)
  • 教学笔记 (23)
  • 未分类 (2)
  • 物理之外 (2)
  • 物理科普 (3)
  • 理论物理极础 (16)
  • 科学史哲 (6)
  • 科学时评 (31)
  • 科研笔记 (12)
  • 软物质物理 (7)
    • Doi高分子物理导论 (3)
  • 软物质科普 (9)

标签

latex PI random walk 世界一流 人物 假说 偏倚随机行走 凝胶 创新 加速度 势能 吉尔伯特 奥斯特 富兰克林 对称性 导体 导数 库仑定律 感应电动势 拉格朗日量 数学 无规行走 标度理论 梯度 泊松括号 法拉第 波 浓度梯度 理想链 电场 相空间 研究生 科学 科学家 简谐振动 能量守恒 蒲慕明 连续带电体 迭代 随机行走 静电屏蔽 非线性方程 高分子 高斯分布 高斯定理

近期评论

  • 蒲慕明所长在中科院神经科学研究所历年年会上的讲话 – Joyful Physics发表在《论研究生教育——蒲慕明所长在神经所2005年会上的讲话》
  • 术索发表在《磁介质》
  • 瞿立建发表在《友链》
  • taho发表在《磁介质》
  • taho发表在《友链》
2022年八月
日 一 二 三 四 五 六
 123456
78910111213
14151617181920
21222324252627
28293031  
« 10月    

文章归档

  • 2017年十月 (1)
  • 2017年四月 (1)
  • 2017年二月 (1)
  • 2016年十二月 (4)
  • 2016年十一月 (3)
  • 2016年十月 (10)
  • 2016年九月 (10)
  • 2016年八月 (10)
  • 2016年七月 (7)
  • 2016年六月 (8)
  • 2016年五月 (9)
  • 2016年四月 (7)
  • 2016年三月 (6)
  • 2016年二月 (7)
  • 2016年一月 (8)
  • 2015年十二月 (6)
  • 2015年十一月 (6)
  • 2015年十月 (10)
  • 2015年九月 (9)
  • 2015年八月 (9)
  • 2015年七月 (7)
  • 2015年六月 (7)

Joyful Physics © 2022