In [43]:

    

function dogleg(f, g, B, x₀, Δ₀, Δₘ, η, ϵ)
    Δ = Δ₀
    x = x₀
    
    while norm(g(x)) > ϵ
        fx = f(x)
        gx = g(x) 
        Bx = B(x) 

        pv = - (gx'gx/(gx'Bx*gx))[1]* gx
        pB = - Bx\gx
        
        if norm(pB) <= Δ
            p = pB
        elseif norm(pv) < Δ
            # Resolver a equação quadrática
            #     ‖pv + (τ-1)(pB-pv)‖² = Δ²
            #     (τ-1)²‖pB-pv‖² + 2pvᵀ(τ-1)(pB-pv) + ‖pv‖² - Δ² = 0
            # ou seja, temos
            #     ax² + 2bx + c = 0
            # onde
            #    x = τ - 1
            #    a = ‖pB-pv‖²
            #    b = pvᵀ(pB-pv)
            #    c = ‖pv‖² - Δ²
            # e
            #    x = (-b + √d)/a
            #    d = b² - ac 
            #      = a(1 - c)
            # com o + porque queremos a solução positiva, já que
            #    x ∈ [0, 1]
            
            a = norm(pB - pv)^2
            b = pv⋅(pB - pv)
            c = norm(pv)^2 - Δ^2
            d = b^2 - a*c
            
            τ = 1 + (-b + √d)/a
            p = pv + (τ - 1)*(pB - pv)
        else
            p = Δ/norm(pv)*pv
        end

        m(p::Vector) = fx + gx'p + 1/2*p'Bx*p |> sum
        
        ρ = (fx - f(x+p))/(m([0;0]) - m(p))
        
        if ρ < 1/4
            Δ = Δ/4
        elseif ρ > 3/4 && norm(p) == Δ
            Δ = min(2Δ, Δₘ)
        end
        
        if ρ > η
            x = x + p
        end
    end
    
    return x
end


    

    
Out[43]:





dogleg (generic function with 1 method)

In [44]:

    

f(x) =       10*(x[2] - x[1]^2)^2 +   (1 - x[1])^2
g(x) = [-40x[1]*(x[2] - x[1]^2)   - 2*(1 - x[1]);
             20*(x[2] - x[1]^2)]

h(x) = [-40*(x[2]-x[1]^2)+80x[1]^2+2    -40x[1];
                    -40x[1]                20  ]
x₀ = [0; -1]


    

    
Out[44]:





2-element Array{Int64,1}:
  0
 -1

In [45]:

    

dogleg(f, g, h, x₀, 1.0, 2.0, 0.5, 1e-5)


    

    
Out[45]:





2-element Array{Float64,1}:
 1.0
 1.0

In [59]:

    

f([1;1])


    

    
Out[59]:





0

In [58]:

    

ϵ = 1e-2
x = linspace(-2, 2, 100)
y = x
contour(x, y, (x, y) -> f([x; y]))


    

    
Out[58]:

In [ ]: