Нелинейная система уравнений Юлия

Я пытаюсь решить большое количество (50) нелинейных одновременных уравнений в Julia. На данный момент я просто пытаюсь сделать эту работу с двумя уравнениями, чтобы получить правильный синтаксис и т. Д. Однако я пробовал множество пакетов / инструментов - NLsolve, nsolve в SymPy и NLOpt в JuMP (где я игнорирую цель функция и просто введите ограничения равенства) - без особого успеха. Думаю, мне стоит сосредоточиться на том, чтобы все работало в одном. Буду признателен за любые советы по выбору пакетов и, если возможно, кода.

Вот как я пытался сделать это в NLsolve (используя его в режиме mcpsolve, чтобы я мог наложить ограничения на переменные, которые я вычисляю - x [1] и x [2] - которые являются уровнями безработицы и поэтому ограничены между нулем и 1) :

using Distributions
using Devectorize
using Distances
using StatsBase
using NumericExtensions
using NLsolve

beta = 0.95                                                                
xmin= 0.73                                                                 
xmax = xmin+1                                                             
sigma = 0.023                                                            
eta = 0.3                                        
delta = 0.01                                                                                                
gamma=0.5                                                                   
kappa = 1                                                                  
psi=0.5
ns=50
prod=linspace(xmin,xmax,ns)
l1=0.7
l2=0.3                                            
wbar=1
r=((1/beta)-1-1e-6 +delta)


## Test code

function f!(x, fvec)

    ps1= wbar + (kappa*(1-beta*(1-sigma*((1-x[1])/x[1]))))
    ps2= wbar + (kappa*(1-beta*(1-sigma*((1-x[2])/x[2]))))

    prod1=prod[1]
    prod2=prod[50]
    y1=(1-x[1])*l1
    y2=(1-x[2])*l2
    M=(((prod1*y1)^((psi-1)/psi))+((prod2*y2)^((psi-1)/psi)))
    K=((r/eta)^(1/(eta-1)))*M

    pd1=(1-eta)*(K^eta)*(M^(-eta))*prod1
    pd2=(1-eta)*(K^eta)*(M^(-eta))*prod2

    fvec[1]=pd1-ps1
    fvec[2]=pd2-ps2
end

mcpsolve(f!,[0.0,0.0],[1.0,1.0], [ 0.3, 0.3])

Я получаю это сообщение об ошибке:

сообщение об ошибке

Любые предложения приветствуются! Я ценю, что формулы довольно уродливые, поэтому дайте мне знать, будут ли полезны дальнейшие упрощения (я пробовал!).


person David Zentler-Munro    schedule 29.06.2015    source источник
comment
Проблема открыта здесь: github.com/EconForge/NLsolve.jl/issues/19   -  person Simon Byrne    schedule 06.09.2016


Ответы (1)


Я думал, что вы задаете начальные условия за пределами диапазона, потому что я попробовал mcpsolve(f!,[0.0,0.0],[0.0,0.0],[0.3, 0.3]), и это сработало.

Однако я пробовал и другие комбинации:

mcpsolve(f!,[0.4,0.4], [0.0,0.0], [0.3, 0.3]) действительно сработало

mcpsolve(f!,[0.4,0.4], [0.3,0.3], [1.0,1.0]) не

mcpsolve(f!,[0.6,0.6], [1.0,1.0], [0.3,0.3]) не

Вы проверили эти значения в своем тесте?

person jmeloc    schedule 31.07.2015