No result from nlsolve

New to Julia
,
1

Below is a MWE of a problem I am working with. First, a function is defined with a system of 4 equations. Second, matrices are initialized. Third, I try to solve the system using nlsolve. After running this MWE, the code stops without error, but nothing is saved in results? Any clue where it might be going wrong? Thanks in advance.

function sys_tHE_AB!(F, Z_AB, k_e, fe, Fe) # k_e: (1000, 10, 11)  fe: (11, 11)  Fe:  (11, 11)  
        ρ = 0.05 # rate of time preference
        δ = 0.05 # exogenous death rate (here also endogenous exit)
        β = sqrt(2*(ρ + δ))+5
        α = -sqrt(2*(ρ + δ))
        disc = ρ + δ
        F[1] =   Z_AB[3]*max(Z_AB[1],0)^α +        k_e*Z_AB[1]/disc - fe/disc - Fe - Z_AB[4]*max(Z_AB[1],0)^β
        F[2] =   Z_AB[3]*max(Z_AB[2],0)^α +        k_e*Z_AB[2]/disc - fe/disc      - Z_AB[4]*max(Z_AB[2],0)^β
        F[3] = α*Z_AB[3]*max(Z_AB[1],0)^(α-1) +    k_e/disc    - β*Z_AB[4]*max(Z_AB[1],0)^(β-1)
        F[4] = α*Z_AB[3]*max(Z_AB[2],0)^(α-1) +    k_e/disc    - β*Z_AB[4]*max(Z_AB[2],0)^(β-1)
        N = [F[1] F[2] F[3] F[4]]     
end

n=20
N=10
k_e=zeros(n,N,N+1) # 20 x 10 x 11
fe=rand(11,11) # 11 x 11
Fe=rand(11,11) # 11 x 11
Z_HE_0=1*rand(20,10,11) # 20 x 10 x 11
Z_EH_0=2*rand(20,10,11) # 20 x 10 x 11
Bo_0=3*rand(20,10,11) # 20 x 10 x 11
init_HE=rand(20,10,11,4) # 20 x 10 x 11 x4

for i in 1:20
        for j in 2:10
                for k in 3:11
                        init_HE[i,j,k,:] = [Z_HE_0[i,j,k] Z_EH_0[i,j,k] 0 Bo_0[i,j,k]] # 1 x 4 vector
                        result = nlsolve((F, Z_AB) -> sys_tHE_AB!, init_HE[i,j,k,:], method = :trust_region, xtol=10e-5, ftol=10e-8, iterations=10, store_trace=true) 
                end
        end
end
2

There is no variable called results in your code. You have a variable called result in your loop, but that is local to the loop. Make sure you define this before you enter the loop, and then either push! into it or preallocate the right length (and type) structure in the first place.

The simplest change would be:

result = []
for i in 1:2
      for j in 2:3
              for k in 3:4
                      init_HE[i,j,k,:] = [Z_HE_0[i,j,k] Z_EH_0[i,j,k] 0 Bo_0[i,j,k]] # 1 x 4 vector
                      push!(result, nlsolve((F, Z_AB) -> sys_tHE_AB!, init_HE[i,j,k,:], method = :trust_region, xtol=10e-5, ftol=10e-8, iterations=10, store_trace=true))
              end
      end
end

which gives

julia> result
8-element Vector{Any}:
 Results of Nonlinear Solver Algorithm
...

This isn’t ideal for performance as result is untyped and push!ing into it can be slow, although it might be a negligible loss in performance relative to the cost of the solve.

1 Like
3

Thank you so much @nilshg!! That is extremely helpful!!

1 Like