非线性常微分方程画分岔图问题

dairyqzy

应用四阶Runge-Kutta方法直接迭代求解非线性微分方程组,迭代步长为0.001.
方程的M文件是：
function dx=shiyan(t,x,flag,d)
a1=4.730040745;a2=7.853204624;a3=0.005;a4=0.2;a5=10;a6=50;a7=0.05;u=7.6;
n=20;p=20;
b12=0.835504011618438;b21=-b12;g11=-12.302618623949989;
g22=-46.050120142194714;h11=-6.151309311974995;h12=-8.380326623605050;
h21=-6.709318600368174;h22=-23.025060071097357;e=-7.544822611986612;
format long
dx=zeros(4,1);
dx(1)=x(3);
dx(2)=x(4);
dx(3)=-(a1^4+(u^2-1/2*a6)*g11)*x(1)-a6*e*x(2)-a3*a1^4*x(3)-2*a4^0.5*u*b12*x(4)+d*n^2*sin(n*t)*(sinh(a1)-sin(a1)-((cosh(a1)-cos(a1))/(sinh(a1)-sin(a1)))*(cos(a1)+cosh(a1)-2))/a1-a7*g11*x(1)*(2*u*sin(p*t)+a4^0.5*p*cos(p*t))+a4^0.5*p*a7*(h11*x(1)+h12*x(2))*cos(p*t)-2*a4^0.5*a7*b12*x(4)*sin(p*t)-a5*(((x(1)^2+2*a3*x(1)*x(3))*g11+(x(2)^2+2*a3*x(2)*x(4))*g22)*x(1)*g11);
dx(4)=-a6*e*x(1)-(a2^4+(u^2-1/2*a6)*g22)*x(2)-2*a4^0.5*u*b21*x(3)-a3*a2^4*x(4)+d*n^2*sin(n*t)*(sinh(a2)-sin(a2)-((cosh(a2)-cos(a2))/(sinh(a2)-sin(a2)))*(cos(a2)+cosh(a2)-2))/a2-a7*g22*x(2)*(2*u*sin(p*t)+a4^0.5*p*cos(p*t))+a4^0.5*p*a7*(h21*x(1)+h22*x(2))*cos(p*t)-2*a4^0.5*a7*b21*x(3)*sin(p*t)-a5*(((x(1)^2+2*a3*x(1)*x(3))*g11+(x(2)^2+2*a3*x(2)*x(4))*g22)*x(2)*g22);
分岔参数为 d
在分岔图中使用的触发条件是1.588146262064604*y(:,3)+6.968647880967183*10^(-11)*y(:,4)从正向趋于0的时候，在图中记录下参数n激励时1.588146262064604*y(:,1)+6.968647880967183*10^(-11)*y(:,2)的值（即以d为横坐标，以1.588146262064604*y(:,1)+6.968647880967183*10^(-11)*y(:,2)的值为纵坐标画图）；
解方程及其画分岔图的程序为：
d=0:0.0001:0.11;
for j=1:length(d)
[t,y]=ode45('shiyan',[0:0.001:100],[-0.001,-0.001,0,0],[],d(j));
    data=y(:,1);
    m=length(data);
    for k=1:m-1
        if 1.588146262064604*y(k,3)+6.968647880967183*10^(-11)*y(k,4)>=0 && 1.588146262064604*y(k+1,3)+6.968647880967183*10^(-11)*y(k+1,4)<0;
        plot(d(j),1.588146262064604*y(k,1)+6.968647880967183*10^(-11)*y(k,2));
        hold on
        else
        end
    end
end
这样做出来的完全对不上，是不是我的分岔程序有问题，还是其他的地方出了问题，希望各位如果有时间帮忙看一下，到底哪里出了问题，我需要改进哪些部分。以后如果画图需要注意什么，因为这些真的只能自己琢磨，没有师兄的指点，有时候真的很头疼，先提前谢谢大家了~~

为了节约计算时间，自己画图时d的步长选择为d=0:0.001:0.11;


我的模型是简谐激励振动条件下直管道的模型，模型为两个2阶非线性常微分方程组，模型是文献上的，应该没错误的，参数是根据文献上得来的，我主要是想通过文献验证来验证编写的程序。模型的求解是通过采用数值模拟的方法对系统运动形态进行仿真研究，主要考察流体平均流速、流速脉动振幅和基础简谐运动激励振幅对系统动态响应的影响。应用四阶Runge-Kutta方法直接迭代求解方程。迭代步长为0.001，而触发条件是：管道中点的速度从正向趋于0的时候，记下管道中点此时的位移。即式子：速度：1.588146262064604*y(:,3)+6.968647880967183*10^(-11)*y(:,4)从正向趋于0；位移：记下1.588146262064604*y(:,1)+6.968647880967183*10^(-11)*y(:,2)的值并描在图上（即纵坐标的值）。


接触非线性时间不是很长，希望各位大大帮帮忙看看，帮我指出错误。。。    在这里先谢谢各位了。

dairyqzy

没人能帮忙看看嘛？

驻印尼大使