A positive root of an equation

Hello dear,

I have a mathemtaical model and one of the dependant variable, M, is a positive root of a nonlinear equation. I am not interested in finding this M but I would like to tell MATLAB that M is a positive root of the equation r(M)(a+P(M))(K-d)-bKM=0. How to do that?

Thank you!

Answer 1

[Torsten]

Without further knowledge about the background of your question:

Maybe by setting (r(M)(a+P(M))(K-d))/(bK) instead of M in your model.

[Fares]

I set M to be equal to (r(M)(a+P(M))(K-d))/(bK) but still MATLAB doesnt recognize M. I am wondering if I can do something close to this:

First I introduce the equation as

F(M)=r(M)(a+P(M))(K-d)-bKM;

and then I write something like this:

M=roots(F(M),M>0)

Any suggestions can be helpful! Thanks!

[Torsten]

You said you were not interested in finding the root. And now you solve for it ?

To give advice, we need more information about the background of your question.

If you want to solve an equation for a variable, usually "vpasolve" or "solve" is the correct command.

If the equation is a polynomial in the unknown, you can use "roots".

To force a variable to be positive, you may use

assume(M,'positive')

or directly when you declare the variable

syms M positive

I assume you are talking about symbolic variables, don't you ?

[Fares]

Thank for your response, Torsten, and I apologize for any inconvenience caused.

I have this equation

(r0*(1+k1*(mu+(rho*M)/(1+M))*(1-k2*(mu+(rho*M)/(1+M)))))*(eta+(epsilon*M)/(1+M))*(K-alpha*sigma*m*(eta+(epsilon*M)/(1+M))*(1+(mu+(rho*M)/(1+M)))-b*(m+alpha)*(pi*M-omega))-b^2*alpha*K*(pi*M-omega) == 0

and I would like to write M explicitly.

[Torsten]

You have a polynomial of degree 6 in M. For polynomials with degree >4, it's not possible to give an analytical expression for their roots. You will have to give values for the parameters and use "roots" to solve for M numerically.

syms r0 k1 mu rho M k2 eta epsilon alpha sigma m eta b K omega
expr = r0*(1+k1*(mu+rho*M/(1+M))*(1-k2*(mu+rho*M/(1+M))))*(eta+epsilon*M/(1+M))*(K-alpha*sigma*m*(eta+epsilon*M/(1+M))*(1+mu+rho*M/(1+M))-b*(m+alpha)*(pi*M-omega))-b^2*alpha*K*(pi*M-omega) 

expr = 

[N,D] = numden(expr)

N = 

D = 

r = coeffs(N,M)

r = 

size(r)

ans = 1×2

1 7

[Fares]

Thanks Torsten for your help. I do not think I need to include all of this in my research. The idea here is that suppose we have this mathematical model (an artificial model):

dx/dt=ax-bxy

dy/dt=bxy-cy

dz/dt=dy-mz

dM/dt=e+ky-wM

Suppose after finding the positive fixed point we ended up with expressions like these:

x=cm/d*bM

y=c/b*(wM-e)

z=(wM-e)/b

but M is a positive root of some polynimoal of degree > 4, let's say the polynomial takes this form:

(r0*(1+k1*(mu+(rho*M)/(1+M))*(1-k2*(mu+(rho*M)/(1+M)))))*(eta+(epsilon*M)/(1+M))*(K-alpha*sigma*m*(eta+(epsilon*M)/(1+M))*(1+(mu+(rho*M)/(1+M)))-b*(m+alpha)*(pi*M-omega))-b^2*alpha*K*(pi*M-omega) == 0

Now, I know how to define the first three positive fixed points in MATLAB but I have no idea how to tell MATLAB that the last positive fixed point is a positive root of this polynomial. Is there a way to do this?

I really appreciate your help.

[Torsten]

Now, I know how to define the first three positive fixed points in MATLAB but I have no idea how to tell MATLAB that the last positive fixed point is a positive root of this polynomial.

For which purpose do you want to use this fact within MATLAB ? Do you need M in further calculations ?

[Fares]

Thanks Torsten for your reply.

