Finding Solution to Inequality in Matlab

I have a system of five equations. As below

syms x y z v w q beta rho theta lambda b tau

eqn1 = x == q*(2*rho/(1-lambda) + beta*y) + (1-q)*((1-theta)*(1-tau)/(1-lambda)+beta*z);
eqn2 = y == q*beta*v + (1-q)*(b*(1-theta)/(1-lambda)+beta*z);
eqn3 = v == q*(2*rho/(1-lambda)+ beta*y) + (1-q)*(b*(1-theta)/(1-lambda)+beta*z);
eqn4 = z == (x + w)/2;
eqn5 = w == q*beta*y + (1-q)*((1-tau)*(1-theta)\(1-lambda)+ beta*z);

I want to find the solution for x, in terms of q, beta, rho, theta, lambda, b and tau.

sol = solve([eqn1, eqn2, eqn3, eqn4, eqn5], [x, y, z, w, v]);
sol.x

After that, I want to find the solution to the following inequality

eqn6 = (1-tau)*(1-theta)/((1-beta)*(1-lambda))> sol.x;

This is where the problem begins. When I run the following command:

solution = solve (eqn6, theta);

I get the following error message:

Warning: Unable to find explicit solution. For options, see help. > In sym/solve (line 317)
In dynamic_model_8 (line 13) 

I get that these equations are complex, but is there a way to get matlaf to find the result?

Thanks in advance.

Answer 1

[Torsten]

I don't understand the result, but there is some output from the symbolic toolbox.

As you can see, the numerator of sol.x is a linear function of theta and the denominator does not contain theta.

Thus the whole thing boils down to determine the zero of a linear function - I'd do it using pencil and paper considering several cases (tau < 1, tau > 1, beta < 1, beta > 1, lambda < 1, lambda > 1) and combinations of these.

I changed

eqn5 = w == q*beta*y + (1-q)*((1-tau)*(1-theta)\(1-lambda)+ beta*z);

to

eqn5 = w == q*beta*y + (1-q)*((1-tau)*(1-theta)/(1-lambda)+ beta*z);

I think it was a typo.

syms x y z v w q beta rho theta lambda b tau
eqn1 = x == q*(2*rho/(1-lambda) + beta*y) + (1-q)*((1-theta)*(1-tau)/(1-lambda)+beta*z);
eqn2 = y == q*beta*v + (1-q)*(b*(1-theta)/(1-lambda)+beta*z);
eqn3 = v == q*(2*rho/(1-lambda)+ beta*y) + (1-q)*(b*(1-theta)/(1-lambda)+beta*z);
eqn4 = z == (x + w)/2;
eqn5 = w == q*beta*y + (1-q)*((1-tau)*(1-theta)/(1-lambda)+ beta*z);
sol = solve([eqn1, eqn2, eqn3, eqn4, eqn5], [x, y, z, w, v]);
sol.x

ans = 

[N,D] = numden(sol.x);
collect(N,theta)

ans = 

eqn6 = (1-tau)*(1-theta)/((1-beta)*(1-lambda)) - sol.x > 0;
solution = solve(eqn6, theta)

solution = 

Answer 2

