Solving non-linear systems of equations
Non-linear systems of equations arise in many different applications. As for optimization, different situations will call for different inputs and optimizations: scalar, static, mutating, iterative solvers, etc. Below, we will cover some important algorithms and special cases.
Scalar solving with 1st order derivatives
Assume that you have a residual system that you want to solve. This means that if your original system to solve is is F(x) = K
for some real number K
, you have made sure to provide a residual function G(x) = F(x) - K
that we can solve for G(x) = 0
instead. Then, you can use Newton's method for root finding as follows:
function F(x)
x^2
end
function FJ(Jx, x)
x^2, 2x
end
prob_obj = NLSolvers.ScalarObjective(
f=F,
fg=FJ,
)
prob = NEqProblem(prob_obj; inplace = false)
x0 = 0.3
res = solve(prob, x0, LineSearch(Newton(), Backtracking()))
with output
julia> res = solve(prob, x0, LineSearch(Newton(), Backtracking()))
Results of solving non-linear equations
* Algorithm:
Newton's method with default linsolve with backtracking (no interp)
* Candidate solution:
Final residual 2-norm: 5.36e-09
Final residual Inf-norm: 5.36e-09
Initial residual 2-norm: 9.00e-02
Initial residual Inf-norm: 9.00e-02
* Stopping criteria
|F(x')| = 5.36e-09 <= 1.00e-08 (true)
* Work counters
Seconds run: 1.41e-05
Iterations: 12
The output reports initial and final residual norms. The convergence check is also reported, and some work counters report the time spent and the number of iterations.
Multivariate non-linear equation solving
Multivariate non-linear equation solving requires writing a VectorObjective
instead of a ScalarObjective
as above. The VectorObjective
is
using NLSolvers, ForwardDiff
function theta(x)
if x[1] > 0
return atan(x[2] / x[1]) / (2.0 * pi)
else
return (pi + atan(x[2] / x[1])) / (2.0 * pi)
end
end
function F_powell!(Fx, x)
if (x[1]^2 + x[2]^2 == 0)
dtdx1 = 0
dtdx2 = 0
else
dtdx1 = -x[2] / (2 * pi * (x[1]^2 + x[2]^2))
dtdx2 = x[1] / (2 * pi * (x[1]^2 + x[2]^2))
end
Fx[1] =
-2000.0 * (x[3] - 10.0 * theta(x)) * dtdx1 +
200.0 * (sqrt(x[1]^2 + x[2]^2) - 1) * x[1] / sqrt(x[1]^2 + x[2]^2)
Fx[2] =
-2000.0 * (x[3] - 10.0 * theta(x)) * dtdx2 +
200.0 * (sqrt(x[1]^2 + x[2]^2) - 1) * x[2] / sqrt(x[1]^2 + x[2]^2)
Fx[3] = 200.0 * (x[3] - 10.0 * theta(x)) + 2.0 * x[3]
Fx
end
function F_jacobian_powell!(Fx, Jx, x)
ForwardDiff.jacobian!(Jx, F_powell!, Fx, x)
Fx, Jx
end
prob_obj = VectorObjective(F=F_powell!, FJ=F_jacobian_powell!)
prob = NEqProblem(prob_obj)
x0 = [-1.0, 0.0, 0.0]
res = solve(prob, copy(x0), LineSearch(Newton(), Backtracking()))
with result
julia> res = solve(prob, copy(x0), LineSearch(Newton(), Backtracking()))
Results of solving non-linear equations
* Algorithm:
Newton's method with default linsolve with backtracking (no interp)
* Candidate solution:
Final residual 2-norm: 8.57e-16
Final residual Inf-norm: 6.24e-16
Initial residual 2-norm: 1.88e+03
Initial residual Inf-norm: 1.59e+03
* Stopping criteria
|F(x')| = 6.24e-16 <= 1.00e-08 (true)
* Work counters
Seconds run: 2.10e-04
Iterations: 33
We again see initial and final residual norms. This time the 2- and Inf-norm are different because there is more than one element in the state to be optimized over. The final 2-norm also satisfies the required threshold, and the run-time and number of iterations are reported.