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 18579 days, 8 hours, 40 minutes ago.