大家看看我的分岔图吧，提点意见

hsfy919

我的程序是求解二阶常微分方程组，我把它化成了状态方程，用变步长龙格库塔法求解，按照论坛里的方法做出的分插图怪怪的，恳请大家指正
我想画y1随b的分岔图，b的变化范围是0.4到1.2

%%%%%%%%%%%%%%%%%%%%%%%%%这是方程表达式，我把b定义为全局变量
function  dy=ODE_fun(t,y)
sigma=9;delta=0.05;a=0.01;n=10;k=10;pai=pi;

global b

e2=y(1)^2+y(3)^2;
e=sqrt(e2);
sss=y(1)/e;
ccc=-y(3)/e;
de=(y(1)*y(2)+y(3)*y(4))/e;
df=(y(1)*y(4)-y(3)*y(2))/e2;
q1=1-2*df;
q2=2*de;
E1=2*e2/((1-e2)*(2+e2));
E2=(pai/2-8/(pai*(2+e2)))./((1-e2)^1.5);
E3=pai*e/(((1-e2)^1.5)*(2+e2));
E4=2*e/((1-e2)*(2+e2));
fr1=sigma*b*(q1*E1+q2*E2);
ft1=sigma*b*(q1*E3+q2*E4);
c=-fr1*sss+ft1*ccc;
d=fr1*ccc+ft1*sss;
e2=(y(9)-delta)^2+y(11)^2;
e=sqrt(e2);
sss=(y(9)-delta)/e;
ccc=-y(11)/e;
de=((y(9)-delta)*y(10)+y(11)*y(12))./e;
df=((y(9)-delta)*y(12)-y(11)*y(10))./e2;
q1=1-2*df;
q2=2*de;
E1=2*e2/((1-e2)*(2+e2));
E2=(pai/2-8/(pai*(2+e2)))/((1-e2)^1.5);
E3=pai*e/(((1-e2)^1.5)*(2+e2));
E4=2*e/((1-e2)*(2+e2));
fr1=sigma*b*(q1*E1+q2*E2);
ft1=sigma*b*(q1*E3+q2*E4);
f=-fr1*sss+ft1*ccc;
g=fr1*ccc+ft1*sss;
dy(1)=y(2);
dy(3)=y(4);
dy(5)=y(6);
dy(7)=y(8);
dy(9)=y(10);
dy(11)=y(12);
dy(2)=1/(b.^2)*c-1/(b.^2)*k*(y(1)-y(5))+1/(b.^2);
dy(4)=1/(b^2)*d-1/(b^2)*k*(y(3)-y(7));
dy(6)=-k*(y(5)-y(1))/(n*(b.^2))-k*(y(5)-(y(9)-delta))/(n*(b^2))+1/(b^2)+a*cos(t);
%dy(6)=-k*(y(5)-y(1))/(n*(b.^2))-k*(y(5)-y(9))/(n*(b^2))+1/(b^2)+a*cos(t);
dy(8)=-k*(y(7)-y(3))/(n*(b^2))-k*(y(7)-y(11))/(n*(b^2))+a*sin(t);
dy(10)=f/(b^2)-k*(y(9)-delta-y(5))/(b^2)+1/(b^2);
%dy(10)=f/(b^2)-k*(y(9)-y(5))/(b^2)+1/(b^2);
dy(12)=g/(b^2)-k*(y(11)-y(7))/(b^2);
dy=[dy(1);dy(2);dy(3);dy(4);dy(5);dy(6);dy(7);dy(8);dy(9);dy(10);dy(11);dy(12)];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%以下为主函数程序
clc
clear
global b f
Y1=[];
f=0.4:0.001:1.2;    %可能步长较小，计算时间长点，可以改变步长算出大概轮廓也行
k=0;
y0=[0.1;0;0.1;0;0.2;0;0;0;0;0;0;0];
for i=1:length(f)
    disp(f(i))
    period=2*pi/f(i);
    b=f(i);
    step=period/100;
    k=k+1;
    tspan=0:step:60*period;  
    [t,y]=ode23tb('ODE_fun',tspan,y0);             %%%%%%%%%%%%%%因为微分方程刚性，故用考虑用23tb求解
    y0=y(end,:);  
    j=1;
    for p=60:200
        tspan=p*period:step:(p+1)*period;
        [t,y]=ode23tb('ODE_fun',tspan,y0);
        Y1(k,j)=y(end,1);   
        j=j+1;
        y0=y(end,:);
    end
end
    bifdata=Y1(:,end-50:end);
    plot(f,bifdata,'r.','markersize',1)
    hold on
%%%%%%%%%%%%%%%%%%%%%%%%%%%以下是我按上面的程序算出来的分岔图，怪怪的，估计是哪个地方不对，恳请指正

hsfy919

希望大家都发表一下自己的意见，讨论一下吧，我初学，谢谢大家了

hsfy919

盼望octo, 咕噜噜, 无水，给予指导啊

xinwilliam

画图的算法不对，你是用离散动力系统的方法来画连续动力系统的分岔图。
对于连续动力系统分岔图的画法，在本版中已经有好多解释了，还是搜索一下就行了。

hsfy919

您好，我是初学，想问一下，数值解法不就是离散化吗？用符号求解不出来的。怎么画连续的图啊，有没有例子，或给个链接，我在论坛上找了，没找到，谢谢

无水1324

你这样子写程序是不是计算的时间很长啊？