JuliaNLSolvers / NLsolve.jl

Julia solvers for systems of nonlinear equations and mixed complementarity problems
Other
324 stars 66 forks source link

Bizarre Bug #241

Closed janrpeters closed 4 years ago

janrpeters commented 4 years ago

I have a piece of illustrative code that works flawlessly:

using NLsolve;

f(x) = [x[2]; -9.81*sin(x[1])/16.0];

for i=1:1000
    x = [-0.5243754571784831; 0.5382977006541599]

    function f_root!(F_root, y)
        h = 0.01;
        w = f(y);
        F_root[1] = (y[1] - x[1])/h - w[1];
        F_root[2] = (y[2] - x[2])/h - w[2];
    end
    z = nlsolve(f_root!, x).zero 
    println("z[",i,"]=",z)
end

Surprisingly, vectorizing it yields a nondeterministic bug:

using NLsolve;

f(x) = [x[2]; -9.81*sin(x[1])/16.0];

for i=1:1000
    x = [-0.5243754571784831; 0.5382977006541599]

    function f_root!(F_root, y)
        h = 0.01;
        w = f(y);
        F_root = (y - x)./h - w;
    end
    z = nlsolve(f_root!, x).zero 
    println("z[",i,"]=",z)
end
antoine-levitt commented 4 years ago

F_root = ... only changes the binding of F_root, not its value. Use F_root .=.

janrpeters commented 4 years ago

Thanks!