Incoherence between exact and finite differences Jacobian

New to Julia
2

Hey @GuyFromBosnia, welcome!

To make it easier for others to help you, it is a good idea to post a reproducible code example that is cleanly formatted (see here for more advice).

I took the liberty of adjusting formatting and making some simplifications, here’s what I arrive at:

using LinearAlgebra

function f(v)
    n = length(v)
    h = 1 / (n + 1)
    y = zeros(n)
    y[1] = -2v[1] + v[2] + (1 / 4) * exp(v[1])
    y[n] = -2v[n] + v[n-1] + (h^2) * exp(v[n])
    for i in 2:n-1
        h = 1 / (i + 1)
        y[i] = ((v[i+1] - 2v[i] + v[i-1]) / h^2) + exp(v[i])
        y[i] *= h^2
    end
    return y
end

function jacobian_exact(f, x)
    n = length(x)
    J = zeros(n, n)
    for i = 1:n
        for j = 1:n
            h = 1 / (i + 1)
            if i == j
                J[i, j] = 2 / h^2 + exp(x[i])
            elseif i == j + 1 || i == j - 1
                J[i, j] = -1 / h^2
            end
        end
    end
    return J
end

function jacobian_approx(f, x, epsilon)
    n = length(x)
    J = zeros(n, n)
    for j in 1:n
        ej = zeros(n)
        ej[j] = 1
        J[:, j] = (f(x + epsilon * ej) - f(x)) / epsilon
    end
    return J
end

function newton_method(f, x0, epsilon, max_iter, method::Symbol=:exact)
    x = copy(x0)
    iter = 0
    while norm(f(x)) > epsilon && iter < max_iter
        if method == :exact
            J = jacobian_exact(f, x)
        elseif method == :approximate
            J = jacobian_approx(f, x, epsilon)
        else
            error("Invalid method specified. Must be :exact or :approximate.")
        end
        x -= J \ f(x)
        iter += 1
    end
    return x
end

epsilon = 1e-6
xr = collect(1:4)
jacobian_exact(f, xr)
jacobian_approx(f, xr, epsilon)

x0 = zeros(4)
exact_root_newton = newton_method(f, x0, epsilon, 100, :exact)
approx_root_newton = newton_method(f, x0, epsilon, 100, :approximate)

From what I can tell, the jacobian_approx function seems to be correct, so I would say there is either an error in jacobian_exact or in the function f itself. What is the function you’re trying to implement?

If you’re new to Julia, you may not be aware of it, but there are plenty of packages that will compute exact derivatives for you much more reliably: see the JuliaDiff webpage for a list

5 Likes