In [1]:

using Plots
gr(size=(400, 300), fmt=:png);

In [2]:

using NLsolve

In [3]:

using ForwardDiff

In [16]:

using LinearAlgebra

In [ ]:

In [5]:

a = Ref{Float64}(2.)

function f!(du, u)
    x, y = u
    du[1] = a[] / (1 + y^2) - x
    du[2] = a[] / (1 + x^2) - y
    nothing
end

Out[5]:

f! (generic function with 1 method)

In [6]:

a[] = 3.
r = nlsolve(f!, [1., 4.], autodiff=:forward)

Out[6]:

Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [1.0, 4.0]
 * Zero: [0.381966, 2.61803]
 * Inf-norm of residuals: 0.000000
 * Iterations: 5
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 6
 * Jacobian Calls (df/dx): 6

In [7]:

r.zero

Out[7]:

2-element Array{Float64,1}:
 0.38196601125010465
 2.6180339887498962

In [ ]:

In [8]:

n = 10000

as = rand(n) * 6 # list of a
result = rand(n, 2) * 6 # use as initial u

@time for i in 1:n
    
    a[] = as[i]
    r = nlsolve(f!, result[i, :], autodiff=:forward)
    
    if r.f_converged
        result[i, :] = r.zero
    end
end

scatter(as, result[:, 2],
    markeralpha=0.2, markersize=2, markerstrokewidth=0,
    legend=:none, xlims=(0, 6), ylims=(0, 6))

  0.338655 seconds (2.22 M allocations: 88.334 MiB, 8.83% gc time)

Out[8]:

In [9]:

scatter(as, result[:, 1], result[:, 2],
    xlims=(0, 6), ylims=(0, 6), zlims=(0, 6),
    markeralpha=0.2, markersize=2, markerstrokewidth=0,
    legend=:none)

Out[9]:

In [ ]:

In [17]:

negative(x) = x < 0
stable(J) = reduce(&, negative.(eigvals(J)))

Out[17]:

stable (generic function with 1 method)

In [18]:

a[] = 4.
@show r = nlsolve(f!, [3., 3.], autodiff=:forward)

@show J = ForwardDiff.jacobian(f!, similar(r.zero), r.zero)
@show eigvals(J)
@show stable(J)

r = nlsolve(f!, [3.0, 3.0], autodiff=:forward) = Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [3.0, 3.0]
 * Zero: [1.3788, 1.3788]
 * Inf-norm of residuals: 0.000000
 * Iterations: 5
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 6
 * Jacobian Calls (df/dx): 6
J = ForwardDiff.jacobian(f!, similar(r.zero), r.zero) = [-1.0 -1.3106; -1.3106 -1.0]
eigvals(J) = [0.310602, -2.3106]
stable(J) = false

Out[18]:

false

In [ ]:

In [21]:

stability = zeros(Int, n)
dummy = zeros(2)

@time for i in 1:n
    
    a[] = as[i]
    J = ForwardDiff.jacobian(f!, dummy, result[i, :])
    
    if stable(J)
        stability[i] = 1
    end
end

scatter(as, result[:, 2], color=stability,
    markeralpha=0.2, markersize=2, markerstrokewidth=0,
    legend=:none, xlims=(0, 6), ylims=(0, 6))

  0.143750 seconds (407.72 k allocations: 67.396 MiB, 24.38% gc time)

Out[21]:

In [ ]: