Skip to content

Joyful Physics

格软物质 推硬物理

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

Todd’s Matlab讲义第4讲:控制误差和条件语句

瞿立建 2016-07-30

误差和残量

数值求解方程\[f(x)=0\]的根,有多种方法测算结果的近似程度。最直接的方法是计算误差。第\[n\]步迭代结果与真值$x^*$的差即为第$n$步迭代的误差:

\begin{equation*}
e_n=x_n-x^*
\end{equation*}

但是,我们一般是不知道真实值$x^*$的,否则,我们也不会费劲去算了。所以,直接计算误差是不可能的,需要我们另辟蹊径。

一个可能的方法是,程序一直运行,直到结果不再变化。这个方法通常还是很管用的。有时候,程序结果不再变化并不意味着$x_n$接近真实值,这时候这个方法就不灵了。

对于牛顿法,我们也可以采用这样的方法:每一步结果有效数字位数加倍。这个方法在程序里难以应用。

在很多情况下,我们不是测算$x_n$逼近$x^*$的程度,而是采用另外一个很实用的方法,计算$x_n$满足方程的程度,即计算$y_n=f(x_n)$与0接近的程度,用$r_n=f(x_n)-0$来表征,这个量称为残量(residual)。大多数情况下,我们只关心$r_n$的绝对值,所以我们只需要考察$\mid r_n \mid=\mid f(x_n)\mid$。

if … end 语句

我们设定一个公差$\mid r_n \mid=\mid f(x_n)\mid$,然后我们通过if … end 语句,将收敛条件包括在牛顿法程序里。程序如下:

function x = mynewton (f,f1 ,x0 ,n,tol )
% Solves f(x) = 0 by doing n steps of Newton ’s method starting at x0.
% Inputs : f -- the function , input as an inline
% f1 -- it ’s derivative , input as an inline
% x0 -- starting guess , a number
% tol -- desired tolerance , prints a warning if |f(x)|> tol
% Output : x -- the approximate solution
x = x0; % set x equal to the initial guess x0
for i = 1:n % Do n times
x = x - f(x)/ f1(x) % Newton ’s formula
end
r = abs (f(x))
if r > tol
warning (’The desired accuracy was not attained ’)
end

程序里,if … end 为条件语句。if 后面的条件 abs(y) > tol 如果满足,则执行后面的语句,如果不满足,程序直接跳至与if 对应的end。

下面我们在命令窗口操练一下。

>> f = inline('x.^3-5','x')
f =
     Inline function:
     f(x) = x.^3-5
>> f1 = inline('3*x.^2','x')
f1 =
     Inline function:
     f1(x) = 3*x.^2
>> mynewton(f,f1,2,5,0.01)
ans =
   1.709975946676697
>> mynewton(f,f1,2,5,1e-10)
ans =
   1.709975946676697
>> 

补充:如果选$n=3$,才会看到tol取得很小时出的警告信息。可能以前的版本算法落后,$n=5$还能看到警告信息。

循环: while … end 语句

前面这一版牛顿法程序可以告诉我们是不是收敛到指定精度的结果,但是我们还是需要输入迭代步数$n$。即便对于不太奇异的函数,如果步数$n$太小,也可能达不到指定的精度,然后就得增大$n$,重新计算。如果步数$n$太大,就会做不必要的迭代,浪费时间。

控制迭代步数的一个方法是,使迭代一直进行,直到残量$\mid r_n \mid =\mid f(x) \mid=\mid y \mid$足够小。在Matlab里,这可以通过 while … end 循环实现:

function x = mynewtontol (f,f1 ,x0 , tol )
x = x0; % set x equal to the initial guess x0
y = f(x);
while abs (y) > tol % Do until the tolerence is reached .
x = x - y/f1(x) % Newton ’s formula
y = f(x)
end

语句 while … end是个循环,与for … end类似,只是前者循环次数不是固定的,会一直进行到条件 abs(y) > tol满足为止。

上述程序的一个缺点是,abs(y)可能永远都不会比tol小,这就会导致程序会一直运行下去,直到我们手动终止。比如把公差设置的非常小:

>> tol = 10^(-100)

然后对函数$f(x)=x^3-5$再运行下该程序。

你可以用快捷键Ctrl-c结束程序。

要避免无限循环,可以再增加一个计数变量,比如i,控制迭代次数的上限。

于是,我们可写出如下程序:

function x = mynewtontol (f,f1 ,x0 , tol )
x = x0; % set x equal to the initial guess x0.
i =0; % set counter to zero
y = f(x);
while abs (y) > tol && i < 1000
% Do until the tolerence is reached or max iter .
x = x - y/f1(x) % Newton ’s formula
y = f(x)
i = i +1; % increment counter
end

练习

4.1 几何级数满足如下关系:

\begin{equation*}
\sum_{i=0}^{\infty}\frac{1}{r^i}=\frac{1}{1-r},\mid r \mid \lt \mid
\end{equation*}

下面是一个计算几何级数的脚本程序,但遗漏了一行,请补充完整,并验证程序是否可行。对于$r=0.5$,迭代步数需要达到多少程序才能收敛?结果与精确值2相差多少?

% Computes a geometric series until it seems to converge
format long
format compact
r = .5;
Snew = 0; % start sum at 0
Sold = -1; % set Sold to trick while the first time
i = 0; % count iterations
while Snew > Sold % is the sum still changing ?
Sold = Snew ; % save previous value to compare to
Snew = Snew + r^i;
i=i +1;
Snew % prints the final value .
i % prints the # of iterations .

在程序尾部增加一行,计算相对误差。计算$r = 0.9,\quad 0.99,\quad 0.999,\quad 0.9999,\quad 0.99999,\quad 0.999999$等情况,并将每个$r$所对应的迭代步数和相对误差填在一个表格内。

4.2 修改第3讲中的练习2的程序,计算当球跳起高度低于$1\mathrm {cm}$时,球走过的距离。要求用while 循环,不能用for 循环和if 语句。

分类目录 Todd讲Matlab 标签 for循环, if语句, while循环, 残量, 牛顿法, 误差, 迭代, 非线性方程
Previous: Todd’s Matlab讲义第3讲:牛顿法和for循环
Next: Todd’s Matlab讲义第5讲:二分法和找根

功能

  • 登录
  • 项目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