(* Example 1 *) <br /> f[x_] := x^2 - 1 ; <br /> fp[x_] := D[f[s], s] /. s -> x ; <br /> Plot[f[x], {x, 0, 2}] ;

[Graphics:HTMLFiles/newton_nb_2.gif]

(* Example 2 *) <br /> f[x_] := x^2 ; <br /> fp[x_] := D[f[s], s] /. s -> x ; <br /> Plot[f[x], {x, -1, 1}] ;

[Graphics:HTMLFiles/newton_nb_4.gif]

(* Example 3 *) <br /> f[x_] := Sign[x] * Abs[x]^(1/3) ; <br /> fp[x_] := D[f[s], s] /. s -> x ; <br /> Plot[f[x], {x, -1, 1}] ;

[Graphics:HTMLFiles/newton_nb_6.gif]

(* initial guess, step size, number of iterations *) <br /> x = 1.5 ; <br /> del = 1 ; <br /> n = 10 ;

(* n steps symbolically computed *) <br /> outlist = Table[<br />          {i, x -= del * f[x]/fp[x]}, <br />          {i, 1, n}] ;

(* n steps numerically computed *) <br /> outlist = Table[<br />          {i, x -= del * f[x]/fp[x] // N}, <br />          {i, 1, n}] ;

Print[TableForm[outlist]] ;
1 0.75`
2 0.375`
3 0.1875`
4 0.09375`
5 0.046875`
6 0.0234375`
7 0.01171875`
8 0.005859375`
9 0.0029296875`
10 0.00146484375`

(* check for Big O rate of convergence *) <br /> p = 2.1 ; <br /> For[i = 1, i <= n, i ++, <br />      Print[Abs[outlist[[i, 2]] - 1] * n^p] ; <br />     ] ;

<br /> (* Check for order/asymptotic rate of convergence *) <br /> (* Example 1 : Quadratic convg, lambda = 1/2 *) <br /> (* Example 2 : Linear convg,      lambda = 1/2 *) <br /> p = 0 ; a = 1 ; <br /> For[i = 1, i <= n - 1, i ++, <br />      <br />      If[(Abs[outlist[[i, 2]] - p] != 0) && <br />               (Abs[outlist[[i + 1, 2]] - p] != 0), <br />              Print[i, ": ", Abs[outlist[[i + 1, 2]] - p]/    <br />                              Abs[outlist[[i, 2]] - p]^a], <br />            Print[i, ": p_n = p to machine accuracy"] <br />         ] ; <br />     ] ;

1 : 0.5`

2 : 0.5`

3 : 0.5`

4 : 0.5`

5 : 0.5`

6 : 0.5`

7 : 0.5`

8 : 0.5`

9 : 0.5`

Converted by Mathematica  (January 8, 2003)