secant.frink

View or download secant.frink in plain text format



// 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]
   while (true)
   {
      x = x1 - (f1 * (x1 - x2)) / (f1-f2)
      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
      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]]


View or download secant.frink in plain text format

This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 18864 days, 15 hours, 53 minutes ago.