Incoherence between exact and finite differences Jacobian

[GuyFromBosnia]

I have this code here where I need to compute the exact jacobian matrix and the approximate jacobian matrix using the newton’s method I have the exact jacobian matrix giving me the correct results but the approximate is giving me the wrong results.I hope that you can help me
function f(v)
n = length(v) # stores the length of v
h = 1/(n+1)
f = zeros(n)
f[1] = -2v[1] + v[2] + (1/4)*exp(v[1])
f[n] = -2v[n] + v[n-1] + (h^2)*exp(v[n])
for i in 2:n-1
local h=1/(i+1)
f[i] = ((v[i+1] - 2v[i] + v[i-1])/h^2) + exp(v[i])
f[i] *= h^2
end
return f
end

epsilon=0.000001

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

function find_jacobian_matrix(f, x, epsilon)
n = length(x)
J = zeros(n, n)
for j in 1:n
ei = zeros(n)
ei[j] = 1
J[:, j] = (f(x + epsilon*ei) - 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 = find_jacobian_matrix(f, x, epsilon)
else
error(“Invalid method specified. Must be :exact or :approximate.”)
end
y = -J \ f(x)
x += y
iter += 1
end
return x
end

initial guess

x0 = [0, 0, 0, 0]

call the function

exact_root = jacobian_exact(f, x0) \ f(x0)
approx_root = find_root(f, x0, epsilon, 100, :approximate)

display the results

println("Exact Root: ", exact_root)
println("Approximate Root: ", approx_root)

[gdalle]

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

[gdalle]

Also, could you maybe modify the title of your question to something more specific, like “Incoherence between exact and finite differences Jacobian”? Basically everyone posting here needs help with a task

[nilshg]

The OPs problem is a homework problem, so he has to solve it without library support. That said a library can of course be helpful to check results.

OP make sure to let people know that you are looking for help with academic coursework so that people responding can tailor their replies appropriately and you’re not violating any academic integrity rules.

[GuyFromBosnia]

The second function is the one I am trying to implement i have fixed the exact jacobian function so I am getting almoast the same ones as approximated but I am not sure if the results are correct