6 views (last 30 days)

Show older comments

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?

John D'Errico
on 8 Sep 2021

Edited: John D'Errico
on 8 Sep 2021

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.

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!