Skip to content

Commit 4d2e96a

Browse files
GerhardHippmannGerhardHippmann
committed
Enable jacobian computation by finite difference approximation (FiniteDiff)
1 parent a2a6798 commit 4d2e96a

File tree

3 files changed

+50
-21
lines changed

3 files changed

+50
-21
lines changed

‎Project.toml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@ version = "0.11.0"
66
[deps]
77
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
88
DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b"
9+
FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
910
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
1011
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
1112
Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"

‎src/NonlinearEquations.jl‎

Lines changed: 20 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@ export mild, high, extreme, quiet, info, verbose, debug, solveNonlinearEquations
1010
using Printf
1111
using LinearAlgebra
1212
import ForwardDiff
13+
import FiniteDiff
1314

1415
@enum Nonlinearity begin
1516
mild
@@ -45,13 +46,14 @@ end
4546
xtol = 1e-6,
4647
maxiter = 50,
4748
quasi = true,
49+
forwardjac = false,
4850
loglevel::Loglevel = info) -> convergence::Bool
4951
5052
Solve nonlinear equation system ``F(x) = 0`` by global Gauss-Newton method with error oriented convergence criterion and adaptive trust region strategy. Optionally Broyden's 'good' Jacobian rank-1 updates are used. In case of underdetermined systems of equations, the 'shortest' least square solution is computed.
5153
5254
`F!` is a C1 continuous function with length ``n`` of input and ``m`` of output vector where ``n >= m`` (determined or underdetermined system of equations). It has to be defined by
5355
F!(F::Vector, x::Vector) -> success::Bool
54-
returning true in case of successful evaluation. Jacobian is computed by [ForwardDiff](https://github.com/JuliaDiff/ForwardDiff.jl) package.
56+
returning true in case of successful evaluation.
5557
5658
`m` is the length of the output vector of `F!`.
5759
@@ -69,6 +71,8 @@ returning true in case of successful evaluation. Jacobian is computed by [Forwar
6971
7072
`quasi` enables switch to local quasi Newton iteration (which avoids costly Jacobian evaluations) in case of advanced convergence.
7173
74+
`forwardjac` enables analytic Jacobian evaluation using [ForwardDiff](https://github.com/JuliaDiff/ForwardDiff.jl) package. Otherwise numeric computation by [FiniteDiff](https://github.com/JuliaDiff/FiniteDiff.jl) is used.
75+
7276
`loglevel` is the log message level. Possible @enum values are `quiet`, `info`, `verbose` and `debug`.
7377
7478
Reference: P. Deuflhard: Newton Methods for Nonlinear Problems. - Affine Invariance and Adaptive Algorithms (section 4.4). Series Computational Mathematics 35, Springer (2004).
@@ -82,6 +86,7 @@ function solveNonlinearEquations!(F!::Function,
8286
xtol = 1e-6,
8387
maxiter = 50,
8488
quasi = true,
89+
forwardjac = false,
8590
loglevel::Loglevel = info)
8691

8792
n = length(x)
@@ -107,7 +112,14 @@ function solveNonlinearEquations!(F!::Function,
107112
dxbar = similar(x, n)
108113
fxk = similar(x, m)
109114
jac = similar(x, m, n)
110-
jcfg = ForwardDiff.JacobianConfig(F!, fxk, x)
115+
if forwardjac
116+
jcfg = ForwardDiff.JacobianConfig(F!, fxk, x)
117+
else
118+
xche = similar(x, n)
119+
fxche1 = similar(x, m)
120+
fxche2 = similar(x, m)
121+
jche = FiniteDiff.JacobianCache(xche, fxche1, fxche2)
122+
end
111123
dxl = similar(x, n)
112124
Pdxbar = similar(x, n)
113125
fxl = similar(x, m)
@@ -197,8 +209,11 @@ function solveNonlinearEquations!(F!::Function,
197209
end
198210

199211
# 1. Step k: Evaluate Jacobian matrix F'(x^k)
200-
# In nleq_err.c fxk is not updated -> inconsistent after return from quasi Gauss-Newton iteration
201-
ForwardDiff.jacobian!(jac, F!, fxk, x, jcfg)
212+
if forwardjac
213+
ForwardDiff.jacobian!(jac, F!, fxk, x, jcfg)
214+
else
215+
FiniteDiff.finite_difference_jacobian!(jac, F!, x, jche)
216+
end
202217
njac = njac + 1
203218

204219
# Scale Jacobian and change sign of it
@@ -321,7 +336,7 @@ function solveNonlinearEquations!(F!::Function,
321336
# Natural/restricted monotonicity test failed
322337
lambda_new = min(mue, 0.5*lambda) # Corrected damping factor (3.48)
323338
if lambda <= lambda_min
324-
lambda = lambda_new # Does not make sense, bug in nleq_err.c?
339+
lambda = lambda_new # Does not make sense, bug?
325340
else
326341
lambda = max(lambda_new, lambda_min)
327342
end

‎test/TestNonlinearEquations.jl‎

Lines changed: 29 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,7 @@ println("... Test NonlinearEquations.jl")
4343
nonlinearity = mild
4444
restricted = false
4545
quasi = true
46+
forwardjac = false
4647

4748
if loglevel >= info
4849
println("Linear problem:")
@@ -55,7 +56,7 @@ println("... Test NonlinearEquations.jl")
5556
println("* scale vector = $fscale")
5657
end
5758

58-
converged = solveNonlinearEquations!(fx!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, loglevel=loglevel)
59+
converged = solveNonlinearEquations!(fx!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, forwardjac=forwardjac, loglevel=loglevel)
5960

6061
if loglevel == debug
6162
println("* solution = $x")
@@ -68,7 +69,7 @@ println("... Test NonlinearEquations.jl")
6869
end
6970

7071

71-
# Determined test from https://www.zib.de/codelib/NewtonLib/nleq_err.tar
72+
# Determined test from https://www.zib.de/codelib/NewtonLib/
7273
@testset "... Test Chebyquad problem of dimensions 2..7" begin
7374

7475
# Test function
@@ -108,11 +109,12 @@ println("... Test NonlinearEquations.jl")
108109
nonlinearity = high
109110
restricted = false
110111
quasi = true
112+
forwardjac = false
111113

112114
for n = 2:nn
113115

114116
if n == 6
115-
continue # NLEQ_ERR does not converge for n = 6, 8, 9
117+
continue
116118
end
117119

118120
if loglevel >= info
@@ -130,7 +132,7 @@ println("... Test NonlinearEquations.jl")
130132
println("* scale vector = $fscale")
131133
end
132134

133-
converged = solveNonlinearEquations!(fx!, n, x, fscale; nonlinearity=high, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, loglevel=loglevel)
135+
converged = solveNonlinearEquations!(fx!, n, x, fscale; nonlinearity=high, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, forwardjac=forwardjac, loglevel=loglevel)
134136

135137
if loglevel == debug
136138
println("* solution = $x")
@@ -179,6 +181,7 @@ println("... Test NonlinearEquations.jl")
179181
nonlinearity = high
180182
restricted = false
181183
quasi = true
184+
forwardjac = false
182185

183186
for ii in 1:2
184187
if loglevel >= info
@@ -196,7 +199,7 @@ println("... Test NonlinearEquations.jl")
196199
println("* scale vector = $fscale")
197200
end
198201

199-
converged = solveNonlinearEquations!(f_2by2!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, loglevel=loglevel)
202+
converged = solveNonlinearEquations!(f_2by2!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, forwardjac=forwardjac, loglevel=loglevel)
200203

201204
if loglevel == debug
202205
println("* solution = $x")
@@ -233,6 +236,7 @@ println("... Test NonlinearEquations.jl")
233236
nonlinearity = high
234237
restricted = false
235238
quasi = true
239+
forwardjac = false
236240

237241
if loglevel >= info
238242
println("MINPACK Powell singular problem:")
@@ -244,7 +248,7 @@ println("... Test NonlinearEquations.jl")
244248
println("* scale vector = $fscale")
245249
end
246250

247-
converged = solveNonlinearEquations!(fx!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, loglevel=loglevel)
251+
converged = solveNonlinearEquations!(fx!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, forwardjac=forwardjac, loglevel=loglevel)
248252

249253
if loglevel == debug
250254
println("* solution = $x")
@@ -280,6 +284,7 @@ println("... Test NonlinearEquations.jl")
280284
nonlinearity = high
281285
restricted = false
282286
quasi = true
287+
forwardjac = false
283288

284289
if loglevel >= info
285290
println("MINPACK Powell badly scaled problem:")
@@ -291,7 +296,7 @@ println("... Test NonlinearEquations.jl")
291296
println("* scale vector = $fscale")
292297
end
293298

294-
converged = solveNonlinearEquations!(fx!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, loglevel=loglevel)
299+
converged = solveNonlinearEquations!(fx!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, forwardjac=forwardjac, loglevel=loglevel)
295300

296301
if loglevel == debug
297302
println("* solution = $x")
@@ -333,6 +338,7 @@ println("... Test NonlinearEquations.jl")
333338
nonlinearity = high
334339
restricted = false
335340
quasi = true
341+
forwardjac = true
336342

337343
if loglevel >= info
338344
println("MINPACK Wood problem:")
@@ -344,7 +350,7 @@ println("... Test NonlinearEquations.jl")
344350
println("* scale vector = $fscale")
345351
end
346352

347-
converged = solveNonlinearEquations!(fx!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, loglevel=loglevel)
353+
converged = solveNonlinearEquations!(fx!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, forwardjac=forwardjac, loglevel=loglevel)
348354

349355
if loglevel == debug
350356
println("* solution = $x")
@@ -390,6 +396,7 @@ println("... Test NonlinearEquations.jl")
390396
nonlinearity = high
391397
restricted = false
392398
quasi = true
399+
forwardjac = false
393400

394401
if loglevel >= info
395402
println("MINPACK helical valey problem:")
@@ -401,7 +408,7 @@ println("... Test NonlinearEquations.jl")
401408
println("* scale vector = $fscale")
402409
end
403410

404-
converged = solveNonlinearEquations!(fx!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, loglevel=loglevel)
411+
converged = solveNonlinearEquations!(fx!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, forwardjac=forwardjac, loglevel=loglevel)
405412

406413
if loglevel == debug
407414
println("* solution = $x")
@@ -423,6 +430,7 @@ println("... Test NonlinearEquations.jl")
423430
nonlinearity = high
424431
restricted = false
425432
quasi = true
433+
forwardjac = true
426434

427435
for n in (6, 9)
428436

@@ -471,7 +479,7 @@ println("... Test NonlinearEquations.jl")
471479
println("* scale vector = $fscale")
472480
end
473481

474-
converged = solveNonlinearEquations!(fx!, n, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, loglevel=loglevel)
482+
converged = solveNonlinearEquations!(fx!, n, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, forwardjac=forwardjac, loglevel=loglevel)
475483

476484
if loglevel == debug
477485
println("* solution = $x")
@@ -517,6 +525,7 @@ println("... Test NonlinearEquations.jl")
517525
nonlinearity = extreme
518526
restricted = true
519527
quasi = false
528+
forwardjac = true
520529
521530
if loglevel >= info
522531
println("MINPACK trigonometric problem:")
@@ -528,7 +537,7 @@ println("... Test NonlinearEquations.jl")
528537
println("* scale vector = $fscale")
529538
end
530539
531-
converged = solveNonlinearEquations!(fx!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, loglevel=loglevel)
540+
converged = solveNonlinearEquations!(fx!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, forwardjac=forwardjac, loglevel=loglevel)
532541
533542
if loglevel == debug
534543
println("* solution = $x")
@@ -561,6 +570,7 @@ println("... Test NonlinearEquations.jl")
561570
nonlinearity = high
562571
restricted = false
563572
quasi = true
573+
forwardjac = false
564574

565575
if loglevel >= info
566576
println("Determined Rosenbrock problem:")
@@ -572,7 +582,7 @@ println("... Test NonlinearEquations.jl")
572582
println("* scale vector = $fscale")
573583
end
574584

575-
converged = solveNonlinearEquations!(fx!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, loglevel=loglevel)
585+
converged = solveNonlinearEquations!(fx!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, forwardjac=forwardjac, loglevel=loglevel)
576586

577587
if loglevel == debug
578588
println("* solution = $x")
@@ -609,6 +619,7 @@ println("... Test NonlinearEquations.jl")
609619
nonlinearity = high
610620
restricted = false
611621
quasi = false
622+
forwardjac = true
612623

613624
if loglevel >= info
614625
println("Underdeterminded Rosenbrock problem:")
@@ -620,7 +631,7 @@ println("... Test NonlinearEquations.jl")
620631
println("* scale vector = $fscale")
621632
end
622633

623-
converged = solveNonlinearEquations!(fx!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, loglevel=loglevel)
634+
converged = solveNonlinearEquations!(fx!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, forwardjac=forwardjac, loglevel=loglevel)
624635

625636
if loglevel == debug
626637
println("* solution = $x")
@@ -663,6 +674,7 @@ println("... Test NonlinearEquations.jl")
663674
nonlinearity = mild
664675
restricted = false
665676
quasi = true
677+
forwardjac = false
666678

667679
if loglevel >= info
668680
println("Slider-crank initial state problem 1:")
@@ -679,7 +691,7 @@ println("... Test NonlinearEquations.jl")
679691
println("* scale vector = $fscale\n")
680692
end
681693

682-
converged = solveNonlinearEquations!(res_z!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, loglevel=loglevel)
694+
converged = solveNonlinearEquations!(res_z!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, forwardjac=forwardjac, loglevel=loglevel)
683695

684696
if loglevel == debug
685697
println("* solution = $x")
@@ -697,7 +709,7 @@ println("... Test NonlinearEquations.jl")
697709

698710

699711
# Underdetermined slider-crank initial state problem 2
700-
@testset "... Test slider-crank initial state problem 1" begin
712+
@testset "... Test slider-crank initial state problem 2" begin
701713

702714
# kinematic parameters
703715
wheelRadius = 0.1
@@ -726,6 +738,7 @@ println("... Test NonlinearEquations.jl")
726738
nonlinearity = high
727739
restricted = false
728740
quasi = true
741+
forwardjac = false
729742

730743
if loglevel >= info
731744
println("Slider-crank initial state problem 2:")
@@ -743,7 +756,7 @@ println("... Test NonlinearEquations.jl")
743756
println("* scale vector = $fscale\n")
744757
end
745758

746-
converged = solveNonlinearEquations!(res!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, loglevel=loglevel)
759+
converged = solveNonlinearEquations!(res!, m, x, fscale; nonlinearity=nonlinearity, restricted=restricted, xtol=tol, maxiter=maxiter, quasi=quasi, forwardjac=forwardjac, loglevel=loglevel)
747760

748761
if loglevel == debug
749762
println("* solution = $x")

0 commit comments

Comments
 (0)