2
$\begingroup$

I found this command in Mathematica documentation to do various types of root finding.

 <<Optimization`UnconstrainedProblems`
 FindRootPlot[{Tan[Pi*x] - 6}, {x, 0, 0.48}]

It output this:

  {{x -> 0.447432}, {"Steps" -> 10, "Residual" -> 11}

I'm wondering if there's a way to get it to show output for all the iterations of the recursion?

$\endgroup$
  • $\begingroup$ Have you seen StepMonitor? $\endgroup$ – Michael E2 Sep 11 '16 at 23:47
  • 1
    $\begingroup$ When I run your example, I get a plot of the steps already, without having to do anything. Don't you get it, too? $\endgroup$ – Michael E2 Sep 11 '16 at 23:48
  • 1
    $\begingroup$ @MichaelE2: I also see what you see. I think he had something more along the lines of your first comment where he actually wants to see a table of the calculations for the Secant method including the error if possible. I wish this would be made easier to be quite honest given that that is how every numerical analysis text shows it. Maybe there is something to that there is too much generalization when it hides the incredible things MMA can do from the average user. $\endgroup$ – Moo Sep 12 '16 at 0:08
  • $\begingroup$ Yes @MichaelE2 I did. I was hoping for a table of numerical values. How do I use StepMonitor? Thank you! $\endgroup$ – Jabernet Sep 12 '16 at 0:27
  • 2
    $\begingroup$ @Moo The docs for FindRoot show how to use the option StepMonitor to collect the steps. For what you suggest, it would be changed to StepMonitor :> Sow[{x, Tan[Pi*x] - 6}]. A table could be formatted with TableForm, Grid or TeXForm. Not sure why there isn't a StepTable function. It doesn't seem nearly as useful as the plot, though, except for teaching or learning about convergence rates, which is a fairly limited application. $\endgroup$ – Michael E2 Sep 12 '16 at 1:40
4
$\begingroup$

Introduction

Here's a utility I might use for generating tables for class, modeled on the NDSolve Method Plugin Framework. It collects the steps and function values with StepMonitor. To this data is applied a post-processor to do with it what it will. The built-in default constructs a table. The processor can be passed various data arguments represented by the left-hand-side strings found in the rules in $stepDataElements and $stepDataInject[] below. In the implementation further down, it is assumed that the data arguments are bundled in a list in first position, followed optionally by further arguments, like this:

processor[{"X", "NF"},..]

This could be made more general with some work, but it suited the applications I had in mind. Note in this setup, the values for each step are found in the columns of the array. Depending on the application, transposing the array may be necessary.

Code

Here are the possible arguments. The rules $stepDataElements give label names for the elements and $stepDataInject[data] injects their values. More elements could be added.

ClearAll[$stepDataInject];
$stepDataElements = {        (* arguments for processor:    *)
   "N" -> "Step",            (*   step number               *)
   "X" -> "X",               (*   input vector X            *)
   "F" -> "Function",        (*   output value/vector F[X]  *)
   "NF" -> "|Function|",     (*   norm of output vector     *)
   "DX" -> "ΔX",             (*   step change X[n]-X[n-1]   *)
   "NDX" -> "|ΔX|",          (*   norm of step (step size)  *)
   "RF" -> "Ratio"};         (*   ratio of norm of output NF[X[n]]/NF[X[n-1]]  *) 
$stepDataInject[data_] := {  (* computations of elements    *)
   "N" :> Range@Length@data,
   "X" :> data[[All, 1]],
   "F" :> data[[All, 2]],
   "NF" :> stepData`norm /@ data[[All, 2]],
   "DX" :> Prepend[Differences@data[[All, 1]], ""],           (* not numeric arrays! *)
   "NDX" :> Prepend[stepData`norm /@ Differences@data[[All, 1]], ""],
   "RF" :> Prepend[Ratios[stepData`norm /@ data[[All, 2]]], ""]};

In addition to the data arguments, the strings "Arguments" and "Elements" may appear in other arguments to the post-processor; they will be replaced by data argument strings and their labels respectively.

Here is stepData function. There are two built-in post-processors, one to generate a table and one just to return the data with the steps in the rows. They are used like this:

StepMonitor -> "Table"[{elem1, elem2,...}, opts]
StepMonitor -> "Data"[{elem1, elem2,...}]

where elem1, elem2 are chosen from the elements "N", "X", "F", etc. The default processor is

StepMonitor -> "Table"[{"X", "NF"}, TableHeadings -> {Automatic, "Elements"}],

which gives a table of the input step and norm of the residual, with numbered rows and columns labeled by the elements in the table.

ClearAll[stepData];
SetAttributes[stepData, HoldAll];
Options[stepData] = {
   StepMonitor -> "Table"[{"X", "NF"}, TableHeadings -> {Automatic, "Elements"}],
   NormFunction -> Norm};

(* processor[data_, args___] -- data must be an array *)
$stepDataProcessors = {
   "Data" -> Transpose, 
   "Table" -> (TableForm[Transpose[#1], ##2] &)};

stepData[FindRoot[f_, vars_, fropts___], opts : OptionsPattern[]] :=
   Module[{ff, monitor, tag, getVars, X, sol, data, res},
   Block[{stepData`norm = OptionValue[NormFunction]},
    getVars[{x_, __?NumericQ}] := x;
    getVars[v : {__List}] := v[[All, 1]];
    monitor = OptionValue[Automatic, Automatic, StepMonitor, Hold];
    X = getVars[vars];
    ff = f /. Equal -> Subtract;

    (* How to avoid double evaluation f?? *)
    res = Reap[
      FindRoot[ff, vars, StepMonitor :> Sow[{X, ff}, tag], fropts],
      tag
      ];
    (data = res[[2, 1]];
      monitor /. 
        Hold[h_[d_, o___]] :> 
         h[d /. $stepDataInject[data], 
          Hold[o] /. {"Arguments" -> d, 
             "Elements" -> d /. $stepDataElements} // 
           ReleaseHold] /. $stepDataProcessors) /; FreeQ[res, FindRoot]
    ]];

