Search code examples
matlabsymbolic-mathsolverequation-solving

From Solve to Fsolve


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

Solution

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