Solve equation that has a complex subexpression

10 views (last 30 days)
Bill Tubbs on 30 Jul 2020
Commented: Star Strider on 2 Dec 2020
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:
Result:
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:
% Confirm solution:
omega_sol =
0.7417
ans =
-1.8367e-40
In summary, is there any way to solve the original expression for omega directly:

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
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:
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);
ans =
0.47764713626095932403299027979129
Star Strider on 2 Dec 2020
I’m not certain what you’re plotting.
Experiment with something like this:
ar = -pi:0.31:pi;
ar2pi = mod(ar+2*pi,2*pi);
.

R2019b

Community Treasure Hunt

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

Start Hunting!