牛顿求根公式（牛顿迭代法）（Octave/matlab 绘图）

贡献者： ACertainUser

图 1：牛顿迭代法求解 $f(x)=0.2x^2=0$ 示意图

　　 Newton 迭代法是求解方程（或函数零点）$f(x)=0$ 的经典方法。它的迭代公式是

\begin{equation} x_{k+1} = x_k-\frac{f(x_k)}{f'(x_k)} \qquad k=0,1,2,...~, \end{equation}

其中 $x_0$ 是人为选定的、对解的初始猜测。如果迭代收敛，那么迭代的结果将收敛于方程 $f(x)=0$ 的解 $x$。 $$\lim_{k\to+\infty} x_k = x~.$$

　　以下的 Octave/matlab 代码粗糙地实现了 Newton 迭代法，同时做出了过程示意图（类似图 1 ）：

clc
clear

f = @(x) 0.2*x.^2; %待求解的 f(x)=0
x(1) = 4;  
%初始猜测。由于octave/matlab的数组序号从1开始，因此我们用x1表示x0

hold on
axis equal

for i = 1:100
  if i ~=1
    dx = 0.001;
    fx = (f(x(i-1)+dx) - f(x(i-1)))/dx; 
    %数值近似x(i-1)处函数的导数值f'(x(i-1))。更好的办法是使用三点牛顿法。

    x(i) = x(i-1) - f(x(i-1))/fx; %Newton 迭代公式
    line([x(i) x(i-1)],[0 f(x(i-1))])
  end

  line([x(i) x(i)],[0 f(x(i))],'color','r')
  scatter(x(i),0);

  if i~=1 && abs(x(i) - x(i-1)) < 0.01 %如果误差已经足够小，则终止迭代
    break;
  end
end

if x(i) < x(1)
  _x = x(i)-0.5:0.01:x(1)+0.5;
else
  _x = x(1)-0.5:0.01:x(i)+0.5;
end

_y = f(_x);
plot(_x, _y,'k')

致读者：小时百科一直以来坚持所有内容免费，这导致我们处于严重的亏损状态。长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。因此，我们请求广大读者热心打赏 ，使网站得以健康发展。如果看到这条信息的每位读者能慷慨打赏 10 元，我们一个星期内就能脱离亏损，并保证在接下来的一整年里向所有读者继续免费提供优质内容。但遗憾的是只有不到 1% 的读者愿意捐款，他们的付出帮助了 99% 的读者免费获取知识，我们在此表示感谢。