Nonlinear Systems in IronPython QuickStart Sample
Illustrates the use of the NewtonRaphsonSystemSolver and DoglegSystemSolver classes for solving systems of nonlinear equations in IronPython.
View this sample in: C# Visual Basic F#
```Python
import numerics
from math import exp, sin, cos
from System import Array, Func
# The equation solver classes reside in the
# Extreme.Mathematics.EquationSolvers namespace.
from Extreme.Mathematics.EquationSolvers import *
# Function delegates reside in the Extreme.Mathematics
# namespace.
from Extreme.Mathematics import *
# Vectors reside in the Extreme.Mathematics.LinearAlgebra
# namespace.
from Extreme.Mathematics.LinearAlgebra import *
# Illustrates solving systems of non-linear equations using
# classes in the Extreme.Mathematics.EquationSolvers namespace
# of Numerics.NET.
#
# Target function
#
# The function we are trying to solve can be provided
# on one of two ways. The first is as an array of
# Func<Vector, double> delegates. See the end of this
# sample for definitions of the methods that are referenced here.
f = Array[Func[Vector, float]]([ \
lambda x: exp(x[0])*cos(x[1]) - x[0]*x[0] + x[1]*x[1], \
lambda x: exp(x[0])*sin(x[1]) - 2*x[0]*x[1] ])
# We can also supply the Jacobian, which is the matrix of partial
# derivatives. We do so by providing the gradient of each target
# function as a delegate that takes a second argument:
# the vector that is to hold the return value. This avoids unnecessary
# creation of new Vector instances.
def df1(x,y):
y[0] = exp(x[0])*cos(x[1]) - 2*x[0]
y[1] = -exp(x[0])*sin(x[1]) + 2*x[1]
return y
def df2(x,y):
y[0] = exp(x[0])*sin(x[1]) - 2*x[1]
y[1] = exp(x[0])*cos(x[1]) - 2*x[0]
return y
df = Array[Func[Vector,Vector,Vector]]([ df1, df2 ])
# The initial values are supplied as a vector:
initialGuess = Vector([0.5, 0.5])
#
# Newton-Raphson Method
#
# The Newton-Raphson method is implemented by
# the NewtonRaphsonSystemSolver class.
solver = NewtonRaphsonSystemSolver(f, df, initialGuess)
# and call the Solve method to obtain the solution:
solution = solver.Solve()
print "N-dimensional Newton-Raphson Solver:"
print "exp(x)*cos(y) - x^2 + y^2 = 0"
print "exp(x)*sin(y) - 2xy = 0"
print " Initial guess: {0:F2}", initialGuess
# The Status property indicates
# the result of running the algorithm.
print " Status:", solver.Status
# The result is also available through the
# Result property.
print " Solution:", solver.Result
print " Function value:", solver.ValueTest.Error
# You can find out the estimated error of the result
# through the EstimatedError property:
print " Estimated error:", solver.EstimatedError
print " # iterations:", solver.IterationsNeeded
print " # evaluations:", solver.EvaluationsNeeded
### Symbolic derivatives are currently not supported ###
## In .NET 4.0, you can use Automatic Differentiation to compute the Jacobian.
## To do so, set the target functions using the SetSymbolicTargetFunctions
## method:
#solver = NewtonRaphsonSystemSolver()
#solver.InitialGuess = initialGuess
#
#solver.SetSymbolicTargetFunctions( \
# lambda x: exp(x[0]) * cos(x[1]) - x[0] * x[0] + x[1] * x[1], \
# lambda x: exp(x[0]) * sin(x[1]) - 2 * x[0] * x[1])
#
#solution = solver.Solve()
#print "Using Automatic Differentiation:"
#print " Solution:", solver.Result
#print " Function value:", solver.ValueTest.Error
#print " # iterations:", solver.IterationsNeeded
# When you don't have the Jacobian for the target functions
# and you don't use Automatic Differentiation, the equation solver
# will use a numerical approximation.
#
# Controlling the process
#
print "Same with modified parameters:"
# You can set the maximum # of iterations:
# If the solution cannot be found in time, the
# Status will return a value of
# IterationStatus.IterationLimitExceeded
solver.MaxIterations = 10
# The ValueTest property returns the convergence
# test based on the function value. We can set
# its tolerance property:
solver.ValueTest.Tolerance = 1e-10
# Its Norm property determines how the error
# is calculated. Here, we choose the maximum
# of the function values:
solver.ValueTest.Norm = Algorithms.VectorConvergenceNorm.Maximum
# The SolutionTest property returns the test
# on the change in location of the solution.
solver.SolutionTest.Tolerance = 1e-8
# You can specify how convergence is to be tested
# through the ConvergenceCriterion property:
solver.SolutionTest.ConvergenceCriterion = ConvergenceCriterion.WithinRelativeTolerance
solver.InitialGuess = initialGuess
solution = solver.Solve()
print " Status:", solver.Status
print " Solution:", solver.Result
# The estimated error will be less than 5e-14
print " Estimated error:", solver.SolutionTest.Error
print " # iterations:", solver.IterationsNeeded
print " # evaluations:", solver.EvaluationsNeeded
#
# Powell's dogleg method
#
# The dogleg method is more robust than Newton's method.
# It will converge often when Newton's method fails.
dogleg = DoglegSystemSolver(f, df, initialGuess)
# Unique to the dogleg method is the TrustRegionRadius property.
# Any step of the algorithm is not larger than this value.
# It is adjusted at each iteration.
dogleg.TrustRegionRadius = 0.5
# Call the Solve method to obtain the solution:
solution = dogleg.Solve()
print "Powell's Dogleg Solver:"
print " Initial guess: {0:F2}", initialGuess
print " Status:", dogleg.Status
print " Solution:", dogleg.Result
print " Estimated error:", dogleg.EstimatedError
print " # iterations:", dogleg.IterationsNeeded
print " # evaluations:", dogleg.EvaluationsNeeded
# The dogleg method can work without derivatives. Care is taken
# to keep the number of evaluations down to a minimum.
dogleg.JacobianFunction = None
# Call the Solve method to obtain the solution:
solution = dogleg.Solve()
print "Powell's Dogleg Solver (no derivatives):"
print " Initial guess: {0:F2}", initialGuess
print " Status:", dogleg.Status
print " Solution:", dogleg.Result
print " Estimated error:", dogleg.EstimatedError
print " # iterations:", dogleg.IterationsNeeded
print " # evaluations:", dogleg.EvaluationsNeeded
```