Hopf分岔点搜索、同宿分岔，高手帮忙看看，是不是程序哪里有问题？
本来是想做一个搜索Hopf分岔点的程序，结果在分岔参数Q1=11.14时，状态方程的解出现NAN的情况，和实例对比，发现在Q1=11.2的时候发生同宿分岔。百思不得其解。对照了实例，没有发现有输入错误的地方。特求助各位达人！望给予指导！用matcont做了仿真，发现Q1=11.14没有什么异常情况，和实例吻合的很好。
实例中，(X,Q1)= 为Hopf点
(X,Q1)= 为Hopf点
(X,Q1)= 为SNB点
///////////////////////////////////////////////////////非线性微分方程组
function F=myfun(t,x,Q1)
F=[x(2);
16.6667*sin(x(3)-x(1)+0.0873)*x(4)-0.1667*x(2)+1.8807;
496.8718*x(4)^2-166.6667*cos(x(3)-x(1)-0.0873)*x(4)-666.6667*cos(x(3)-0.2094)*x(4)-93.3333*x(4)+33.3333*Q1+43.3333;
-78.7638*x(4)^2+26.2172*cos(x(3)-x(1)-0.0124)*x(4)+104.8689*cos(x(3)-0.1346)*x(4)+14.5229*x(4)-5.2288*Q1-7.0327];
///////////////////////////////////////////////////一下是主程序
%搜索Hopf分岔点
tic
A0=eye(4);
A=j*100*pi*A0;
h=1e-2;
t0=0; %开始时间
tf=15;%终了时间
x0=; %初值
Q1=10:0.01:12;%分岔参数取值范围
for k=1:length(Q1) % length(Q1)为Q1的点数
kk=Q1(k); %每一个Q1值赋给kk
for t=t0:h:tf %四阶龙格库塔法解方程组
k1=myfun(t,x0,kk);
k2=myfun(t+h/2,x0+h/2*k1,kk);
k3=myfun(t+h/2,x0+h/2*k2,kk);
k4=myfun(t+h,x0+h*k3,kk);
x=x0+h/6*(k1+2*k2+2*k3+k4);
x0=x;
end
x1(k)=x(1);
x2(k)=x(2);%把每一个Q1值下的解保存
x3(k)=x(3);
x4(k)=x(4);
J=[0,1,0,0;
-16.6667*cos(x(3)-x(1)+0.0873)*x(4),-0.1667,16.6667*cos(x(3)-x(1)+0.0873)*x(4),16.6667*sin(x(3)-x(1)+0.0873);
-166.6667*sin(x(3)-x(1)-0.0873)*x(4),0,166.6667*sin(x(3)-x(1)-0.0873)*x(4)+666.6667*sin(x(3)-0.2094)*x(4),496.8718*2*x(4)-166.6667*cos(x(3)-x(1)-0.0873)-666.6667*cos(x(3)-0.2094)-93.3333;
26.2172*sin(x(3)-x(1)-0.0124)*x(4),0,-26.2172*sin(x(3)-x(1)-0.0124)*x(4)-104.8689*sin(x(3)-0.1346)*x(4),-78.7638*2*x(4)+26.2172*cos(x(3)-x(1)-0.0124)+104.8689*cos(x(3)-0.1346)+14.5229];
%构造In,m
if det(A-J)<=1e-3&& myfun(t,x,Q1)<=1e-3 %判断是否为Hopf分岔点
XHopf=x;
disp('Q1=')
disp(Q1(k))
disp('Hopf Bifurcation XHopf=')
disp(XHopf)
end
end
plot(Q1,x1,'-b',Q1,x4,'-r',Q1,x3,'-y');
toc
////////////////////////////////////////////////////////////////////
Hopf分岔点判据：(见《基于分岔理论的电力系统电压稳定性分析》一书P12)
F(X,Q1)=0 && det(j*2*pi*f*A)-J(X,Q1))=0
J为雅克比矩阵
上式中，A=eye(4);
//////////////////////////////////////////////////////////////////
附图：
红色为x(4) 蓝色为x(1) 黄色为x(3) 这是用matcont仿真出来的图形：
你是用A MODEL OF VOLTAGE COLLAPSE IN ELECTRIC POWER SYSTEMS它做的吧。
页:
[1]