Root finding Secant method

Question

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?

Answer 1

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}]]

Modified "Table" post-processor:

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

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"}]]

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]]

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"}]]