How to define this variable on MATLAB

I have this positive equilibrium (x,l,y,v,M) where

x = (q/b)*(eta+(eps*M)/(1+M))*(1+mu+((rho*M)/(1+M)))

l = (q/(alp*sig))*(zeta*M-omega0)

y = (1/sig)*(zeta*M-omega0)

v = (b/sig)*(1/(eta+(eps*M)/(1+M)))*(zeta*M-omega0)

M is a positive root of the following equation

r0*(1+k1*(mu+(rho*M)/(1+M))*(1-k2*(mu+(rho*M)/(1+M))))*(eta+(eps*M)/(1+M))*(K-alp*sig*q*(eta+(eps*M)/(1+M))*(1+mu+(rho*M)/(1+M))-b*(q+alp)*(zeta*M-omega0))-b^2*alp*K*(zeta*M-omega0)=0

Does anyone have an idea of how to enter these informaion on MATLAB, espically the last three lines?

Answer 1

[Suvansh Arora]

In order to solve this equation with 5 unknows, follow the MATLAB answers article mentioned below:

How to solve 5 equations with 5 unknowns - MATLAB Answers - MATLAB Central

[Yakoob]

Thanks Suvansh for your reply.

I am not trying to solve these equations. I am looking for the right command that allows MATLAB to recognize that M* is a positive root of the equation above. Please note that

r0=0.05;

k1=0.5;

mu=0.5;

rho=0.5;

k2=0.5;

eta=0.05;

eps=0.25;

K=1;

alp=0.1;

sig=0.05

q=0.5;

b=0.1;

zeta=4;

omega0=1;

[Steven Lord]

If you are able to specify that, what are you going to do with that variable and your other equations? Knowing how you're planning to use that information can help us offer suggestions for the best way to achieve your ultimate goal.

[Steven Lord]

If you want to assign a sym object into an array and operate symbolically, preallocate the array as a sym array rather than as a double. If you then wanted to pass that into a function that only accepts double arrays you may need to convert at that point (which may involve using subs to substitute values for the symbolic variables.)

A = zeros(3, 4, 'sym')

A = 

syms x
A(2, 2) = x^2 % this works

A = 

Compare with:

B = zeros(3, 4)

B = 3×4

0 0 0 0 0 0 0 0 0 0 0 0

B(2, 2) = x^2 % this throws an error

Unable to perform assignment because value of type 'sym' is not convertible to 'double'.

Caused by:
Error using symengine
Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for variables.

Answer 2

[Paul]

Hi @Yakoob

We can use the Symbolic Math Toolbox as one option

Set up the equations

r0=0.05;
k1=0.5;
mu=0.5;
rho=0.5;
k2=0.5;
eta=0.05;
eps=0.25;
K=1;
alp=0.1;
sig=0.05;
q=0.5;
b=0.1;
zeta=4;
omega0=1;

syms M positive

x = (q/b)*(eta+(eps*M)/(1+M))*(1+mu+((rho*M)/(1+M)));
l = (q/(alp*sig))*(zeta*M-omega0);
y = (1/sig)*(zeta*M-omega0);
v = (b/sig)*(1/(eta+(eps*M)/(1+M)))*(zeta*M-omega0);

eqn = r0*(1+k1*(mu+(rho*M)/(1+M))*(1-k2*(mu+(rho*M)/(1+M))))*(eta+(eps*M)/(1+M))*(K-alp*sig*q*(eta+(eps*M)/(1+M))*(1+mu+(rho*M)/(1+M))-b*(q+alp)*(zeta*M-omega0))-b^2*alp*K*(zeta*M-omega0)==0

eqn = 

eqn = simplify(lhs(eqn)) == 0

eqn = 

Try solve

Msol1 = solve(eqn,M)

Msol1 = 

Use vpa to get the root of the polynomial

vpa(Msol1)

ans = 

2.1265555666349343861231037078585

Verify the solution