I am trying to calculate numerically the stability boundaries using the code I have. In the code, you have to insert the poitive fixed point. I am getting stuck in this part as I do not know how to tell MATLAB that M cannot be written explicitly, it is just a poisitve root of some polynomial.

[Torsten]

And the stability boundaries are calculated symbolically ?

Otherwise why don't you calculate the positive root and hand it to the code that calculates the stability boundary ?

format long
r = [-1 2 3 16 5 6 7];
sol = roots(r)

sol = 
  3.921379649150494 + 0.000000000000000i
 -0.858569973575940 + 1.771560634548524i
 -0.858569973575940 - 1.771560634548524i
  0.248532487550623 + 0.771366979108634i
  0.248532487550623 - 0.771366979108634i
 -0.701304677099852 + 0.000000000000000i

sol = sol(real(sol)>0 & abs(imag(sol))<1e-6)

sol = 
   3.921379649150494

[Fares]

Thanks Torsten for your continuous support.

To study the stability boundary, we fix all parameters except two. Varying these two might tell us the stability regions of our fixed point. The problem is that the polynomial has one of these two parameters. This parameter has to be there symbolically.

[Torsten]

If you don't get passed numerical values for all the parameters of this polynomial that enable you to return a numerical root, you won't succeed.

But somewhere in the stability boundary code, the additional symbolic parameter must become numeric in order to compute the boundary, doesn't it ?

[Fares]

Thank you for your cooperation, Torsten!

The only thing mentions about the symbolic parameter, sigma, is this:

NN2=20;

sigma=linspace(0.01,5,NN2);

Does it become numeric with this?

[Torsten]

Yes, it's numeric (double). So calculating numerical roots for the polynomial should work, shouldn't it ?

NN2=20;
sigma=linspace(0.01,5,NN2)

sigma = 1×20

0.0100 0.2726 0.5353 0.7979 1.0605 1.3232 1.5858 1.8484 2.1111 2.3737 2.6363 2.8989 3.1616 3.4242 3.6868 3.9495 4.2121 4.4747 4.7374 5.0000

class(sigma)

ans = 'double'

[Fares]

Thank you for your help Torsten!

So this is what I have now:

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;

m=0.5;

b=0.1;

pi=4;

omega0=1;

sigma=linspace(0.01,5,20);

syms M

assume(M,'positive')

eqn = r0*(1+k1*(mu+(rho*M)/(1+M))*(1-k2*(mu+(rho*M)/(1+M))))*(eta+(eps*M)/(1+M))*(K-alp*sigma*m*(eta+(eps*M)/(1+M))*(1+mu+(rho*M)/(1+M))-b*(m+alp)*(pi*M-omega0))-b^2*alp*K*(pi*M-omega0)==0;

Would it be possible to have all positive real roots of this equation?

[Torsten]

format long
syms M sigma
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;
m=0.5;
b=0.1;
pi=4;
omega0=1;
expr = r0*(1+k1*(mu+(rho*M)/(1+M))*(1-k2*(mu+(rho*M)/(1+M))))*(eta+(eps*M)/(1+M))*(K-alp*sigma*m*(eta+(eps*M)/(1+M))*(1+mu+(rho*M)/(1+M))-b*(m+alp)*(pi*M-omega0))-b^2*alp*K*(pi*M-omega0);
[N,D] = numden(expr);
r = coeffs(N,M);

sigma_num = linspace(0.01,5,20);
for i=1:numel(sigma_num)
  ri = double(roots(subs(r,sigma,sigma_num(i))));
  r_num(i) = ri(real(ri)>0 & abs(imag(ri))<1e-6);
end
plot(sigma_num,r_num)

[Fares]

Thanks Torsten, this is very helpful. In case I want only one positive real root of the equation, what command should I use?

[Torsten]

Did you have a case where you got more than one solution ?

The above code gave only one positive root for each value of sigma.

[Torsten]

There seems to be more than one root that satisfies that it is real and positive. And which one should be taken ? The smallest ? The largest ? The one that is nearest to pi ? :-)

[Fares]

According to my analytical work, M is just a positive real root. So I think it does not matter the location of M as long as it is real and positive.

However, when I solve the mathmatical model with these selection of parameters, M always tends to a value that is very near to omega0/pi but greater than omega0/pi. I would say we take the real positive root nearest to omega0/pi but greater than omega0/pi. Many thanks!

