如何绘制变参数根轨迹(针对复杂情况，无法分离出开环增益k*时)

最近阅读论文遇到了一个很棘手的问题，如何绘制变参数根轨迹方程，下面通过合工大博士论文来进行分析:
问题


给出闭环特征方程，只有J在0.5-8之间变化，其余的参数都是常量，论文中已经给出，下面就是要绘制当J变化时系统的根轨迹曲线，以此来表征系统在参数变化时的稳定性。
分析
传统方式：
通过化简，整成1+J*G(S)=0，这是自控说的。绘制等效开环的根轨迹，限制一下开环增益J，就能够得到闭环的根轨迹。**但是！**你会发现开环增益并不是单纯的J，算式太复杂，无法正常分离！
本质出发：
咱们的任务就是求解当J在一个范围内变化时闭环的根轨迹。说白了，就是求解特征方程在某一参数下的根！！
有了这个思路，利用solve求解数学方程函数，提取出方程的解，利用scatter(real(x),imag(x))绘制特征根的散点图，这就是在J参数聂动时闭环的根轨迹！
代码块

%合工大石荣亮博士论文
%J变化时参数根轨迹
clear all;
clc;
syms s;  %syms用于定义未知的变量
U=220;E=220;
theta=0.02;
X=0.1;W0=314;D=2;
m=3.1e-5;n=7.1e-5;
kpe=U*theta/X;
kpd=E*U/X;
kqe=(2*E-U)/X;
kqd=U*E/X;
for J=0.5:0.1:8
a=(2+n*kqe)*(D*W0+1/m)/(J*W0);
b=kpd/(J*W0)+(n*kqe+1)*(D*W0+1/m)^2/(J*W0)^2;
c=(kpd+n*kqe*kpd-n*kpe*kqd)*(D*W0+1/m)/(J*W0)^2;
x=solve(s^3+a*s^2+b*s+c,s);   %求解出的s是一个关于J的函数
y=vpa(x,2);                   %对求解的结果保留两位小数
e=y(1);f=y(2);g=y(3);         %取出计算的结果
scatter(real(e),imag(e),'rx');%画出散点图
hold on;
scatter(real(f),imag(f),'b');
scatter(real(g),imag(g),'k');
hold on;
grid on;
end
title('J变化时参数根轨迹');
xlabel('实部:Re');
ylabel('虚部:Im');
axis([-60,0,-15,15]);

结果

可以看出，该方法能够实现参数聂动时无法分离k*的根轨迹绘制！
最后
MATLAB菜鸟一枚，代码还有较大的提升空间，欢迎大家与我交流讨论！