2015-06-29 10 views
5

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

Вот как я пытался это сделать в 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]) 

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

error message

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

+0

Проблема открыта здесь: https://github.com/EconForge/NLsolve.jl/issues/19 –

ответ

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]) не

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