Solving a system of equations with arbitrary Parameter value

I have the following system of equations:

syms x0 y0 omega0 E0 Pm0
% x0, y0 omega0 E0 are variables
% Pm0 is a parameter
assume([x0 y0 omega0 E0], 'real');
assumeAlso(x0*x0 + y0*y0 == 1);

% I want to solve the following system of equations for my variables
% parametrized by Pm0. Then substitue different values for Pm0 to get
% solutions for each Pm0 value.
eqn1 = -y0*omega0 + x0*(1-x0*x0-y0*y0) == 0;
eqn2 = x0*omega0 + y0*(1-x0*x0-y0*y0) == 0;
eqn3 = -3.33*E0*y0 - 0.66*omega0 + 3.33 * Pm0 == 0;
eqn4 = 0.5*x0 - E0 + 0.5 == 0;

% NOTE: I need to calculate the solns for thousands of parameter values.

What I have tried:

% create a vector of the equations to feed into solve();
eqn_final = [eqn1, eqn2, eqn3, eqn4];
% create a vector of variables to solve for
variable_final = [x0, y0, omega0, E0];

solution = solve(eqn_final, variable_final, "ReturnConditions", true);

The Problem:

This gives me a structure "solution". which has parameters and conditions.

How do I interpret these and finally substitute different values for Pm0 to get solutions to the above system of equations for different values of the parameter?

Answer 1

