Solve equation that has a complex subexpression

I want to solve the following equation for omega:

where

So I tried this:

syms s omega G(s)
G(s) = 10/(s*(1+s)*(1+0.2*s));
% Try to find omega that satisfies the equation:
solve(angle(subs(G(s),s,omega*j))-deg2rad(-135),omega,'Real',true)

Result:

Error using mupadengine/feval_internal (line 172)

No complex subexpressions allowed in real mode.

Error in solve (line 293)

sol = eng.feval_internal('solve', eqns, vars, solveOptions);

Although there is an imaginary number in the expression, the decision variable is real and the expression evaluates to a real number (due to angle) so I don't see why it should have a problem solving this.

Obviously, I can think of other ways to solve the problem, but it would be nice to just use angle on the whole transfer function.

% Get solution a different way:
omega_sol = solve(-pi/2-atan2(omega,1)-atan2(omega,5)-deg2rad(-135),omega)
% Confirm solution:
subs(angle(subs(G(s),s,omega*j))-deg2rad(-135),omega,omega_sol)

omega_sol =

0.7417

ans =

-1.8367e-40

In summary, is there any way to solve the original expression for omega directly:

angle(subs(G(s),s,omega*j)) == deg2rad(-135)

Answer 1

[Star Strider]

Solving for the tangent of the phase angle, rather than using the arctangent of the transfer function, appears to produce the correct result:

syms s omega G(s)
assume(omega > 0)
G(s) = 10/(s*(1+s)*(1+0.2*s));
G = subs(G, s, 1j*omega)
OMG = solve(imag(G)/real(G) == tan(deg2rad(-135)), omega)
vpaOMG = vpa(OMG)

producing:

vpaOMG =
0.74165738677394138558374873231655

.

[Bill Tubbs]

Thanks. This solves it. B.t.w. to avoid redefining G, this also works: OMG = solve(imag(G(j*omega))/real(G(j*omega)) == tan(deg2rad(-135)), omega).

[Star Strider]

As always, my pleasure!

I thought about using ‘1j*omega’ as a function argument, however went with subs because that was in your original code, and there was some reason you specifically used it.

Answer 2

[Bill Tubbs]

I just discovered that you can also solve this numerically with vpasolve:

syms s omega G(s)
assume(omega > 0)
G(s) = 10/(s*(1+s)*(1+0.2*s));
% Try to find omega that satisfies the equation:
vpasolve(angle(G(omega*j)) == deg2rad(-135),omega)

ans =

0.74165738677394138558374873231655

This is a more robust solution as it can handle more complex functions such as this:

G(s) = exp(-4*s)/(1+s);
vpasolve(angle(G(omega*j)) == deg2rad(-135),omega)

ans =

0.47764713626095932403299027979129

[Bill Tubbs]

Actually, there is a problem with using angle. It wraps the angle so that it is always in the range -180˚ to 180˚:

It got the right answer in the above case but it could easily have found any number of wrong solutions. So I don't think using angle is a good solution here.

[Star Strider]

The unwrap function may be helpful here.

[Bill Tubbs]

Thanks, yes I looked at that but I think it requires a vector of consecutive angles as input. Here, we are passing the function to vpasolve so I don't think unwrap could be used.

[Star Strider]

I’m not certain what you’re plotting.

Experiment with something like this:

ad = -180:20:180;
ad360 = mod(ad+360,360);

ar = -pi:0.31:pi;
ar2pi = mod(ar+2*pi,2*pi);

.