subs(lhs(eqn),M,vpa(Msol1))

ans = 

Alternatively, we can extract the numerator of the LHS and use roots, there is only ony positive solution:

Msol2 = roots(sym2poly(numden(lhs(eqn))))

Msol2 = 6×1

2.1266 -1.2184 -1.0161 -0.9851 -0.7685 -0.2448

If the equation to solve was more complicated then just a sixth order polynmial, we could use vpasolve

Msol3 = vpasolve(eqn,M,[0 inf])

Msol3 = 

2.1265555666349343861231037078585

Or if desired to just get a numerical solutiion, use fsolve. The function to solve could have been defined directly without going through the symbolic stuff, but since we already have the symbolic stuff, just use it

fun = matlabFunction(lhs(eqn));
Msol4 = fsolve(fun,1)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

Msol4 = 2.1285

fun(Msol4)

ans = -1.2185e-05

[Yakoob]

Thanks Paul for your response.

What if one of the paramters, let's say sig, takes this form: sig=linspace(0.01,5,20)?

[Paul]

r0=0.05;
k1=0.5;
mu=0.5;
rho=0.5;
k2=0.5;
eta=0.05;
eps=0.25;
K=1;
alp=0.1;
%sig=0.05;
q=0.5;
b=0.1;
zeta=4;
omega0=1;
syms M sigma positive
x = (q/b)*(eta+(eps*M)/(1+M))*(1+mu+((rho*M)/(1+M)));
l = (q/(alp*sigma))*(zeta*M-omega0);
y = (1/sigma)*(zeta*M-omega0);
v = (b/sigma)*(1/(eta+(eps*M)/(1+M)))*(zeta*M-omega0);
eqn = r0*(1+k1*(mu+(rho*M)/(1+M))*(1-k2*(mu+(rho*M)/(1+M))))*(eta+(eps*M)/(1+M))*(K-alp*sigma*q*(eta+(eps*M)/(1+M))*(1+mu+(rho*M)/(1+M))-b*(q+alp)*(zeta*M-omega0))-b^2*alp*K*(zeta*M-omega0)==0;
[num,den]=numden(lhs(eqn))

num = 

den = 

sigval = linspace(0.01,5,20);
for ii = 1:numel(sigval)
    Msol(ii) = vpa(solve(subs(num,sigma,sigval(ii))));
end
% verify solutions
vpa(subs(lhs(eqn),{M sigma},{Msol sigval})).'

ans = 

plot(sigval,double(Msol),'-x')

[Yakoob]

Thanks Paul for your help!

[Yakoob]

Dear Paul, I just noticed something here. If vpa(Msol1)=2.1265555666349343861231037078585 is a root of the equation above, then when substituting we must get zero. However, I get tis number -8.6736e-19! Could you please clairify this? Many thanks!

[Steven Lord]

Can you prove to me that 2.1265555666349343861231037078585 is the exact root of your system? Are you sure it's not 2.12655556663493438612310370785851 or 2.12655556663493438612310370785849 and is only displayed as 2.1265555666349343861231037078585 because vpa approximated it to 31 decimal places?

Is -8.6736e-19 "close enough" to 0 to be considered zero? For most purposes I'd say yes.

If we look at the Wikipedia page giving some examples of orders of magnitude, the closest example to 8.6e-19 is "The probability of matching 20 numbers for 20 in a game of keno is approximately 2.83 × 10−19." That would require you to choose the same set of 20 numbers from 1-80 as the house. That's not very likely at all.

p = 1./nchoosek(80, 20)

Warning: Result may not be exact. Coefficient has a maximum relative error of 2.2204e-16, corresponding to absolute error 785.

p = 2.8286e-19

[Paul]

Unless solve can find an exact answer, any other solution, whether in double precision or variable precision arithmetic, is going to be an approximation. It's just a matter of how accurate that approximation is and what you're willing to accept. I guess I'm saying the same thing as Steven.