hsfy919 发表于 2010-2-21 11:48

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%这是方程表达式,我把b定义为全局变量
functiondy=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=;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%以下为主函数程序
clc
clear
global b f
Y1=[];
f=0.4:0.001:1.2;    %可能步长较小,计算时间长点,可以改变步长算出大概轮廓也行
k=0;
y0=;
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;
    =ode23tb('ODE_fun',tspan,y0);             %%%%%%%%%%%%%%因为微分方程刚性,故用考虑用23tb求解
    y0=y(end,:);
    j=1;
    for p=60:200
      tspan=p*period:step:(p+1)*period;
      =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 发表于 2010-2-21 15:01

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

hsfy919 发表于 2010-2-21 23:27

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

xinwilliam 发表于 2010-2-23 20:06

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

hsfy919 发表于 2010-2-23 22:47

回复 地板 xinwilliam 的帖子

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

无水1324 发表于 2010-3-10 16:25

你这样子写程序是不是计算的时间很长啊?
页: [1]
查看完整版本: 大家看看我的分岔图吧,提点意见