2

In Wellin's Programming with Mathematica book, here's one of his implementations of the Newton method, where the iteration runs until the error tolerance is reached.

Mathematica graphics

Mathematica graphics

Clear[findRoot]
findRoot[expr_, {var_, init_}, ϵ_] :=
 Module[{xi = init, fun = Function[fvar, expr]},
  While[Abs[fun[xi]] > ϵ,
   xi = N[xi - fun[xi]/fun'[xi]]];
  {var -> xi}]

findRoot[x^3 - 2, {x, 2.0}, 0.0001]
(* {x -> 2.} *)

As you could see the result is clearly wrong. I think it's because of the presence of fvar in the body of Function, which was never defined. I think he meant to use var. I tried that, and it works, but there was a warning that "The variable name has been used twice in a nested scoping construct, in a way that is likely to be an error", with var highlighted in red.

Mathematica graphics

Should I be concerned about this warning? In what circumstances would that be an issue? Please feel free to come up with examples of your own.

Edit: Here's a way that I found that avoid the scoping warning:

Mathematica graphics

Is this a better way? I think the main reason of using Function in Wellin's case and defining a function within the Module here is to make a regular expression (which can't take in variable directly) become a function that could take in a value (xi in this function). What's the best way to do it? This is closely related to this question here.

Community
  • 1
seismatica
  • 5,101
  • 1
  • 22
  • 33
  • 5
    Simple errors and typos in books are best taken up with the author/publisher... – rm -rf Aug 12 '14 at 23:29
  • The thing is if it's really a typo (which I'm pretty sure it is), then what is he really trying to do? If he meant to use var instead of the typo fvar, then was he aware of the scoping warning. As a result, my question also asks if it's a good idea to construct the function that way, using Function. In other words, do you or other MMA experts on here consider that a good practice? I should have worded my title a little bit better. – seismatica Aug 12 '14 at 23:48
  • Yeah, the title and perhaps most of the question... To answer the question in your comment, yes, people do write code that way if necessary. The warning exists because more often than not, this is a mistake made by inexperienced users and is probably not what they want. If you know what you're doing, you're free to ignore the warning. There have been related questions such as this on syntax warnings. – rm -rf Aug 13 '14 at 00:05
  • my wag older mma versions didnt flag the issue with a warning (dont know how old that book is) – george2079 Aug 13 '14 at 00:47
  • 1
    For what it's worth, this seems to clear up the redness: fun = Function[Evaluate[var], expr]. – Daniel Lichtblau Aug 13 '14 at 15:55
  • There's something deeper wrong with the proffered code: the stopping criterion is that the magnitude of $f$ at the current approximation $x_i$ be sufficiently small. But such $x_i$ might be quite far from the desired root if the graph of $f$ is flat near there. One should substitute, or append, a rule that $|x_i - x_{i-1}|$ be sufficiently small. – murray Aug 14 '14 at 16:21

2 Answers2

7

"Is this a better way?"

Yes, but I think there are alternatives:

Nest

newton1[fun_, xi_, n_] := With[{f = fun/D[fun, x]}, Nest[# - f /. x -> # &, 2., n]]

newton1[x^3 - 2, 2., 10]

1.25992

NestList

newton2[fun_, xi_, n_] := With[{f = fun/D[fun, x]}, NestList[# - f /. x -> # &, 2., n]]

ListLinePlot[newton2[x^3 - 2, 2., 10], Mesh -> All, MeshStyle -> Directive[PointSize[Medium], Red]]

enter image description here

For the next functions we define

 f = # / D[#, x] & [x^3 - 2];

NestWhileList

NestWhileList[# - f /. x -> # &, 2., Abs[#1 - #2] > 0.0001 &, 2]

{2., 1.5, 1.2963, 1.26093, 1.25992, 1.25992}

NestWhile

NestWhile[# - f /. x -> # &, 2., Abs[#1 - #2] > 0.0001 &, 2]

1.25992

FixedPointList

FixedPointList[# - f /. x -> # &, 2.]

{2., 1.5, 1.2963, 1.26093, 1.25992, 1.25992, 1.25992, 1.25992}

FixedPoint

FixedPoint[# - f /. x -> # &, 2.]

1.25992

Fold (or FoldList)

Fold[# - f /. x -> # &, 2., Range@10]

1.25992

Your last solution could also be written as follows:

findRoot[expr_, {var_, init_}, e_] :=
 Module[{xi = init, fun},
  fun[x_] := expr / D[expr, var] /. var :> x;
  While[ Abs @ fun @ xi > e, xi -= fun @ xi ];
  {var -> xi}]

Here is a timetable computed over 1000 runs with initial value 0.1:

enter image description here

eldo
  • 67,911
  • 5
  • 60
  • 168
  • In timing results, why did you switch to initial approximation of 0.1 whereas the various methods preceding that all started at 2.0? – murray Aug 14 '14 at 17:52
  • @murray Because 0.1 is farer away from 1.25... than 2.0, and I wanted to see how the different methods would cope with this. – eldo Aug 14 '14 at 20:16
3

Here is one approach to using a more reasonable stopping criterion. For clarity we separately define a function that effects one step of the Newton-Raphson procedure. The first argument will be a function rather than a function expression.

  newtonStep[{f_, x_}] := {f, x - f[x]/f'[x]}

  newton[f_, start_, \[Delta]_] := 
   Last /@ NestWhileList[newtonStep, {f, start}, Abs[Last@(#1 - #2)] >= \[Delta] &, 2]

For example:

  f[x_] := x^3 - 2
  newton[f, 2., 10.^-12] // NumberForm[#, 15] &
(* {2., 1.5, 1.2962962962962963, 1.2609322247417485, 1.2599218605659261,   
    1.2599210498953948, 1.2599210498948732} *)

Of course one may in effect use a "function expression" as the first argument by using a pure function. For example:

  newton[#^3 - 2 &, 2., 10.^-12]
murray
  • 11,888
  • 2
  • 26
  • 50