Solver Functions Support
The Fractal Science Kit fractal generator Solver functions allow you to produce Newton fractals and other root-finding method based fractals. The Fractal Science Kit fractal generator Solver functions support the following root-finding methods: Newton's method, Schroder's method, Halley's method, Whittaker's method, Cauchy's method, the Super-Newton method (also called Householder's method), the Super-Halley method, the Chebyshev-Halley family of methods, and the C-Iterative family of methods. Developers can set up a program to examine the fractals associated with these methods with just a few lines of code.
Solver.Apply(method, methodArg, z, f0, f1, f2)
Solver.ApplyToPolynomial(method, methodArg, coef[], z)
Solver.Apply applies the given method and methodArg to the value z and returns a value closer (hopefully) to a root of the associated function. The values of f0, f1, and f2, are the function evaluated at the value z and the 1st and 2nd derivative also evaluated at z.
Solver.ApplyToPolynomial is more efficient and easier to set up but can only be used for polynomial based equations. This function computes the 1st and 2nd derivatives of the polynomial function defined by the coefficients in the coef[] array, and applies the given root-finding method to the value z. The number of elements in the coef[] array defines the degree of the polynomial (e.g., a 5th degree polynomial would have 6 values) and any terms missing in the equation must have a 0 coefficient. The coefficients are ordered from high to low in the array; i.e., the coefficient associated with the term with the highest degree is in the array at index 0. The degree of the polynomial must be greater than 0. Solver.ApplyToPolynomial is highly optimized.
Example:
comment:
5th degree polynomial fractal equation based on Newton's method.
f(z) = z^5 - 5/3*(K+L+M)*z^4 + 10/3*(K*L+K*M+L*M)*z^3 -
10*K*L*M*z^2 + c
f'(z) = 5*z^4 - 20/3*(K+L+M)*z^3 + 10*(K*L+K*M+L*M)*z^2 -
20*K*L*M*z
f''(z) = 20*(z-K)*(z-L)*(z-M)
Critical points: K, L, M
Set the Classic Controller to Newton and use the
Julia preview window to explore the Julia fractals.
global:
FSK.OverrideValue("RootDetection", True)
FSK.OverrideValue("MaxPower", Solver.MaxPower(SolverMethod))
Complex coef[] = 1, -5/3*(K+L+M), 10/3*(K*L+K*M+L*M), -10*K*L*M,
0, c
initialize:
if (~IsJulia) {
coef[5] = c
switch (criticalPoint) {
case CriticalPoints.K: z = K
case CriticalPoints.L: z = L
case CriticalPoints.M: z = M
}
}
iterate:
z = Solver.ApplyToPolynomial(SolverMethod, SolverMethodArg,
coef[], z)
properties:
#include SolverMethodOptions
option K {
type = Complex
caption = "K"
default = -0.5
}
option L {
type = Complex
caption = "L"
default = 0.5i
}
option M {
type = Complex
caption = "M"
default = -0.5i
}
enum CriticalPoints {
K, "K"
L, "L"
M, "M"
}
option criticalPoint {
type = criticalPoints
caption = "Critical Point"
default = CriticalPoints.K
}
This Fractal Equation defines the coefficients of the base equation in an array called coef[] and simply calls Solver.ApplyToPolynomial to apply the root-finding method SolverMethod to the value z.
In this example, the degree of the polynomial is 5 since the number of elements in the array is 6, and the coefficient of the next to last term (i.e., the coefficient for the z term) is 0.
Notice how we set coef[5]=c if we are processing a Mandelbrot fractal. This is required for Mandelbrot fractals since c is initialized to the pixel location just prior to calling the initialize section so the assignment in the global section is not correct with respect to coef[5]. This is not required for Julia fractals since c is set to the Julia Constant before the program runs and is never changed. We could move the initialization of the coef[] array out of the global section and into the initialize section to avoid this complexity and thereby simplify the code but this would be less efficient since the initialize section is executed at the beginning of every iteration while the global section is executed once during program compilation. In fact, since IsJulia is a constant, the compiler will evaluate the if statement during Program Optimization and remove the entire if block for Julia fractals, and at runtime, the initialize section is not even executed!
SolverMethod is a value from the SolverMethods enum defined in the built-in macros:
#define SolverMethods
enum SolverMethods {
Newton,
"Newton's method"
Schroder, "Schroder's
method"
Halley,
"Halley's method"
Whittaker, "Whittaker's method"
Cauchy,
"Cauchy's method"
SuperNewton, "Super-Newton method"
SuperHalley, "Super-Halley method"
ChebyshevHalley, "Chebyshev-Halley family"
CIterative, "C-Iterative family"
}
#end
The SolverMethodArg is used to specify an argument to 2 of the root-finding methods. For the ChebyshevHalley, the argument is the family parameter Theta and for the CIterative, the argument is the family parameter C. The ChebyshevHalley parameter Theta is related to several of the other methods as follows:
Theta = 0 ->
Super-Newton Method
Theta = 0.5 -> Halley's Method
Theta = 1 -> Super-Halley Method
Theta = Infinity -> Newton's Method (use 1e8 to approximate Infinity)
The SolverMethod and SolverMethodArg options are defined by including the SolverMethodOptions in the properties section using the statement:
#include SolverMethodOptions
SolverMethodOptions is defined in the built-in macros but is given here for reference:
#define SolverMethodOptions
#include SolverMethods
option SolverMethod {
type = SolverMethods
caption = "Method"
details = "Root finding method"
default = SolverMethods.Newton
}
option SolverMethodArg {
type = Float
caption = "Argument"
details = "Theta (Chebyshev-Halley) or C (C-Iterative)"
default = 1
enabled = SolverMethod = SolverMethods.ChebyshevHalley || \
SolverMethod
= SolverMethods.CIterative
}
#end
Including the SolverMethodOptions in the properties section adds the following properties to the page:
The different root-finding methods each look best with different settings for the Max Power property and you can use the function Solver.MaxPower to override the value in the global section using the statement:
FSK.OverrideValue("MaxPower", Solver.MaxPower(SolverMethod))
Example:
comment:
Mandelbrot fractal based on Newton's method for
finding roots applied to:
Sin(z) - c = 0
Set the Classic Controller to 'Gradient Map - Newton' and
use the Julia preview window to explore the Julia fractals.
Notes:
f(z) = Sin(z) - c
f'(z) = Cos(z)
f''(z) = -Sin(z)
global:
FSK.OverrideValue("RootDetection", True)
FSK.OverrideValue("MaxPower", Solver.MaxPower(SolverMethod))
iterate:
sz = Sin(z)
cz = Cos(z)
f0 = sz - c
f1 = cz
f2 = -sz
z = Solver.Apply(SolverMethod, SolverMethodArg, z, f0, f1, f2)
properties:
#include SolverMethodOptions
This example uses Solver.Apply to implement the Fractal Equation. Like Solver.ApplyToPolynomial, the 1st and 2nd arguments are user specified options that define the method and argument for the fractal. However, Solver.Apply also requires you to evaluate the function and the function's 1st and 2nd derivative at the current point z and pass these 3 values as the last 3 arguments as illustrated above.
For details on Newton's method see:
From Wikipedia, the free encyclopedia.
Eric W. Weisstein. "Newton's Method."
From MathWorld--A Wolfram Web Resource.
For details on many of the other algorithms, I refer you to the paper:
"Review of some iterative root-finding methods from a dynamical point of view"
by Sergio Amat, Sonia Busquier, and Sergio Plaza.