Matlab gives wrong solution of a nonlinear system of equations.

I had bigger (6 equations) system of ODE and reduced it to the follwing system to check if we cand find equilibruim points .

The code:

syms x y1 y2 real;
assume (x>=0,y1>=0,y2>=0);
syms m1 m2 epsilon1 epsilon2 positive;
syms omega11 omega12 omega 21 omefa22 positive;
syms b12 p12 gamma12 beta1 beta2 positive;
F2=x*(1-x-y1/(m1*y1+x)-y2/(m2*y2+x))==0;
F6=y1*(-epsilon1*(1+omega11*y1+omega12*y2)-(b12*y2)/(p12*y2+y1)+(beta1*x)/(m1*y1+x))==0;            
F7=y2*(-epsilon2*(1+omega21*y1+omega22*y2)+(gamma12*y1)/(p12*y2+y1)+(beta2*x)/(m2*y2+x))==0;
eqns=[F2,F6,F7];

diary('SolutionOfEquation_Subsystem_NonDim_RDFR_Condition.txt')
S = solve(eqns,[x,y1,y2], 'returnconditions', true);
x=latex(S.x)
y1=latex(S.y1)
y2=latex(S.y2)
S.x
S.y1
S.y2

However Matlab gave one solution in term of z and others in terms of parameters ,what does that mean and how I can use the results? Is (z,z1,z2) a free solution?

Is it possible to get z's in terms of parameters ?

The solution is in the attached file.

I also tried to check if the solution is satisfies the equations but find out that just the first and second do and others no.

[darova]

Use ode45 or dsolve

[F.O]

@darova I am solving the steady state system not ODE system here. Analytical solution using dsolve is not possible because I tried.

[Walter Roberson]

The poster is not looking for results for any one set of specific numeric values for the 13 parameters: the poster is looking for the symbolic solution in 13 parameters and is looking for all of the solutions, not just one.

Unfortunately, some of the solution branches require solving polynomial of degree 5 and other branches require solving polynomial of degree 7.

It would perhaps be acceptable to state some of the variables in terms of the roots of high degree polynomials, but the poster also has a need to characterize all the conditions under which the solutions are valid. That is not exactly the same as solving the polynomial, but it does require characterizing the conditions under which the roots of the polynomial are real valued which appears to be approachable (with some work) for numeric coefficients but not so much for symbolic coefficients https://mathoverflow.net/questions/20946/criteria-to-determine-whether-a-real-coefficient-polynomial-has-real-root (symbolic eigenvalues are not realistically calculable for 5x5 or larger especially for coefficients as complicated as would be required here)

Answer 1

[Walter Roberson]

The answer to the question is, NO, it is not possible to get z in terms of the parameters.

... At least not if, like me, you only have 32 gigabytes of memory.

Your system generates many different branches of conditional cases. Working all of the cases through takes a lot of time to figure out how to proceed, and takes even more memory. I ended up having to kill the process after it started swapping my system.

I suspect you just might be able to get further if you were to remove some of the constraints to positive, as that would reduce the number of different cases to consider. You would then weed through through the results after you had the general framework.

[F.O]

When I remove the conditions on variables and parameters , I dont get answer at all . I get error and numerical values in the solution.

The solution when removing the conditions is in the attached file.

The error message

Warning: Solutions might be lost. 
> In symengine
  In mupadengine/evalin_internal (line 104)
  In mupadengine/feval_internal (line 167)
  In solve (line 301)
  In Analysis_Subsystem_NonDim_RDFR_ZeroInflow_NoCond (line 17) 
Warning: Unable to solve symbolically. Returning a numeric solution using vpasolve. 
> In solve (line 304)
  In Analysis_Subsystem_NonDim_RDFR_ZeroInflow_NoCond (line 17) 

[Walter Roberson]

