Newton-Raphson算法简介及其R实现

本文简要介绍了Newton-Raphson方法及其R语言实现并给出几道练习题供参考使用。 下载PDF格式文档（Academia.edu）

• Newton-Raphson Method
  Let $f(x)$ be a differentiable function and let $a_0$ be a guess for a solution to the equation $$f(x)=0$$ We can product a sequence of points $x=a_0, a_1, a_2, \dots $ via the recursive formula $$a_{n+1}=a_n-\frac{f(a_n)}{f'(a_n)}$$ that are successively better approximation of a solution to the equation $f(x)=0$.
• R codes
  
  There are 4 parameters in this function:
  • f is the function you input.
  • tol is the tolerance (default $1e-7$).
  • x0 is the initial guess.
  • N is the default number (100) of iterations.
  The process will be end up until either the absolute difference between two adjacent approximations is less than tol, or the number of iterations reaches N.
• Examples
  Generally speaking, the "guess" is important. More precisely, according to Intermediate Value Theorem we can find two values of which function value are larger and less than 0, respectively. Then choosing the one, which first derivative is larger than another, as the initial guess value in the iterative formula. This process will guarantee the convergence of roots. Let's see some examples.
  • Example 1
    Approximate the fifth root of 7.
    Solution:
    Denote $f(x)=x^5-7$. It is easily to know that $f(1)=-6 < 0$ and $f(2)=25 > 0$. Additionally, $f'(1)=5 < f'(2)=80$, so we set the initial guess value $x_0=2$. By Newton-Raphson method we get the result is 1.47577316159. And $$f(1.47577316159)\approx 1.7763568394e-15$$ which is very close to 0. R codes is below:

    # Example 1
    f = function(x){x^5 - 7}
    h = 1e - 7
    df.dx = function(x){(f(x + h) - f(x)) / h}
    df.dx(1); df.dx(2)
    # [1] 5.0000009999
    # [1] 80.0000078272
    app = newton(f, x0 = 2)
    app
    # [1] 1.68750003057 1.52264459615 1.47857137506 1.47578373325 1.47577316175
    # [6] 1.47577316159
    f(app[length(app)])
    # [1] 1.7763568394e-15

  • Example 2
    The function $f(x)=x^5-5x^4+5x^2-6$ has a root between 1 and 5. Approximate it by Newton-Raphson method.
    Solution:
    We try to calculate some values first. $f(1)=-5, f(2)=-34, f(3)=-123, f(4)=-182, f(5)=119$, so there should be a root between 4 and 5. Since $f'(4)=40 < f'(5)=675$, hence $x_0=5$ is a proper initial guess value. By Newton-Raphson method we get the result is 4.79378454069 and $$f(4.79378454069)\approx -2.84217094304e-14$$ which is a desired approximation. R codes is below:

    # Example 2
    f = function(x){x^5 - 5 * x^4 + 5 * x^2 - 6}
    x = c(1 : 5)
    f(x)
    # [1]   -5  -34 -123 -182  119
    h = 1e-7
    df.dx = function(x){(f(x + h) - f(x)) / h}
    df.dx(4); df.dx(5)
    # [1] 40.0000163836
    # [1] 675.000053008
    app = newton(f, x0 = 5)
    app
    # [1] 4.82370371755 4.79453028339 4.79378501861 4.79378454069 4.79378454069
    f(app[length(app)])
    # [1] -2.84217094304e-14

  • Example 3
    A rectangular piece of cardboard of dimensions $8\times 17$ is used to make an open-top box by cutting out a small square of side $x$ from each corner and bending up the sides. Find a value of $x$ for which the box has volume 100.
    Solution:
    Firstly, building the model. $V(x)=x(8-2x)(17-2x)=100$, that is, we want to find the root of equation $$f(x)=x(8-2x)(17-2x)-100=0\Leftrightarrow f(x)=4x^3-50x^2+136x-100=0$$ We know that $0 < x < 4$ and hence try to calculate some non-negative integers: $$f(0)=-100, f(1)=-10, f(2)=4, f(3)=-34, f(4)=-100$$ Note that there are two intervals may have roots: $(1, 2)\cup (2,3)$. Since $$f'(1)=48 > f'(2)=-16 > f'(3)=-56$$ so we set the initial guess values $x_0=1$ and $x'_0=2$ (i.e. there are two separate iteration procedures). By using Newton-Raphson method we obtain the result are 11.26063715644 and 2.19191572127 respectively. Both of them are quite accurate. R codes is below:

    # Example 3
    f = function(x){4 * x^3 - 50 * x^2 + 136 * x - 100}
    x = c(0 : 4)
    f(x)
    # [1] -100  -10    4  -34 -100
    h = 1e-7
    df.dx = function(x){(f(x + h) - f(x)) / h}
    df.dx(1); df.dx(2); df.dx(3)
    # [1] 47.9999962977
    # [1] -16.0000024607
    # [1] -56.0000012229
    app1 = newton(f, x0 = 1)
    app2 = newton(f, x0 = 2)
    app1; app2
    # [1] 1.20833334940 1.25768359879 1.26062673622 1.26063715631 1.26063715644
    # [1] 2.24999996155 2.19469026652 2.19192282154 2.19191572132 2.19191572127
    f(app1[length(app1)]); f(app2[length(app2)])
    # [1] 2.84217094304e-14
    # [1] -2.84217094304e-14