Error in code for MATlab

clear all
syms s w 

G = 1/((s)*(s+1)*(s+2)); %transfer function
G_w = subs(G,s,j*w); 
W= [-1000:0.1:1000]; %[min_range:step size:max_range]
nyq = eval(subs(G_w,w,W)); 

x = real(nyq)
y = imag(nyq)

plot(x,y)

I can't seem to run this code and it keeps displaying error in line 100++ where I've only less than 20 lines.

This are the errors shown:

Error using symengine (line 59)
Division by zero.

Error in sym/subs>mupadsubs (line 139)
G = mupadmex('symobj::fullsubs',F.s,X2,Y2);

Error in sym/subs (line 124)
    G = mupadsubs(F,X,Y);

Error in nyquist2 (line 8)
nyq = eval(subs(G_w,w,W)); %replace W with w in equation G_w

Does any expert know the cause of the problem?

[Image Analyst]

I've fixed the formatting for you this time, but for next time, please read this link to learn how to format.

[Niam]

Thank you for the help in formatting, was new and didn't know.

My apologies.

Answer 1

[Star Strider]

If you have the Control System Toolbox, this works:

s = tf('s');
G = 1/((s)*(s+1)*(s+2)); %transfer function

figure(1)
pzmap(G)

figure(2)
rlocus(G)

figure(3)
nyquist(G)

EDIT — I didn’t notice the ‘nyquist’ tag at first. Added rlocus and nyquist calls.

[Niam]

Hi there,

As i need the wmin,wmax range to be included so i can test on the speed of plotting the nyquist, is there any possible way to add that in ?

Thank you for your reply !

[Star Strider]

My pleasure!

From this excerpt from the documentation for the nyquist function:

• nyquist(sys,w) explicitly specifies the frequency range or frequency points to be used for the plot. To focus on a particular frequency interval, set w = {wmin,wmax}. To use particular frequency points, set w to the vector of desired frequencies.

So, the short answer is: yes.

See the documentation for the nyquist function for a full description and all the options.

[Niam]

Thanks once again, i think i got it.

But i'm still figuring out my initial coding, do you by chance understand whats wrong ? Sorry for my lack of knowledge in it !

[Star Strider]

Part of the problem with your initial coding is the subtle fact that the Laplace variable ‘s’ is ‘sigma ± j*w’. (This is necessary because of the symmetry of the complex plane.) You then have to supply a range for both ‘sigma’ (the real part) and ‘w’ (the complex frequency).

Your code, slightly revised:

syms s w sigma
G = symfun(1/((s)*(s+1)*(s+2)), s); %transfer function
H = expand(G(sigma+1i*w) * G(sigma-1i*w));
H = simplify(H, 'Steps',10);
Hfun = matlabFunction(H)

producing:

Hfun = (sigma,w) 1.0./(sigma.^2.*w.^2.*1.8e1+sigma.^3.*w.^2.*1.2e1+sigma.^2.*w.^4.*3.0+sigma.^4.*w.^2.*3.0+sigma.*w.^2.*1.2e1+sigma.*w.^4.*6.0+sigma.^2.*4.0+sigma.^3.*1.2e1+sigma.^4.*1.3e1+sigma.^5.*6.0+sigma.^6+w.^2.*4.0+w.^4.*5.0+w.^6);

that you can use outside of the Symbolic Math Toolbox. I leave the rest to you. If you use meshgrid with vectors for ‘sigma’ and ‘w’ and then plot it with surf or mesh. you should be able to see the entire complex surface with the poles. (There are zeros only at ±Inf in your transfer function.) I have not done the plot.