[John D'Errico]

Simple enough. I'll just work with eqn6 here.

syms tau beta lambda theta Solx
eqn6 = (1-tau)*(1-theta)/((1-beta)*(1-lambda))> Solx;

solve(eqn6,theta,'returnconditions',true)

ans = struct with fields:
         theta: ((x + (tau - 1)/((beta - 1)*(lambda - 1)))*(beta - 1)*(lambda - 1))/(tau - 1)
    parameters: x
    conditions: Solx < x & tau ~= 1

[Dyuman Joshi]

Also, as every other division operator used is a backslash or rdivide, ./, there seems to be a typo here -

%                                             v
qn5 = w == q*beta*y + (1-q)*((1-tau)*(1-theta)\(1-lambda)+ beta*z);

[Torsten]

@John D'Errico

sol.x is a function of theta.

Answer 3

[Walter Roberson]

Turn inequalities into equalities by adding a positive (or non-negative) variable representing how much more one side is compared to the other.

syms x y z v w q beta rho theta lambda b tau

eqn1 = x == q*(2*rho/(1-lambda) + beta*y) + (1-q)*((1-theta)*(1-tau)/(1-lambda)+beta*z);
eqn2 = y == q*beta*v + (1-q)*(b*(1-theta)/(1-lambda)+beta*z);
eqn3 = v == q*(2*rho/(1-lambda)+ beta*y) + (1-q)*(b*(1-theta)/(1-lambda)+beta*z);
eqn4 = z == (x + w)/2;
eqn5 = w == q*beta*y + (1-q)*((1-tau)*(1-theta)\(1-lambda)+ beta*z);
sol = solve([eqn1, eqn2, eqn3, eqn4, eqn5], [x, y, z, w, v]);
sol.x

ans = 

syms Excess positive
eqn6 = (1-tau)*(1-theta)/((1-beta)*(1-lambda)) - Excess == sol.x

eqn6 = 

solution = solve (eqn6, theta)

solution = 

sympref('AbbreviateOutput', 0)

ans = logical
   1

collect(simplify(solution, 'steps', 50), Excess)

ans = 

[Sam Chak]

Wow! I wonder how to write such super-long equations in the journal paper or thesis.

[Walter Roberson]

I notice that the collect() is not working inside the sqrt()

syms x Excess
collect((Excess^2 + 5*Excess + sin(x)*Excess), Excess)

ans = 

collect(sqrt(Excess^2 + 5*Excess + sin(x)*Excess), Excess)

ans = 

[Sam Chak]

Not only sqrt(). It also doesn't work inside the sin() function.

Perhaps the collect() only works on the symbolic polynomial expression, and the user-defined symbolic polynomial function, but not on transcendental expressions or functions.

syms a b x y son f(x, y)

theChildren   = a^2*x*y + a*b*x^2 + a*x*y + x^2;
coeffs_xy     = collect(a^2*x*y + a*b*x^2 + a*x*y + x^2, [x y])

coeffs_xy = 

childrenOfSin = collect(sin(theChildren), [x y])    % doesn't work inside a function

childrenOfSin = 

%% Further test using tanh(x) = exp(2*x) - 1)/(exp(2*x) + 1)
f(x, y) = a^2*x*y + a*b*x^2 + a*x*y + x^2;
collect(tanh(f(x, y)), [x y])

ans = 

expr    = (exp(2*f(x, y)) - 1)/(exp(2*f(x, y)) + 1)

expr = 

out     = collect(expr, [x y])

out = 

simplify(out, 'Steps', 50)

ans = 

%% Other tests
whosDaughters = collect(theChildren, x^2)

whosDaughters = 

whosSons      = collect(theChildren, x*y)   % doesn't work as intended

whosSons = 

whosSons      = subs(collect(subs(whosDaughters, x*y, son), son), son, x*y)

whosSons = 

groupChildren = children(whosSons)

groupChildren = 1×2 cell array
    {[x^2*(a*b + 1)]}    {[x*y*(a^2 + a)]}

[Walter Roberson]

Looks like it does not descend into function calls except possibly _power calls with integer coefficients.

But notice also the coeffs3 result of collect on ^2 factors down to (a*b+1)^2 but that does not happen when you collect on ^-2 or 1/^2

syms a b x y son f(x)
theChildren   = a^2*x*y + a*b*x^2 + a*x*y + x^2;
coeffs1     = collect(theChildren, [x y])

coeffs1 = 

coeffs2 = collect(f(theChildren), [x y])

coeffs2 = 

coeffs3 = collect(theChildren^2, [x y])

coeffs3 = 

coeffs4 = collect(1/(theChildren^2), [x y])

coeffs4 = 

coeffs5 = collect(theChildren^(-2), [x y])

coeffs5 = 

1/coeffs3

ans = 

ans - coeffs4

ans = 

simplify(ans)

ans = 

0

[Sam Chak]

Thanks @Walter Roberson. I guess that your tests further reinforce that the 'collect()' function only works on strictly polynomials, , but not on rational functions, even though they are quotients of two polynomials, .

[Walter Roberson]

If you look at coeffs4 the denominator is collected with respect to x: it just isn't as factored as for the coeffs3 case. It did not fail to collect: it just didn't fully factor after collecting, even though fully factoring is not a documented outcome in the first place.