3
$\begingroup$

This is an example of procedural programming borrowed from Paul Wellin's book. It's a simple implementation of the Newton-Rhapson method.

I overloaded the findRoot function to accommodate different kinds of arguments. It works fine in the following first case:

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

f[x_] := x^2 - 2;

findRoot[f,{x,2},.0001]
{x-> 1.41422}

But when I overload it to accommodate expressions, like the following for example, it doesn't evaluate.

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

findRoot[x^2 - 2 == 0, {x, 2.0}, 0.0001]
{x - >2.}

How do I make Mathematica evaluate it numerically all the way through?

$\endgroup$
  • 2
    $\begingroup$ Try fun = Function[fvar, Evaluate[expr - val]] -- This technique used here and here $\endgroup$ – Michael E2 Mar 25 at 2:25
  • $\begingroup$ @Michael E2 Still doesn't work but thanks. $\endgroup$ – Thadeu Freitas Filho Mar 25 at 4:49
  • $\begingroup$ I think there's a typo: var should replace fvar in your/my codes. Also Nest (or FixedPoint) could be used instead of While, if you haven't seen it. $\endgroup$ – Michael E2 Mar 25 at 13:33
4
$\begingroup$

It's a little tricky because Mathematica renames the variable given in the 1st argument of Function to localize it. So the renaming has to be taken into account when substituting the desired variable into Function. The unevaluated locally defined fun is then passed to the previously defined version of findRoot and all goes well.

findRoot[expr_ == val_, {var_, init_}, ϵ_] :=
  Module[{fun},
    fun = Function[u, expr - val] /. u$ -> var;
    findRoot[Unevaluated[fun], {var, init}, ϵ]]

Note: I pass Unevaluated[fun] to the previously defined version of findRoot because that function requires it's 1st argument to be a symbol.

findRoot[x^2 - 2 == 0, {x, 2.0}, 0.0001]

{x -> 1.41422}

Another useful test is to write the equation as x^2 == 2.

findRoot[x^2 == 2, {x, 2.0}, 0.0001]

{x -> 1.41422}

$\endgroup$
  • $\begingroup$ What's the purpose of the $ in u$? $\endgroup$ – Thadeu Freitas Filho Mar 26 at 4:36
  • 1
    $\begingroup$ @ThadeuFreitasFilho. It's there because, behind the scenes, Mathematica changes the u in Function[u, expr - val] to Function[u$, expr - val], so that is what the replacement has to match. That is what I was discussing in the 1st paragraph of the answer. $\endgroup$ – m_goldberg Mar 26 at 4:51
  • $\begingroup$ @ThadeuFreitasFilho. I have updated my answer with an explanation of why I pass Unevaluated[fun] to the original findRoot. $\endgroup$ – m_goldberg Mar 26 at 5:37
2
$\begingroup$

From my comment, with the OP's typo, which I did not notice, fixed:

findRoot[expr_ == val_, {var_, init_}, \[Epsilon]_] := 
 Module[{xi = init,
   fun = Function[var, Evaluate[expr - val]]},
  While[Abs[fun[xi]] > \[Epsilon],
   xi = N[xi - fun[xi]/fun'[xi]]];
  {var -> xi}]

findRoot[x^2 == 2, {x, 1.}, 10^-8]
(*  {x -> 1.41421}  *)

If you want the variable var to be protected from evaluation inside findRoot, just as it is in FindRoot, you can do the following:

ClearAll[findRoot];
SetAttributes[findRoot, HoldAll];
findRoot[expr_ == val_, {var_, init_}, \[Epsilon]_] := 
 Module[{xi = init,
   fun = Function @@ Hold[var, expr - val]},
  While[Abs[fun[xi]] > \[Epsilon],
   xi = N[xi - fun[xi]/fun'[xi]]];
  {var -> xi}];

Then the following finds the root, just like FindRoot, but just like FindRoot, x is evaluated after the solution is returned. This could be prevented by returning {HoldPattern[var] -> xi} instead of {var -> xi}, but FindRoot does not do it.

x = 2;
findRoot[x^2 == 2, {x, 1.}, 10^-8]
(*  {2 -> 1.41421}  *)

As @m_goldberg's answer shows, it's good to have one core routine and have all the interfaces call it. Here's another way but putting the core routine in an "internal" function:

ClearAll[findRoot, iFindRoot];
iFindRoot[fun_, dfun_, init_, \[Epsilon]_] := Module[{xi = init},
   While[Abs[fun[xi]] > \[Epsilon],
    xi = N[xi - fun[xi]/dfun[xi]]];
   xi];
findRoot[fun : _Symbol | _Function, {var_, 
    init_?NumericQ}, \[Epsilon]_] :=
  {var -> 
    iFindRoot[fun, fun', init, \[Epsilon]]};
findRoot[fun : _Symbol | _Function, {init_?
     NumericQ}, \[Epsilon]_] :=
  {iFindRoot[fun, fun', 
    init, \[Epsilon]]};
findRoot[expr_ == val_, {var_, init_?NumericQ}, \[Epsilon]_] := 
  Module[{fun},
   fun = Function[var, Evaluate[expr - val]];
   {var -> iFindRoot[fun, fun', init, \[Epsilon]]}];
$\endgroup$
  • $\begingroup$ I noticed that in both answers, your function starting with ClearAll[findRoot]; and m_goldberg's, you preclude fun from evaluating. However, seems to me for different reasons (I could be wrong!). m_goldberg uses Unevaluated to pass an unevaluated function to the first defined version of findRoot. You use HoldAll and Hold to prevent var from evaluating inside findRoot. Why is this prevention important? $\endgroup$ – Thadeu Freitas Filho Mar 26 at 4:27
  • 1
    $\begingroup$ @ThadeuFreitasFilho Some Mathematica functions effectively localize their argument variables so that the remain symbols. For instance, x = 2; FindRoot[x^2 == 2, {x, 1.}] will find the right root even though x^2 == 2 would evaluate to False if evaluated when x equaled 2. The use of Hold in the 2nd alternative definition prevents x and the sides of the equation from evaluating, using Function to localize the argument var. It's also important the findRoot has the attribute HoldAll. See also reference.wolfram.com/language/tutorial/… $\endgroup$ – Michael E2 Mar 26 at 10:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.