求助：如何画二阶微分方程的分岔图？

aliaofly

各位学长，小弟刚学MATLAB，求解一非线性微分方程后需要画分岔图，可小弟不会啊！
解出来的是位移x1关于时间的函数，但怎么画出位移关于阻尼a、频率c的分岔图呢？
微分方程为：dx1=x2
                     dx2=-2*a*x2-[1+k1*cos(98*c*t)+k2*cos(4*c*t)]*f+p0+p1*cos(c*t)
其中f为一非线性函数
下面是程序：
function [x1,x2]=vdp0(x1,x2)
x1=zeros(1,201);
x2=zeros(1,201);
dx1=x2;
f=zeros(1,201);
a=0.08;
c=0.5;
k1=0.1;
k2=0.1;
p0=0.1;
p1=0.2;
x1(1)=0;
x2(1)=1;
for i=2:1:201   %   下面用四阶龙格库塔法解二阶微分方程
    k1=2.25e-3*(i-2)
    t=0.01*(i-1)
if  x1(i-1)>1
    f(i-1)=x1(i-1)-1;
elseif  x1(i-1)<-1
            f(i-1)=x1(i-1)+1;
else
        f(i-1)=0;
        end
K1=-2*a*x2(i-1)-[1+k1*cos(98*c*t)+k2*cos(4*c*t)]*f(i-1)+p0+p1*cos(c*t);
K2=-2*a*(x2(i-1)+0.005*K1)-(1+k1*cos(98*c*(t+0.005))+k2*cos(4*c*(t+0.005)))*f(i-1)+p0+p1*cos(c*(t+0.005));
K3=-2*a*(x2(i-1)+0.005*K2)-(1+k1*cos(98*c*(t+0.005))+k2*cos(4*c*(t+0.005)))*f(i-1)+p0+p1*cos(c*(t+0.005));
K4=-2*a*(x2(i-1)+0.01*K3)-(1+k1*cos(98*c*(t+0.01))+k2*cos(4*c*(t+0.01)))*f(i-1)+p0+p1*cos(c*(t+0.01));
x2(i)=x2(i-1)+0.01*(K1+2*K2+2*K3+K4)/6;
k1=x2(i-1);
k2=x2(i-1)+0.005*k1;
k3=x2(i-1)+0.005*k2;
k4=x2(i-1)+0.01*k3;
x1(i)=x1(i-1)+0.01*(k1+2*k2+2*k3+k4);
end

咕噜噜

你的程序用的是最大值法，最大值法论坛里面有相关的程序，你找找看
分叉图的画法论坛里面也有很多，你参考一下

aliaofly

这几天一直在看分岔版块的，我把程序修改了一下 ，但一直有错误，折腾了很长一段时间了，老板像催魂一样像我要结果，郁闷啊，还恳请小咕多多指教！
function xd=fun(t,x,flag,w)
a=0.08;
k1=0.1;
k2=0.1;
p0=0.1;
p1=0.2;
  if  x(1)>1
    f=x(1)-1;
elseif  x(1)<-1
    f=x(1)+1;
else
        f=0;
        end
xd=[x(2);
    -2*a*x(2)-(1+k1*cos(98*w*t)+k2*cos(4*w*t))*f+p0+p1*cos(w*t)];

w=0:0.01:1;
for i=1:length(w)
T=2*pi/w(i);
[t,x]=ode45(@fun,[0:T/200:200*T],[0,0],[],w(i));
plot(w(i),x(36000:200:end,1),'k.'),hold on;
end
其中方程为：x'' +2*a*x' +(1+k1*cos(98*w*t)+k2*cos(4*w*t))*f =p0+p1*cos(w*t)

yanzi

[t,x]=ode45(@fun,[0:T/200:200*T],[0,0],[],w(i));
plot(w(i),x(36000:200:end,1),'k.'),hold on
一共只有20000个点 ，怎么消除瞬态有36000个点？

yibanlixue

200*200不是40000个点么

大漠苍狼

用专业软件,winpp,一下子就可以画了,不用编程的

无水1324

winpp可以求任意系统吗？

octopussheng

winpp？？第一次看到，是否有介绍？

shanyun

WinPP软件是近年出现的一种研究非线性系统动力学的重要软件,它功能强大,简单实用,是设计用来求形式为常微分方程(组) 、时滞微分方程(组) 、微分代数方程(组) 、Volterra积分微分方程(组) 、离散动力系统、边值问题等动力系统的解的程序. 通过创建一个扩展名为". ode"的文件,这个文件包含了WinPP对动力系统进行数值模拟所必需的方程、参数、函数等元素,并且可以在文件中定义画图参数. 其界面是基于Win95 /NT的. 除了变量名称和维数不能改变、不退出不能使用新的文件以外,模拟的各个方面都是可以改变的. 它具有画动力系统的相图、时间历程、Poincare截面、寻找平衡点、利用单参数区域积分和双参数区域积分来发现初值在某个区域内的动力学现象、提取数据进行数值解和解析解的对比、保存图像以加入Latex、Word文档等强大功能, 可以进行稳定性、混沌、分岔等动力学分析、验证[1].

[1] 郭基凤,裴利军. WinPP软件在一类捕食- 被捕食系统动力学研究中的应用[ J ]. 信阳师范学院学报:自然科学版, 2007, 20 (4) : 482-485.

无水1324

“除了变量名称和维数不能改变”这个是什么意思，它只能做2维的？

[ 本帖最后由 无水1324 于 2009-8-10 09:48 编辑 ]

shanyun

winpp我也没用过，前面介绍是从一片论文里摘取的...

尝试过下载，没成功...   似乎所有下载连接都指向http://www.math.pitt.edu/~bard/bardware/   但我打不开

无水1324

试了一下，可以打开的

shanyun

能发给我么？ shanyunh@QQ.com   谢谢！

jasonpq

winpp出图的确方便。。。

bohua1208

回复 4 # yanzi 的帖子

消除瞬态是什么意思，麻烦解释一下，最近看了一篇文章，说积分1400个周期，消除前300个周期瞬态解是什么意思？