2017-02-15 63 views
1

У меня есть три уравнения в трех неизвестных, которые я хотел бы решить. Я задаю уравнения с symbolic toolbox. Я знаю, что могу использовать функцию solve, чтобы спросить Matlab, чтобы найти мне числовое решение. Однако с 3 уравнениями в 3 неизвестных, Matlab должен иметь возможность найти аналитическое решение (fsolve). Я просто не уверен, как изменить код, чтобы я мог использовать fsolve вместо solve.От Solve to Fsolve

Ниже мой код:

очистить все

syms Kl Kh alpha nu w phi delta P beta zh zl Ezh Ezl 

nu1 = (1/(1-nu)); 

f1 = ((zl * (Kl^alpha))^nu1 + (zh * (Kh^alpha))^nu1) * nu^(nu*nu1) * (w^(-nu*nu1)) - w/phi + delta*(Kl + Kh)*P 
f2 = Kh - (((1-beta*(1-delta))*P * (w^(nu1*nu))*(nu^(nu*nu1)))/(beta*alpha* (Ezh)^nu1))^((1-nu)/(alpha+nu-1)) 
f3 = Kl - (((1-beta*(1-delta))*P * (w^(nu1*nu))*(nu^(nu*nu1)))/(beta*alpha* (Ezl)^nu1))^((1-nu)/(alpha+nu-1)) 

f1 = subs(f1, {alpha, beta, nu, phi,delta, zh, zl, Ezh, Ezl, P}, {0.27, 0.96, 0.60, 2.15,0.065,1.11687642219068,0.895354204038589,1.07811003137331,0.934120594855956, 0.95}) 
f2 = subs(f2, {alpha, beta, nu, phi,delta, zh, zl, Ezh, Ezl, P}, {0.27, 0.96, 0.60, 2.15,0.065,1.11687642219068,0.895354204038589,1.07811003137331,0.934120594855956, 0.95}) 
f3 = subs(f3, {alpha, beta, nu, phi,delta, zh, zl, Ezh, Ezl, P}, {0.27, 0.96, 0.60, 2.15,0.065,1.11687642219068,0.895354204038589,1.07811003137331,0.934120594855956, 0.95}) 



S = solve([f1 == 0, f2 == 0, f3 == 0],... 
    [w, Kh, Kl], 'ReturnConditions', true); 
+2

'fsolve' используются, чтобы найти * * численное решение, а не аналитические. 'solve', в символической математической панели инструментов, может найти аналитические решения (если они существуют) в некоторых случаях, но другие мудрые будут использовать математику с переменной точностью, чтобы найти числовую. Вы ищете числовое решение? – horchler

+0

Да, мои извинения. Я имел в виду численное решение. – phdstudent

ответ

0

В то же время я нашел решение.

Вот код завершения:

function SSfunction = SSfunction(x) 


syms alphaa nu phi delta p betaa zh zl ezh ezl 

nu1 = (1/(1-nu)); 

f1 = ((zl * (x(3)^alphaa))^nu1 + (zh * (x(2)^alphaa))^nu1) * nu^(nu*nu1) * (x(1)^(-nu*nu1)) - x(1)/phi - delta*(x(3) + x(2))*p; 
f2 = x(2) - ((betaa*alphaa*(ezh^(nu1)) * (nu^(nu*nu1)))/((1-betaa*(1-delta))*p* (x(1)^(nu*nu1))))^((1-nu)/(1-alphaa-nu)); 
f3 = x(3) - ((betaa*alphaa*(ezl^(nu1)) * (nu^(nu*nu1)))/((1-betaa*(1-delta))*p* (x(1)^(nu*nu1))))^((1-nu)/(1-alphaa-nu)); 

f1 = subs(f1, {alphaa, betaa, nu, phi,delta, zh, zl, ezh, ezl, p}, {0.27, 0.96, 0.60, 2.15,0.065,1.11687642219068,0.895354204038589,1.07811003137331,0.934120594855956, 0.950}); 
f3 = subs(f1, {alphaa, betaa, nu, phi,delta, zh, zl, ezh, ezl, p}, {0.27, 0.96, 0.60, 2.15,0.065,1.11687642219068,0.895354204038589,1.07811003137331,0.934120594855956, 0.950}); 
f2 = subs(f1, {alphaa, betaa, nu, phi,delta, zh, zl, ezh, ezl, p}, {0.27, 0.96, 0.60, 2.15,0.065,1.11687642219068,0.895354204038589,1.07811003137331,0.934120594855956, 0.950}); 



SSfunction(1) = eval(f1) 
SSfunction(2) = eval(f2) 
SSfunction(3) = eval(f3) 


end 

x0 = [1,2,0.7]; 
fun = @SSfunction; 
x = fsolve(fun,x0) 
+0

Два комментария. Во-первых, как правило, плохой практикой называют вывод функции такой же, как и сама функция ('SSfunction'). Кроме того, это решение неэффективно, потому что символические выражения преобразуются в плавающую точку каждый раз, когда оптимизация 'fsolve' называется целевой функцией. В очень особых случаях этот метод может быть полезен, если вы потеряете значимость или другую численную проблему в уравнениях. Это, похоже, не так. – horchler

1

Использование matlabFunction для преобразования символьных выражений в векторизованного числовой функции, которые могут быть непосредственно используемой fsolve:

... 
f = matlabFunction([f1;f2;f3],'Vars',{[w;Kh;Kl]}); 

w0 = 1; 
Kh0 = 1; 
Kl0 = 1; 
x0 = [w0;Kh0;Kl0]; 
x = fsolve(f,x0) 

Это будет порядок больше, чем использование самих символических выражений в fsolve. Для более высокой скорости, вы также можете просто избавиться от символической математики полностью векторизаций вашей функции вручную:

alpha = 0.27; 
beta = 0.96; 
nu = 0.6; 
phi = 2.15; 
delta = 0.065; 
zh = 1.11687642219068; 
zl = 0.895354204038589; 
Ezh = 1.07811003137331; 
Ezl = 0.934120594855956; 
P = 0.95; 
nu1 = (1/(1-nu)); 

f = @(w,Kh,Kl)[((zl * (Kl.^alpha))^nu1 + (zh * (Kh.^alpha))^nu1) * nu^(nu*nu1) .* (w.^(-nu*nu1)) - w/phi + delta*(Kl + Kh)*P; 
       Kh - (((1-beta*(1-delta))*P * (w.^(nu1*nu))*(nu^(nu*nu1)))/(beta*alpha* (Ezh)^nu1))^((1-nu)/(alpha+nu-1)); 
       Kl - (((1-beta*(1-delta))*P * (w.^(nu1*nu))*(nu^(nu*nu1)))/(beta*alpha* (Ezl)^nu1))^((1-nu)/(alpha+nu-1))]; 

w0 = 1; 
Kh0 = 1; 
Kl0 = 1; 
x0 = [w0;Kh0;Kl0]; 
x = fsolve(@(x)f(x(1,:),x(2,:),x(3,:)),x0)