方程数值求解
下面几讲，我们将聚集如下方程的解法：
egin{equation}
f(x)=0
ag{3.1}label{3.1}
end{equation}
在微积分课程中，我们知道，许多优化问题最终归结为求解上述形式的方程，其中(f)为你要求极值的函数(F)的导数。在工程问题中，函数(F)来源多种多样，有公式、微分方程的解、实验和模拟等。
牛顿迭代
我们把方程eqref{3.1}的解记为(x^*)。方程的解法有三种：对分法、割线法和牛顿法。这三种方法都需要猜测接近(x^*)的初始值。
我们这里讨论牛顿法。牛顿法的基础是将函数在某点做线性化近似，即
egin{equation}
f(x) approx f(x_0) + f'(x_0)(x-x_0)
ag{3.2}label{3.2}
end{equation}
我们要求出满足(f(x)=0)的(x)，因此我们令上式左边(f(x))为0，求出(x)：
egin{equation}
xapprox x_0-frac{f(x_0)}{f'(x_0)}
ag{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})，数列满足如下关系
egin{equation}
x_{i+1}approx x_i-frac{f(x_i)}{f'(x_i)}
ag{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$为初值，用纸和计算器进行牛顿迭代，计算每次迭代的误差和相对误差。保留足够多的小数，以免错以为收敛。将结果填在表内。