# secant.frink

``` // Secant method solver. // // This function finds a root of the function f. //  In other words, it returns the value of x for which f[x] = 0. // // The arguments are: //   f: A function that takes a single argument //   x1, x2:  initial guesses //   maxDelta:  the maximum error in x secant[f, x1, x2, maxDelta = 1e-14] := {    f1 = f[x1]    f2 = f[x2]    x = undef    while (true)    {       diff = f1 - f2       if diff == 0          return x1              x = x1 - (f1 * (x1 - x2)) / diff       // println[x]       if abs[x - x1] < maxDelta          return x       x2 = x1       x1 = x       f2 = f1       f1 = f[x]    } } // This uses the secant method to invert the function y = f[x]. // This will essentially find an inverse function for f[x] and return a value //  of x for which f[x] = y. //  other parameters: //   x1,x2:  initial guesses that hopefully bound the desired result. //   maxDelta:  maximum error in y //   TODO:  Use interval techniques to make this more rigorous and powerful? //   TODO:  Automatically make guesses for x1 and x2?  Somehow? secantInvert[f, y, xmin, xmax, maxDelta = 1e-14] := {    x1 = xmin    x2 = xmax    y1 = f[x1]    y2 = f[x2]    xnew = (x2-x1)/2 + x1    while true    {       ydiff = y2 - y1 //      println["ydiff is \$ydiff, x1 is \$x1, x2 is \$x2"]       if ydiff == 0 y   // Degenerate case to avoid dividing by zero.          return xnew    // This may not be always a correct solution?              invSlope = (x2-x1) / ydiff       xnew = x1 + (y - y1) invSlope       if xnew < xmin          xnew = xmin       if xnew > xmax          xnew = xmax              ynew = f[xnew] //      println["xnew=\$xnew\tynew=\$ynew"]       if ynew == 0 y    // Degenerate case to avoid dividing by zero.          return xnew      // This may not be always a correct solution?              if abs[(ynew - y) / ynew] < maxDelta          return xnew       y2 = y1       y1 = ynew       x2 = x1       x1 = xnew    } } // This uses the secant method to invert the function y = f[x], assuming that // y is an angle.  This prevents some overcorrections when angles are // negative. // This will essentially find an inverse function for f[x] and return a value //  of x for which f[x] = y. //  other parameters: //   x1,x2:  initial guesses that hopefully bound the desired result. //   maxDelta:  maximum error in y //   TODO:  Use interval techniques to make this more rigorous and powerful? //   TODO:  Automatically make guesses for x1 and x2?  Somehow? secantInvertAngle[f, y, xmin, xmax, maxDelta = 0.015 arcsec] := {    x1 = xmin    x2 = xmax    y1 = f[x1]    y2 = f[x2]    xnew = (x2-x1)/2 + x1    while true    {       y1e = y1 - y       y2e = y2 - y       if y1e > 180 deg          y1 = y1 - circle       if y2e > 180 deg          y2 = y2 - circle       if y1e < -180 deg          y1 = y1 + circle       if y2e < -180 deg          y2 = y2 + circle              ydiff = y2 - y1              // did we wrap around the circle?       if abs[y1 - y2] > 180 degrees          if (y1 < y2)             y2 = y2 - circle          else             y1 = y1 - circle //      println["ydiff is \$ydiff, x1 is \$x1, x2 is \$x2, y1 is \$y1, y2 is \$y2"]       if ydiff == 0 y   // Degenerate case to avoid dividing by zero.          return xnew    // This may not be always a correct solution?              invSlope = (x2-x1) / ydiff       xnew = x1 + (y - y1) invSlope       if xnew < xmin          xnew = xmin       if xnew > xmax          xnew = xmax              ynew = f[xnew] //      println["xnew=\$xnew\tynew=\$ynew"]       if ynew == 0 y    // Degenerate case to avoid dividing by zero.          return xnew      // This may not be always a correct solution?              if abs[(ynew - y) / ynew] < maxDelta          return xnew       y2 = y1       y1 = ynew       x2 = x1       x1 = xnew    } } // Minimize a function using the secant method.  This doesn't really work yet. secantMinimize[f, xmin, xmax, minStepX] := {    x1 = xmin    x2 = xmax    y1 = f[x1]    y2 = f[x2]    while true    {       println["x1=\$x1\t x2=\$x2"]       diff = x2-x1       if diff == 0          return f[x1]              slope = (y2-y1)/diff       xnew = x1 + slope (x1+x2)/2       ynew = f[xnew]       println["ynew=\$ynew\txnew= \$xnew"]       if (abs[x2-x1] < minStepX)          return ynew              y2 = y1       y1 = ynew       x2 = x1       x1 = xnew       if x1 > x2          [x1, x2] = [x2, x1]       if x1 < xmin          x1 = xmin       if x2 > xmax          x2 = xmax    } } // Sample root-finding: //   Define a procedure block that represents the equation // (this is just a function without a name, or think of it //  as a reference to a function.) //f = { |x| ln[x] - 1} //println["Solution: " + secant[f, 1, 3]] // Sample inverse-finding: //   Find an inverse for the following function. //   The call below finds a value x such that log[x]=2 //     in other words, calculates 10^2 // f = { |x| log[x] } // println[secantInvert[f, 2, 1, 200, 1e-20]] // Example secant method using arbitrary precision to calculate sin[x] /* use ArbitraryPrecision.frink use pi2.frink digits = 75 setPrecision[digits] f ={|x| arbitrarySin[x, 75]} y = secant[f, 3, 4, 10^-digits] println["Solution is \$y"] println["Difference from pi is " + (y - Pi.getPi[digits])] */```