方程数值求解
下面几讲，我们将聚集如下方程的解法：
\begin{equation}
f(x)=0
\tag{3.1}\label{3.1}
\end{equation}
在微积分课程中，我们知道，许多优化问题最终归结为求解上述形式的方程，其中\(f\)为你要求极值的函数\(F\)的导数。在工程问题中，函数\(F\)来源多种多样，有公式、微分方程的解、实验和模拟等。
牛顿迭代
我们把方程\eqref{3.1}的解记为\(x^\*\)。方程的解法有三种：对分法、割线法和牛顿法。这三种方法都需要猜测接近\(x^\*\)的初始值。
我们这里讨论牛顿法。牛顿法的基础是将函数在某点做线性化近似，即
\begin{equation}
f(x) \approx f(x_0) + f'(x_0)(x-x_0)
\tag{3.2}\label{3.2}
\end{equation}
我们要求出满足\(f(x)=0\)的\(x\)，因此我们令上式左边\(f(x)\)为0，求出\(x\)：
\begin{equation}
x\approx x_0-\frac{f(x_0)}{f'(x_0)}
\tag{3.3}\label{3.3}
\end{equation}
把所得\(x\)记为\(x_1\)，将\(x_1\)带入上式，得\(x_2\)，再将\(x_2\)带入上式，得\(x_3\)，……，最后得数列\(\{x_0,x_1,x_2,\cdots\}\)，数列满足如下关系
\begin{equation}
x_{i+1}\approx x_i-\frac{f(x_i)}{f'(x_i)}
\tag{3.4}\label{3.4}
\end{equation}
如果\(f(x)\)在\(x^\*\)附近没有特殊的奇异性，且\(x_0\)接近\(x^\*\)，上述数列将很快收敛到\(x^\*\)。
循环：for ... end 语句
应用牛顿法，我们需要按照\eqref{3.4}式多次反复计算。这在程序中可由循环来实现。下面我们看一个for循环的简单例子：
function S = mysum (n)
% gives the sum of the first n integers
S = 0; % start at zero
% The loop :
for i = 1:n % do n times
S = S + i; % add the current integer
end % end of the loop
这个函数是为了计算前\(n\)个整数的和。在命令窗口运行这个函数
>> mysum(100)
得到从1加到100的和。
下面我们用for循环实现牛顿法。程序如下：
function x = mynewton (f,f1 ,x0 ,n)
% 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
% n -- the number of steps to do
% Output : x -- the approximate solution
format long % prints more digits
format compact % makes the output more compact
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 , prints x too
end
在命令窗口，定义内联函数\(f(x)=x^3-5\)：
>> f = inline('x.^3-5','x')
并定义其导数\(f1\)：
>> f1 = inline('3*x.^2','x')
运行牛顿法，
>> mynewton(f,f1,2,4)
ans =
1.709975946676833
对这个函数应用牛顿法。方程的根是\(5^{1/3}\)。根据你运行的结果，评估一下，程序给出的结果与真实值相差多少？程序收敛，需要\(n\)最小为多少？
收敛
当\(f'(x^\*)\)为非零有限值，且\(x_0\)很接近\(x^\*\)时，牛顿法会很快收敛。下面我们看看啥时候会出问题。
对于\(f(x)=x^{1/3}\)，可知\(x^\*=0\)，\(f'(x^\*)=\infty\)。如果你依次输入和运行以下命令
>> f = inline('x^(1/3)','x')
>> f1= inline('x^(-2/3)/3','x')
>> mynewton(f,f1,0.1,10)
得到的结果为复数：
ans =
1.023999999999999e+02 - 1.617981985233248e-13i
再看一个例子，\(f(x)=x^2\)，\(f'(x)=2x\)，\(x^*=0\)，\(f'(x^*)=0\)，依次输入和运行如下命令：
>> f = inline('x^2','x')
>> f1 = inline('2*x','x')
>> mynewton(f,f1,1,10)
能得到结果，但是，收敛很慢。
>> mynewton(f,f1,1,10)
ans =
9.765625000000000e-04
>> mynewton(f,f1,1,100)
ans =
7.888609052210118e-31
如果初始值\(x_0\)与\(x^\*\)偏离太远，式\eqref{3.2}不成立，按式\eqref{3.4}迭代1次的结果\(x_1\)可能远离也可能接近方程的根\(x^\*\)。继续迭代下去，会出现两种情况：
- \(x_n\)收敛到\(x^\*\)
- 迭代发散，不会收敛到\(x^\*\)
练习
3.1 用牛顿法计算\(f(x)=x^5-7\)，初值为\(x_0=2\)。保证程序收敛（结果不再变化），\(n\)最小为多少？计算误差和残量。误差和残量等于0吗？
3.2 考虑一球从2米高处下落，碰到一坚硬地面反弹，碰撞恢复系数为0.9。写一个注释良好的脚本程序，计算球与地面第\(n\)次碰撞时，球走过的路程。碰撞多少次后，球走过的距离不再变化。
3.3 对于函数$ f(x) = x^3 − 4\(，以\)x_0=2$为初值，用纸和计算器进行牛顿迭代，计算每次迭代的误差和相对误差。保留足够多的小数，以免错以为收敛。将结果填在表内。