14 views (last 30 days)

Show older comments

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)

Star Strider
on 30 Jul 2020

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

.

Star Strider
on 30 Jul 2020

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.

Bill Tubbs
on 2 Dec 2020

Edited: Bill Tubbs
on 2 Dec 2020

Star Strider
on 2 Dec 2020

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);

.

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

Start Hunting!