# Solve equation that has a complex subexpression

조회 수: 6(최근 30일)
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:
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 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표시숨기기 이전 댓글 수: 1
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 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:
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
##### 댓글 수: 4표시숨기기 이전 댓글 수: 3
Star Strider 2020년 12월 2일
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!