how to find the roots of a symbolic system of 4 coupled algebraic equations with a variable parameter?

Dear Matlab community,

I have the following system of 4 coupled nonlinear algebraic equations:

--------------

syms x y u v

u = 1/J*((E2-Ed+g*(x^2+y^2))*x + gamma2*y/2);

v = 1/J*((E2-Ed+g*(x^2+y^2))*y + gamma2*x/2);

x = 1/J*((E1-Ed+g*(u^2+v^2))*u + gamma1*v/2 +F);

y = 1/J*((E1-Ed+g*(u^2+v^2))*v + gamma1*u/2);

--------------

E1, E2, gamma1, gamma2, g, J, and F are just numbers.

I would like to find x, y, u, v, for variable Ed.

having defined u and v first, x and y stand already as two coupled equations for 2 unknowns, from which I could deduce u and v later. My first thought was to make a polynomial for either x or y, find the roots, and plug in the result in the other equation, and so on.

However, I get an error message when trying to use sym2poly: Input has more than one symbolic variable.

I also tried using the function solve with x and y in there, but that didnt seem to work either.

Your help will be much appreciated.

Answer 1

[Walter Roberson]

syms x y u v Ed

Eq1 = u == 1/J*((E2-Ed+g*(x^2+y^2))*x + gamma2*y/2);

Eq2 = v == 1/J*((E2-Ed+g*(x^2+y^2))*y + gamma2*x/2);

Eq3 = x == 1/J*((E1-Ed+g*(u^2+v^2))*u + gamma1*v/2 +F);

Eq4 = y == 1/J*((E1-Ed+g*(u^2+v^2))*v + gamma1*u/2);

solve([Eq1,Eq2,Eq3,Eq4],[x,y,u,v])

The approach you are considering will not work in the form stated. Your x and y look like simple coupled equations in two unknowns, but when you examine the u and v you find that they are each defined in terms of x and y, so x and y are really being defined partly in terms of each other, not just in terms of two outside simple unknowns. You have four coupled equations in four unknowns, not two in two unknowns.

The closest you might get with your approach would be to initialize x and y, use that to calculate u and v, use those to refine x and y, and keep looping around until the system had settled down to particular values. But if I recall correctly, coupled equations of the form you show can be chaotic, so you might never settle to particular values to find the stationary point.

You could try to use stepwise elimination. Solve Eq1 for x, substitute that x into Eq2, solve that for y, solve the x and y into Eq3, solve for u, and so on. However, the expression for y is longish and involves a number of copies of RootOf() of a messy expression that has u and v buried in it. At that point you are effectively stuck: solving for u in an expression that complicated is not going to work. Opps, you also run into the difficulty that there are three solutions for solving Eq1 for x.

To express a point specifically: sym2poly() is strictly for polynomials, and you are not working with polynomials: notice you will have terms that include x to a power multiplied by y to a power. If you follow the Tips section for sym2poly() it points to coeffs() to allow you to extract symbolic coefficients. But that isn't going to get you much further. Probably you want to look at matlabFunction() and work iteratively -- though as I noted above it isn't certain that you would be able to find the stationary point that way.

Also: keep in mind when you are trying to find the solution to multiple simultaneous nonlinear equations that the precision of calculation of each of the floating point quantities becomes important. Your E1, E2, gamma1, gamma2, g, J, and F are floating point numbers, and when floating point numbers interact with symbolic numbers you are almost certain to run into unintended effects. For stability you should convert each of the values to symbolic form, preferably rational, before using them in the symbolic equations. For example,

syms rJ rE2 rg rgamma2
rJ = sym(J,'f');
rE2 = sym(E2,'f');
rg = sym(E2,'f');
rgamma2 = sym(gamma2,'f');

Eq1 = u == 1/rJ*((rE2-Ed+rg*(x^2+y^2))*x + rgamma2*y/2);

[Said Rodriguez]

Hi Walter. Thank you very much for your quick response. You mention several problems I have but unfortunately I couldnt extract a solution from your response.

A bit more details about the problem: the equations describe the stationary state of a nonlinear system. Its dynamics can be chaotic, but I am looking at stationary state(s) only. The issue here is that for variable Ed, or variable F, (driving frequency or amplitude) the system displays multistability: 2 or more stable outputs for a single input. I believe that an iterative method starting from 1 guess is not a good idea because I will only arrive to one of several solutions, and I need to get all solutions. If I am able to cast the equations in the form of a polynomial, I could find the roots and in this way get all solutions; or so I thought. That is why i was trying to go that way.

I tried using the coeffs function, but that didnt work either. I am not sure I understand how to work it out with matlabFunction(). Could you please give a more specific advice?

Below I include the new code following your advice from previous message:

----------------------

syms x y u v

% these are just numbers

E1 = 0; E2 = 0; gamma1 = 0.01; gamma2 = gamma1; g=0.1*gamma1; J = 3*gamma1; F=gamma1; Ed = 0;

% convert numbers to symbolic form

syms rE1 rE2 rgamma1 rgamma2 rg rJ rF rEd

rE1 = sym(E1,'f'); rE2 = sym(E2,'f'); rgamma1 = sym(gamma1,'f'); rgamma2 = sym(gamma2,'f'); rg = sym(g,'f'); rJ = sym(J,'f'); rF = sym(F,'f'); rEd = sym(Ed,'f');

% here the equations:

Eq1 = u == 1/rJ*((rE2-rEd+rg*(x^2+y^2))*x + rgamma2*y/2);

Eq2 = v == 1/rJ*((rE2-rEd+rg*(x^2+y^2))*y + rgamma2*x/2);

Eq3 = x == 1/rJ*((rE1-rEd+rg*(u^2+v^2))*u + rgamma1*v/2 +rF);

Eq4 = y == 1/rJ*((rE1-rEd+rg*(u^2+v^2))*v + rgamma1*u/2);

solve([Eq1,Eq2,Eq3,Eq4],[x,y,u,v])

-----

When running the above script, I get '[ empty sym ]'

Thanks a lot for you help, and I will highly appreciate any insight you may have.

[Walter Roberson]

empty sym as a result means that as far as it could tell there is no solution.

When I threw some random-ish numbers in for the coefficients I was able to demonstrate that there was no real-valued solution for those particular coefficients. I will have a look later with the values you have given here.

It would help if you could specify whether any of the values are permitted to be imaginary valued?

[Said Rodriguez]

Thanks again. I am only interested in real solutions for the variables u, v, x, and y. Actually, these four variables correspond to the real and imaginary parts of 2 coupled equations for complex variables A and B.

The original equations are ,

B = ((E1-Ed-i*gamma1/2+g*A*conj(A))*A + F)/J;

A = ((E2-Ed-i*gamma2/2+g*B*conj(B))*B)/J;

I need to solve for A and B, which are in general complex. To arrive to the 4 equations at the beginning of the post, I wrote A = u+i*v and B = x+i*y. Then I equated real parts for A and B (2 equations), and then I equated imaginary parts for A and B (2 equations). This is how I ended with 4 equations and 4 real unknowns, namely u, v, x and y. My plan was to write these 4 equations as a single polynomial for any one of those variables by successive substitution. However, that proved to be impossible for me so far. The complication in getting a polynomial for A and B directly is in the A*conj(A)*A and B*conj(B)*B terms. I think I can not use such functions as a coefficient in a polynomial whose roots I want to find, can I?

Now that I explained the underlying problem a bit more, perhaps you can think of another way to solve this which is not along the lines discussed in previos posts?

Thank you!