JuliaIntervals / IntervalRootFinding.jl

Library for finding the roots of a function using interval arithmetic
https://juliaintervals.github.io/IntervalRootFinding.jl/
Other
129 stars 25 forks source link

Newton problem with Lagrange points #35

Closed dpsanders closed 6 years ago

dpsanders commented 6 years ago
using StaticArrays

const G = 1.0
const ω = 1 
const μ1 = 0.3 # 1/2
const μ2 = one(μ1) - μ1
const xc1 = -μ2
const xc2 = μ1

fgx(x, y, xc, mu) = G * mu * (x - xc) / ((x - xc)^2+y^2)^(3/2)
fgy(x, y, xc, mu) = G * mu * y / ((x - xc)^2+y^2)^(3/2)

f1(X) = ( (x, y) = X; ω^2*x - fgx(x, y, xc1, μ1) - fgx(x, y, xc2, μ2) )
f2(X) = ( (x, y) = X; ω^2*y - fgy(x, y, xc1, μ1) - fgy(x, y, xc2, μ2) )

f(X) = SVector(f1(X), f2(X))

X = IntervalBox(-10..11, -10..11)

rts = roots(f, X, Bisection);

julia> roots(f, rts, Newton)
7-element Array{IntervalRootFinding.Root{IntervalArithmetic.IntervalBox{2,Float64}},1}:
 Root([1.1232, 1.12321] × [0, 0], :unique)
 Root([-0.200001, -0.199999] × [0.866025, 0.866026], :unique)
 Root([0.299407, 0.300049] × [-0.000518799, 0.000122071], :unknown)
 Root([-0.28613, -0.286129] × [0, 0], :unique)
 Root([-0.700348, -0.699707] × [-0.000518799, 0.000122071], :unknown)
 Root([-0.200001, -0.199999] × [-0.866026, -0.866025], :unique)
 Root([-1.25674, -1.25673] × ∅, :unique)

The last root should have 0 (zero) instead of empty set...

dpsanders commented 6 years ago

Fixed by #55:

julia> roots(f, rts, Newton, 1e-15)
7-element Array{IntervalRootFinding.Root{IntervalArithmetic.IntervalBox{2,Float64}},1}:
 Root([1.1232, 1.12321] × [-1.47863e-28, 1.49292e-28], :unique)
 Root([-0.200001, -0.199999] × [0.866025, 0.866026], :unique)
 Root([0.299999, 0.300001] × [-5.60064e-16, 7.40446e-17], :unknown)
 Root([-0.200001, -0.199999] × [-0.866026, -0.866025], :unique)
 Root([-0.28613, -0.286129] × [-1.54137e-29, 1.54013e-29], :unique)
 Root([-0.700001, -0.699999] × [-5.60064e-16, 7.40446e-17], :unknown)
 Root([-1.25674, -1.25673] × [-4.99325e-29, 4.95997e-29], :unique)