Solve equation that has a complex subexpression

11 次查看（过去 30 天）
Bill Tubbs
Bill Tubbs 2020-7-30
评论： Star Strider 2020-12-2
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
Star Strider 2020-7-30
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
.
  2 个评论
Bill Tubbs
Bill Tubbs 2020-7-30
编辑：Bill Tubbs 2020-7-30
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
Star Strider 2020-7-30
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.

请先登录，再进行评论。

更多回答（1 个）

Bill Tubbs
Bill Tubbs 2020-12-2
编辑：Bill Tubbs 2020-12-2
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
  4 个评论
显示 2更早的评论
Bill Tubbs
Bill Tubbs 2020-12-2
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
Star Strider 2020-12-2
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);
.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numeric Solvers 的更多信息

标签

产品


版本

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by