Download or view 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]
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])]
*/
Download or view 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 19305 days, 18 hours, 5 minutes ago.