[Torsten]

I just saw that you didn't copy my code correctly.

It must read

[N,D] = numden(expr);

instead of

[H,D] = numden(expr);

Does it work out properly now or are there still changes needed ?

[Fares]

Dear Torsten,

Thank you for your help.

I just noticed something in the code you posted here in one of the comments, which is

format long
syms M sigma
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;
m=0.5;
b=0.1;
pi=4;
omega0=1;
expr = r0*(1+k1*(mu+(rho*M)/(1+M))*(1-k2*(mu+(rho*M)/(1+M))))*(eta+(eps*M)/(1+M))*(K-alp*sigma*m*(eta+(eps*M)/(1+M))*(1+mu+(rho*M)/(1+M))-b*(m+alp)*(pi*M-omega0))-b^2*alp*K*(pi*M-omega0);
[N,D] = numden(expr);
r = coeffs(N,M);

sigma_num = linspace(0.01,5,20);
for i=1:numel(sigma_num)
  ri = double(roots(fliplr(subs(r,sigma,sigma_num(i)))));
  r_num(i) = ri(real(ri)>0 & abs(imag(ri))<1e-6);
  res(i) = double(subs(expr,[sigma M],[sigma_num(i) r_num(i)]));
end
res.'

ans = 20×1

1.0e-17 * -0.116957647785477 0.029945010337059 0.107969523137403 -0.005425686328019 0.071367757116093 -0.009157905701059 -0.048611801937400 -0.118476961985516 0.114341480692949 -0.068606849757453

plot(sigma_num,r_num)

The thing is that substituting sigma_num(1) and r_num(1) in the equation

eqq = r0*(1+k1*(mu+(rho*r_num(1))/(1+r_num(1)))*(1-k2*(mu+(rho*r_num(1))/(1+r_num(1)))))*(eta+(eps*r_num(1))/(1+r_num(1)))*(K-alp*sigma_num(1)*m*(eta+(eps*r_num(1))/(1+r_num(1)))*(1+mu+(rho*r_num(1))/(1+r_num(1)))-b*(m+alp)*(pi*r_num(1)-omega0))-b^2*alp*K*(pi*r_num(1)-omega0)

should give us zero since r_num(1) is a root of the equation. However I get a different value

eqq =

0.006632993397833

I am not sure what is wrong here! Could you please clairify this?

[Torsten]

I changed the code above so that it now gives correct results.

I didn't notice that "coeffs" gives the vector of polynomial coefficients in the reverse order as "roots" expects them. Therefore the "fliplr" command had to be added.

Sorry for that.

[Fares]

Thanks Torsten for your replay.

No need to be sorry. Many thanks for your continuous support!

Answer 2

[John D'Errico]

syms M r0 epsilon rho k1 K b m alpha omega mu k2 eta sigma
P = (r0*(1+k1*(mu+(rho*M)/(1+M))*(1-k2*(mu+(rho*M)/(1+M)))))*(eta+(epsilon*M)/(1+M))*(K-alpha*sigma*m*(eta+(epsilon*M)/(1+M))*(1+(mu+(rho*M)/(1+M)))-b*(m+alpha)*(pi*M-omega))-b^2*alpha*K*(pi*M-omega) == 0

P = 

If you multiply by (M+1)^2, then regardless if I assume that m is not just a typo where you intended to write M in one place, this would appear to be effectively a 6th degree polynomial in M as long as M~=1.

If all of those parameters have known values, then you can just use vpasolve to return the 6 roots. If you want to find a solution in terms of symbolc parameters all as unknowns, then you are wasting your time. Abel-Ruffini (long before any of us was born) showed that general polynomials of degree higher than 4 will have no solutions in algebraic form. (Except for certain trivial cases.)

Of course, if you wish, you always can try this:

Msol = solve(P,M,'maxdegree',6,'returnconditions',true)

Msol = struct with fields:
             M: [6×1 sym]
    parameters: [1×0 sym]
    conditions: [6×1 sym]

Msol.M

ans = 

Which is MATLAB's way of telling you that as much as it wants to solve the problem, it cannot resolve the issue of mathematical impossibility.

[Fares]

Thanks John for your help!