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

误差和残量

数值求解方程\(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


练习

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\)所对应的迭代步数和相对误差填在一个表格内。

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