Search code examples
juliaiterationlinear-algebraequation-solvinglinear-equation

Solve a system of N equations with N unknowns using Julia


I have :

  • a set of N locations which can be workplace or residence
  • a vector of observed workers L_i, with i in N
  • a vector of observed residents R_n, with n in N
  • a matrix of distance observed between all pair residence n and workplace i
  • a shape parameter epsilon

Setting N=3, epsilon=5, and

d = [1 1.5 3 ; 1.5 1 1.5 ; 3 1.5 1] #distance matrix
L_i = [13  69  18] #vector of workers in each workplace
R_n = [27; 63; 10]

I want to find the vector of wages (size N) that solve this system of N equations, enter image description here

with l all the workplaces.

Do I need to implement an iterative algorithm on the vectors of workers and wages? Or is it possible to directly solve this system ?

I tried this,

w_i = [1 ; 1 ; 1]
er=1
n =1
while  er>1e-3
    L_i =   ( (w_i ./ d).^ϵ ) ./ sum( ( (w_i ./ d).^ϵ), dims=1) * R 
    er = maximum(abs.(L .- L_i))
    w_i = 0.7.*w_i + 0.3.*w_i.*((L .- L_i) ./ L_i)
    n = n+1
end

Solution

  • If L and R are given (i.e., do not depend on w_i), you should set up a non-linear search to get (a vector of) wages from that gravity equation (subject to normalising one w_i, of course).

    Here's a minimal example. I hope it helps.

    # Call Packages
    using Random, NLsolve, LinearAlgebra
    
    # Set seeds
    Random.seed!(1704)
    
    # Variables and parameters
    N    = 10
    R    = rand(N)
    L    = rand(N) * 0.5
    d    = ones(N, N) .+ Symmetric(rand(N, N)) / 10.0 
    d[diagind(d)] .= 1.0
    ε    = -3.0
    
    # Define objective function
    function obj_fun(x, d, R, L, ε) 
        
        # Find shares
        S_mat  = (x ./ d).^ε
        den    = sum(S_mat, dims = 1)
        s      = S_mat ./ den
    
        # Normalize last wage
        x[end] = 1.0
    
        # Define loss function
        loss   = L .- s * R
    
        # Return 
        return loss
        
    end
    
    # Run optimization
    x₀   = ones(N) 
    res  = nlsolve(x -> obj_fun(x, d, R, L, ε), x₀, show_trace = true)
    
    # Equilibrium vector of wages
    w    = res.zero