From Solve to Fsolve

Question

I have three equations in three unknowns that I would like to solve. I am specifying the equations with symbolic toolbox. I know I can use solve function to ask matlab to find me a numeric solution. However, with 3 equations in 3 unknowns, matlab should be able to find an analytical solution (fsolve). I am just not sure how to change the code so that I can use fsolve instead of solve.

Below my code:

clear all

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);

`fsolve` is used to find a *numerical* solution, not analytical ones. `solve`, in the symbolic math toolbox, can find analytical solutions (if they exist) in some cases, but other wise will use variable precision math to find a numeric one. Are you looking for a numeric solution? — horchler, Feb 15 '17 at 21:08

horchler · Answer 1 · 2017-02-16T15:01:06.013

Use matlabFunction to convert your symbolic expressions to a vectorized numeric function that can be directly used by fsolve:

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

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

This will be an order of magnitude faster than using the the symbolic expressions themselves in fsolve. For more speed, you could also just get rid of symbolic math entirely by vectorizing your function manually:

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)

score 0 · Accepted Answer · answered Feb 16 '17 at 13:27

In the meantime I found the solution.

Here's the code for completion:

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)

Two comments. First, it generally a bad practice to name the output of a function the same as the function itself (`SSfunction`). Also, this solution is inefficient because the symbolic expressions are converted to floating point each time the `fsolve` optimization called the objective function. In very special cases, this method might be useful if you have a loss of significance or other numerical issue in the equations. That doesn't appear to be the case here. — horchler, Feb 16 '17 at 15:07

From Solve to Fsolve

2 Answers2