[John D'Errico]

Look carefullly at what is returned.

solution = solve(eqn_final, variable_final, "ReturnConditions", true);

solution
solution = 
  struct with fields:

            x0: [3×1 sym]
            y0: [3×1 sym]
        omega0: [3×1 sym]
            E0: [3×1 sym]
    parameters: [1×1 sym]
    conditions: [3×1 sym]

So there are THREE distinct solutions possible that it found.

solution.x0
ans =
     -1
      1
2*x - 1

So x0 is one of the unknowns. In the first two solutions, it is a fixed value.

Looking at the conditions return, we see why.

solution.conditions
ans =
                                                                               Pm0 == 0
                                                                               Pm0 == 0
in(x, 'real') & Pm0 ~= 0 & Pm0^2 + 4*x^4 == 4*x^3 & (in(Pm0, 'real') | x == 0 | x == 1)

Ah, when Pm0 is exactly zero, the problem reduces. But MATLAB does not know that Pm0 will not be zero. and when Pm0 is zero, the problem becomes considerably simpler.

So what happens when pm0 is some non-zero value? Anything but zero. Then we have solution number 3, the one you will care about.

So all seems simple. Just pick solution #2, and substiture the value of Pm0. But is it simple? Does your problem have a solution?

solution.x0(3)
ans =
2*x - 1

solution.y0(3)
ans =
(4*(- x^3 + x^2))/Pm0

solution.omega0(3)
ans =
0

solution.E0(3)
ans =
x

This tells us, that for ANY value of Pm0, there are infinitely many solutions. We can arbitrarily choose some real number for x. What is x? x just happens to be the value of E0. It seems you can arbitrarily choose the value for E0.

Could we see that by some simple manipulations of those equations? It gets complicated, because when we try this:

simplify(eqn2 - eqn1)
ans =
x0 + y0 == 0 | omega0 == 0

So that tells us the solution has EITHER omega0 == 0, or it reduces to a simple linear equation in terms of x0 and y0. If omega0 is non-zero, then what happens? Now we have a problem. No non-trivial solution will exist at all, because then we must have x0 == -y0.

That means we can assume omega0 must be zero. Then both of equations 1 and 2 effectively reduce to a unit circle in the parameters x0 and y0, thus the equation

x0^2 + y0^2 == 1

So a circle with center at the origin, and radius 1.

But then we have eqn4, which gives us x0 as a function of E0. And eqn3, which allows us to compute y0, in terms of Pm0 and E0.

subs(eqn3,omega0,0)
ans =
(333*Pm0)/100 - (333*E0*y0)/100 == 0

Effectively, this allows us to solve for y0 as

y0 == Pm0/E0

Similarly, eqn 4 allows us to solve for x0 as

x0 = 2*E0 - 1

But since we already know that x0 and y0 must lie on the unit circle, we must have

subs(x0^2 + y0^2 == 1,[x0,y0],[2*E0-1,Pm0/E0])
ans =
(2*E0 - 1)^2 + Pm0^2/E0^2 == 1

Multiply by E0^2,

simplify(subs(x0^2 + y0^2 == 1,[x0,y0],[2*E0-1,Pm0/E0])*E0^2)
ans =
4*E0^4 + Pm0^2 == 4*E0^3 | E0 == 0

So as long as E0 is non-zer0, that would seem to have 4 solutions, as effectively a 4th degree polynomial in E0. Some of those solutions may be complex of course. And that would depend on the value of Pm0. For some values of Pm0, no solution at all will exist. For example. when Pm0 is 1, we would have the polynomial:

fplot(subs(4*E0^4 + Pm0^2 - 4*E0^3,Pm0,1))

You can see this function NEVER crosses zero. Therefore E0 has no real solutions, and when Pm0 is 1, the problem never has a solution. But when Pm0 has a different value, might the polynomial have a solution?

4*E0^4 + Pm0^2 - 4*E0^3 == 0

As it turns out, since Pm0 is squared there, the constant term in that function is ALWAYS non-negative. The most extreme case is where Pm0 is 0, and then we will have a triple root at E0 == 0, and a single root at E0 == 1.

For some other cases of Pm0, when Pm0 a value between 0 and 1, then we will have exactly 2 real solutions for E0.

solve(subs(4*E0^4 + Pm0^2 - 4*E0^3,Pm0,.5),'maxdegree',4)
ans =
                                                                                                  1/2
1/(9*((11^(1/2)*1728^(1/2))/1728 + 19/216)^(1/3)) + ((11^(1/2)*1728^(1/2))/1728 + 19/216)^(1/3) + 1/6

vpa(ans)
ans =
                               0.5
0.91964337760708056627592628232664

So your system of equations depends on the value of Pm0, but the number of solutions will vary.

[atharva aalok]

I have made the following change:

assume([x0 y0 omega0 E0], 'real');
assumeAlso(x0*x0 + y0*y0 == 1);
% THIS IS NEW: FURTHER ASSUMPTION
assume(Pm0 ~= 0)

Now that problem is gone.

Regarding what you said about there being only 3 independent equations:

I also have the assumption: assumeAlso(x0*x0 + y0*y0 == 1)

Hence I should not get infinitely many solutions.

For example:

syms x0 y0 omega0 E0
assume([x0 y0 omega0 E0], 'real');
% assume(Pm0 ~= 0); I will give Pm0 a definite value...
assumeAlso(x0*x0 + y0*y0 == 1);

Pm0 = 0.4;

eqn1 = -y0*omega0 + x0*(1-x0*x0-y0*y0) == 0;
eqn2 = x0*omega0 + y0*(1-x0*x0-y0*y0) == 0;
eqn3 = -3.33*E0*y0 - 0.66*omega0 + 3.33 * Pm0 == 0;
eqn4 = 0.5*x0 - E0 + 0.5 == 0;

[xFP, yFP, omegaFP, EFP] = solve([eqn1, eqn2, eqn3, eqn4], [x0 y0 omega0 E0]);
disp(double([xFP, yFP, omegaFP, EFP]));

This does give 2 unique solutions:

>> trial
   -0.1858    0.9826         0    0.4071
    0.9078    0.4193         0    0.9539

Also the third condition is:

in(x, 'real') & Pm0 ~= 0 & Pm0^2 + 4*x^4 == 4*x^3 & (in(Pm0, 'real') | x == 0 | x == 1)
% Here it seems to me that Pm0 and x are related by the eqn above. Not
% independent.

So why is the original question I posted giving me infinitely many solutions for any Pm0 value?

Or in other words, how do I make Pm0 also a symbolic variable and when I substitute Pm0 = 0.4 I get the above 2 solutions?

[atharva aalok]

I have added the final version which needs to be overcome to solve the problem as a separate post as its a bit elaborate.

Thank you for all your help

Atharva(Undergraduate in Aerospace Engineering)

[atharva aalok]

I have made some progress and this final problem remains:

I will try to quickly explain the problem in brief again and then pose the problem.

% Creating symbolic variables
% x0 y0 omega0 E0 are physical quantities
% Pm0 is a parameter generally in the range [0 1]
syms x0 y0 omega0 E0 Pm0
assume([x0 y0 omega0 E0], 'real');
assume(Pm0 ~= 0);
% x0 and y0 are related by
assumeAlso(x0*x0 + y0*y0 == 1);

% Equations for the dynamical system, (equating to 0 to find equilibrium
% points)
eqn1 = -y0*omega0 + x0*(1-x0*x0-y0*y0) == 0;
eqn2 = x0*omega0 + y0*(1-x0*x0-y0*y0) == 0;
eqn3 = -3.33*E0*y0 - 0.66*omega0 + 3.33 * Pm0 == 0;
eqn4 = 0.5*x0 - E0 + 0.5 == 0;
 
% solving the equations
eqn_final = [eqn1, eqn2, eqn3, eqn4];
variable_final = [x0, y0, omega0, E0];

a = solve(eqn_final, variable_final, "ReturnConditions", true);

% PmSol is a vector that contains different values for the parameter Pm0
PmSol = linspace(0.2, .5, 100);
% For each value of Pm0 from PmSol we need to solve the equations and find
% equilibrium pts.

% Then we need store the solutions for x0 to E0 for each Pm value.
% create vector of Fixed Point locations for each Pm val in PmSol
Fixed_points_locations = ones(2,length(PmSol)*4)*0;
% For each Pm value we will get either 2 solutions or no solutions
% if we get 2 solutions then we will store a vector like:
% [x0_1, y0_1 omega0_1 E0_1; x0_1, y0_1 omega0_1 E0_1] in above matrix. If
% no solutions then we will store zeros(2 by 4) in above matrix.



% THE PROBLEM:
% THE AIM: the goal is to make the steps from here on out as fast as
% possible. This will become a part of a bigger problem needs to be
% lightning fast!


a

a = struct with fields:
            x0: [1×1 sym]
            y0: [1×1 sym]
        omega0: [1×1 sym]
            E0: [1×1 sym]
    parameters: [1×1 sym]
    conditions: [1×1 sym]

a.parameters

ans = 

x

a.conditions

ans = 

Pm_val = PmSol;
% substitute all Pm values at once in a.conditions to get equation for
% parameter
subs_new = subs(a.conditions, Pm0, Pm_val);  
subs_new(1,1)

ans = 

% Now all we need to do is solve each equation in subs_new (each of which
% looks like the above) get the parameter values, combine that with Pm_val
% and substitue in the expressions for x0 to E0 to get equilibrium points
% for all values of Pm

% THE ISSUE: HOW TO VECTORIZE THIS? (the following for loop)

for j = 1:length(PmSol)
    % solve the each equation in subs_new for the parameter
    parameter_val = solve(subs_new(j));
    % make a vector of same dimensions as parameter vals with value
    % Pm_val(j)
    Pm_new = parameter_val * 0 + Pm_val(j);
    
    % substitute parameter_val and Pm_new to get Equilibrium points
    T1 = {[parameter_val], [Pm_new]};
    T2 = {a.parameters, Pm0};
    xFP = subs(a.x0, T2, T1);
    yFP = subs(a.y0, T2, T1);
    omegaFP = subs(a.omega0, T2, T1);
    EFP = subs(a.E0, T2, T1);
    
    [rows, cols] = size([xFP, omegaFP, EFP]);
    % if no solutions add a zero matrix for the equilibrium point locations
    if rows == 0
        Fixed_points_locations(1:2,3*(j-1)+1:3*(j-1)+3) = ones(2, 3)*0;
        continue;
    end
    Fixed_points_locations(1:2,3*(j-1)+1:3*(j-1)+3) = [xFP, omegaFP, EFP];

end

Answer 2

[atharva aalok]

I have made some progress and this final problem remains:

I will try to quickly explain the problem in brief again and then pose the problem.

% Creating symbolic variables
% x0 y0 omega0 E0 are physical quantities
% Pm0 is a parameter generally in the range [0 1]
syms x0 y0 omega0 E0 Pm0
assume([x0 y0 omega0 E0], 'real');
assume(Pm0 ~= 0);
% x0 and y0 are related by
assumeAlso(x0*x0 + y0*y0 == 1);

% Equations for the dynamical system, (equating to 0 to find equilibrium
% points)
eqn1 = -y0*omega0 + x0*(1-x0*x0-y0*y0) == 0;
eqn2 = x0*omega0 + y0*(1-x0*x0-y0*y0) == 0;
eqn3 = -3.33*E0*y0 - 0.66*omega0 + 3.33 * Pm0 == 0;
eqn4 = 0.5*x0 - E0 + 0.5 == 0;
 
% solving the equations
eqn_final = [eqn1, eqn2, eqn3, eqn4];
variable_final = [x0, y0, omega0, E0];

a = solve(eqn_final, variable_final, "ReturnConditions", true);

% PmSol is a vector that contains different values for the parameter Pm0
PmSol = linspace(0.2, .5, 100);
% For each value of Pm0 from PmSol we need to solve the equations and find
% equilibrium pts.

% Then we need store the solutions for x0 to E0 for each Pm value.
% create vector of Fixed Point locations for each Pm val in PmSol
Fixed_points_locations = ones(2,length(PmSol)*4)*0;
% For each Pm value we will get either 2 solutions or no solutions
% if we get 2 solutions then we will store a vector like:
% [x0_1, y0_1 omega0_1 E0_1; x0_1, y0_1 omega0_1 E0_1] in above matrix. If
% no solutions then we will store zeros(2 by 4) in above matrix.



% THE PROBLEM:
% THE AIM: the goal is to make the steps from here on out as fast as
% possible. This will become a part of a bigger problem needs to be
% lightning fast!


a

a = struct with fields:
            x0: [1×1 sym]
            y0: [1×1 sym]
        omega0: [1×1 sym]
            E0: [1×1 sym]
    parameters: [1×1 sym]
    conditions: [1×1 sym]

a.parameters

ans = 

x

a.conditions

ans = 

Pm_val = PmSol;
% substitute all Pm values at once in a.conditions to get equation for
% parameter
subs_new = subs(a.conditions, Pm0, Pm_val);  
subs_new(1,1)

ans = 

% Now all we need to do is solve each equation in subs_new (each of which
% looks like the above) get the parameter values, combine that with Pm_val
% and substitue in the expressions for x0 to E0 to get equilibrium points
% for all values of Pm

% THE ISSUE: HOW TO VECTORIZE THIS? (the following for loop)

for j = 1:length(PmSol)
    % solve the each equation in subs_new for the parameter
    parameter_val = solve(subs_new(j));
    % make a vector of same dimensions as parameter vals with value
    % Pm_val(j)
    Pm_new = parameter_val * 0 + Pm_val(j);
    
    % substitute parameter_val and Pm_new to get Equilibrium points
    T1 = {[parameter_val], [Pm_new]};
    T2 = {a.parameters, Pm0};
    xFP = subs(a.x0, T2, T1);
    yFP = subs(a.y0, T2, T1);
    omegaFP = subs(a.omega0, T2, T1);
    EFP = subs(a.E0, T2, T1);
    
    [rows, cols] = size([xFP, omegaFP, EFP]);
    % if no solutions add a zero matrix for the equilibrium point locations
    if rows == 0
        Fixed_points_locations(1:2,3*(j-1)+1:3*(j-1)+3) = ones(2, 3)*0;
        continue;
    end
    Fixed_points_locations(1:2,3*(j-1)+1:3*(j-1)+3) = [xFP, omegaFP, EFP];

end

Answer 3

[Ines Mir]