One of the simpler solutions is

  x = 0
  
  y1 = 1/4*(2^(1/2)*signum(omega11*omega22-omega12*omega21)*signum(-gamma12*omega11*
p12^2*epsilon1+b12*omega21*p12*epsilon2+gamma12*omega12*p12*epsilon1-b12*
omega22*epsilon2)*((-gamma12*omega11*p12*epsilon1+b12*omega21*epsilon2)*signum(
-gamma12*omega11*p12*epsilon1+b12*omega21*epsilon2)*(b12*((omega21*p12-2*
omega22)*omega11+omega12*omega21)*epsilon2-p12*gamma12*omega11*epsilon1*(
omega11*p12-omega12))*((-(omega11*p12-omega12)*(omega11*p12-omega21*p12-omega12
+omega22)*epsilon1+b12*((omega21*p12-2*omega22)*omega11+omega12*omega21))*
epsilon2-gamma12*omega12*epsilon1*(omega11*p12-omega12))*((((omega11-omega21)*
p12-omega12+omega22)^2*epsilon2^2+2*gamma12*(((omega11+omega21)*omega12-2*
omega11*omega22)*p12-omega12*(omega12-omega22))*epsilon2+gamma12^2*omega12^2)*
epsilon1^2-2*b12*((omega21*(omega11-omega21)*p12+omega12*omega21-2*omega22*(
omega11-1/2*omega21))*epsilon2+2*gamma12*omega11*omega22-gamma12*omega21*
omega12)*epsilon2*epsilon1+b12^2*omega21^2*epsilon2^2)^(1/2)+b12^2*((omega11*
p12-omega21*p12-omega12+omega22)^2*((omega21^2*p12^2-2*omega21*omega22*p12+2*
omega22^2)*omega11^2-2*omega11*omega12*omega21*omega22+omega12^2*omega21^2)*
epsilon1^2-2*b12*((omega21*p12-2*omega22)*omega11+omega12*omega21)*((omega21^2*
p12^2-2*omega21*omega22*p12+2*omega22^2)*omega11^2-omega21*(omega21^2*p12^2-2*
omega21*omega22*p12+2*omega12*omega22+omega22^2)*omega11+omega12^2*omega21^2)*
epsilon1+b12^2*((omega21*p12-2*omega22)*omega11+omega12*omega21)^2*omega21^2)*
epsilon2^4-2*epsilon1*b12*(p12*omega11^2*(omega21*p12-omega22)*(omega11*p12-
omega12)*(omega11*p12-omega21*p12-omega12+omega22)^2*epsilon1^2-2*b12*(p12^2*(
omega21*p12-2*omega22)*(omega21*p12-omega22)*omega11^4-(p12^3*omega21^3-1/2*
omega21^2*(omega12+5*omega22)*p12^2-p12*omega12*omega21*omega22+omega22^2*(
omega12+2*omega22))*p12*omega11^3-3/2*((omega12+7/3*omega22)*omega21^2*p12^2-2/
3*omega21*omega22*(omega12+4*omega22)*p12+2/3*omega22^2*(omega12-omega22))*
omega12*omega11^2+2*omega12^2*omega21*(p12^2*omega21^2+1/4*omega21*(omega12-5*
omega22)*p12+1/2*omega22*(omega12-omega22))*omega11-1/2*omega21^2*omega12^3*(
omega12-omega22))*epsilon1+b12^2*((omega21*p12-2*omega22)*omega11+omega12*
omega21)*((omega21^2*p12^2-2*omega22^2)*omega11^2+2*omega11*omega12*omega21*
omega22-omega12^2*omega21^2))*gamma12*epsilon2^3+epsilon1^2*gamma12^2*(p12^2*
omega11^2*(omega11*p12-omega12)^2*(omega11*p12-omega21*p12-omega12+omega22)^2*
epsilon1^2-2*b12*(omega11*p12-omega12)*(p12^2*(omega21*p12-2*omega22)*omega11^3
+(-omega21^2*p12^3+2*omega12*omega21*p12^2+3*omega22^2*p12)*omega11^2-3*omega12
*(omega21*(omega12+2/3*omega22)*p12-2/3*omega22*(omega12-omega22))*omega11+2*
p12*omega12^2*omega21^2)*p12*omega11*epsilon1+b12^2*(p12^2*(omega21^2*p12^2+2*
omega21*omega22*p12-6*omega22^2)*omega11^4-2*p12*omega12*omega22*(omega21*p12-4
*omega22)*omega11^3+2*omega12^2*(omega21^2*p12^2-3*omega21*omega22*p12+omega22^
2)*omega11^2-2*omega21*omega22*omega12^3*omega11+omega12^4*omega21^2))*epsilon2
^2-2*epsilon1^3*(omega11*p12-omega12)*gamma12^3*p12*(-p12*(omega11*p12-omega12)
^2*(omega12-omega22)*epsilon1+b12*(p12^2*omega11^2*omega22-2*p12*omega12*
omega22*omega11+2*(p12*omega21-1/2*omega22)*omega12^2))*omega11^2*epsilon2+p12^
2*gamma12^4*omega11^2*omega12^2*epsilon1^4*(omega11*p12-omega12)^2)^(1/2)-
signum(-gamma12*omega11*p12*epsilon1+b12*omega21*epsilon2)*(omega11*p12+omega12
)*(-gamma12*omega11*p12*epsilon1+b12*omega21*epsilon2)*((((omega11-omega21)*p12
-omega12+omega22)^2*epsilon2^2+2*gamma12*(((omega11+omega21)*omega12-2*omega11*
omega22)*p12-omega12*(omega12-omega22))*epsilon2+gamma12^2*omega12^2)*epsilon1^
2-2*b12*((omega21*(omega11-omega21)*p12+omega12*omega21-2*omega22*(omega11-1/2*
omega21))*epsilon2+2*gamma12*omega11*omega22-gamma12*omega21*omega12)*epsilon2*
epsilon1+b12^2*omega21^2*epsilon2^2)^(1/2)-((p12^2*omega11^2+(-omega21*p12^2-3*
omega22*p12)*omega11+3*omega12*(p12*omega21-1/3*omega12+1/3*omega22))*epsilon2+
gamma12*omega12*(omega11*p12+omega12))*gamma12*p12*omega11*epsilon1^2+((p12*(
omega21*p12-2*omega22)*omega11^2-(omega21*p12-omega22)*(omega21*p12-2*omega12+2
*omega22)*omega11+omega21*omega12*(omega21*p12+omega12-omega22))*epsilon2+(
omega11*p12+omega12)*gamma12*((omega21*p12+2*omega22)*omega11-omega12*omega21))
*b12*epsilon2*epsilon1-b12^2*omega21^2*epsilon2^2*(omega11*p12+omega12))/
epsilon1/epsilon2/(omega11*omega22-omega12*omega21)/(-p12*gamma12*(omega11*p12-
omega12)*epsilon1+b12*epsilon2*(omega21*p12-omega22))/omega11

