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">ArbitraryPrecision.frink
use pi2.frink">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 19965 days, 18 hours, 32 minutes ago.