Examples

stepData[FindRoot[{Tan[Pi*x] - 6}, {x, 0, 0.48}]]

Mathematica graphics

Modified "Table" post-processor:

stepData[FindRoot[Exp[x] == 2, {x, 0.3, 0.4}],
 StepMonitor -> 
  "Table"[{"X", "NDX", "NF"}, TableHeadings -> {Automatic, "Elements"}]]

Mathematica graphics

Just the raw data:

stepData[FindRoot[Exp[x] == 2, {x, 0.3, 0.4}],
 StepMonitor -> "Data"[{"N", "X", "NF"}]]
(*
  {{1, 0.757956, 0.13391},
   {2, 0.683303, 0.0195927},
   {3, 0.692831, 0.000632041},
   {4, 0.693149, 3.11688*10^-6},
   {5, 0.693147, 4.92602*10^-10},
   {6, 0.693147, 4.44089*10^-16}}
*)

Multivariate problem:

stepData[FindRoot[{Tan[Pi*x] - 6 + y, y^2 - x}, {{x, 0.3, 0.4}, {y, 0.5, 0.8}}],
 StepMonitor -> 
  "Table"[{"X", "NDX", "NF"}, TableHeadings -> {Automatic, "Elements"}]]

Mathematica graphics

Custom postprocessor, comparing step size with the residual:

ClearAll[myFmt];
myFmt[data_?MatrixQ, args_, opts_] := 
  ListLogPlot[data, opts, Frame -> True, 
   FrameLabel -> {"Step", 
     Row[args /. {"NF" -> "Residual", "NDX" -> "Step error"}, 
      ", "]}];
stepData[FindRoot[Exp[x] == 2, {x, 3, 4}, WorkingPrecision -> 50],
 StepMonitor -> myFmt[{"NDX", "NF"}, "Arguments", PlotRange -> All]]

Mathematica graphics

Plotting one versus the other (note the gymnastics needed to make the data the first argument):

stepData[FindRoot[x^2 == 2, {x, 3, 4}, WorkingPrecision -> 50],
 StepMonitor :> (ListLogLogPlot@*Transpose)[{"NDX", "NF"}]]

Mathematica graphics

$\endgroup$
  • $\begingroup$ a very interesting utility, +1. $\endgroup$ – rcollyer Sep 13 '16 at 13:47

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.