y2 = 1/2*(signum(-gamma12*omega11*p12*epsilon1+b12*omega21*epsilon2)*(-gamma12*
omega11*p12*epsilon1+b12*omega21*epsilon2)*((((omega11-omega21)*p12-omega12+
omega22)^2*epsilon2^2+2*(((omega11+omega21)*omega12-2*omega11*omega22)*p12-
omega12*(omega12-omega22))*gamma12*epsilon2+gamma12^2*omega12^2)*epsilon1^2-2*
b12*((omega21*(omega11-omega21)*p12+omega12*omega21-2*omega22*(omega11-1/2*
omega21))*epsilon2+2*gamma12*omega11*omega22-gamma12*omega21*omega12)*epsilon2*
epsilon1+b12^2*omega21^2*epsilon2^2)^(1/2)+gamma12*((omega11*(omega11-omega21)*
p12+(-omega12-omega22)*omega11+2*omega12*omega21)*epsilon2+gamma12*omega11*
omega12)*p12*epsilon1^2-b12*((omega21*(omega11-omega21)*p12-2*omega11*omega22+
omega21*(omega12+omega22))*epsilon2+gamma12*(omega11*omega21*p12+2*omega11*
omega22-omega12*omega21))*epsilon2*epsilon1+b12^2*omega21^2*epsilon2^2)/(-p12*
gamma12*(omega11*p12-omega12)*epsilon1+b12*epsilon2*(omega21*p12-omega22))/
epsilon1/epsilon2/(omega11*omega22-omega12*omega21)

This then needs to be analyzed to produce the full set of conditions under which y1 and y2 are non-negative. That kind of analysis is necessary because one small step away from the above solution is another solution involving sqrt(-1) times an expression raised to 1/2, and you cannot rule out the solution because there are conditions under which the expression can potentially be negative, giving another sqrt(-1) from the ^(1/2) and the two sqrt(-1) would then cancel out potentially leading to a valid real-valued solution. Working through all the various co-relations between variables to give a real-valued solution for the case is tedious and not realistically workable using the MATLAB Symbolic Toolbox.

And then you get to one of the other solution branches, and you soon get situations where the solutions involve roots of an equation of degree 5 or 7, and it is not workable to figure out the combinations of parameters that lead to real-valued solutions for those.

As I told you before, you have too many parameters to come out with a solution that is going to be in any sense understandable.

[F.O]

I agree that parameters make the analysis difficult but som authors analysed similiar models and got solution.

I already reduced my system from 7 to 3 but still dffcult to do any analysis when choosng Arditi-Ginzberg functional response.

I just wonder how the author in the follwing paper found solution to his system because it is not very different from the one I work wirh.

Take a look at this

diffusive one-prey and two-competing-predator system with a ratio-dependent functional response- I, long time behavior and stability of equilibria

[Walter Roberson]

How many parameters does the author have to work with in that paper?

[F.O]

10 after nondimensionallizing

[F.O]

I tried to replace (b12*y2)/(p12*y2+y1) by just b12*y2 and (gamma12*y1)/(p12*y2+y1) by gamma12*y1 i.e using linear functional response but still got parametrized solution in case of imposing condition and error in case not imposing any condition.

Answer 2

[James Browne]

Looks like you have quite a bit going on there. I am not very familiar with doing symbolic maths with MATLAB, but I can tell you that you have a syntax error in the first line of your code.

Here is a link to the help page for the syms function:

https://www.mathworks.com/help/symbolic/syms.html

Note that the input variables are separated by spaces in all of the examples, not commas. I am not really sure what this is causing MATLAB to do, but it does not seem to be working correctly and you can tell this by looking at your input arguments on line 1. They should all be highlighted in purple, like the examples on the help page, but only the variable "x" is highlighted, along with your variable type selection.

I would try deleting your commas and use spaces to separate your input aruments in the syms function on line 1.

[F.O]

it is a typo here sorry. In my atual code I don't have comma. I will edit the question.