非线性常微分方程画分岔图问题
<h1 style="color:# 麦片财富积分
新手, 积分 5, 距离下一级还需 45 积分
dairyqzy 发表于
应该是可以的啊，我看过有程序这么做过：动力学方程为
M 函数文件：dbfun.m
function ydot=dbfun(t,y,fl ...
能不能帮忙看看啊
论坛优秀回答者
<h1 style="color:#6 麦片财富积分
关注者： 11
能不能把你的结果（图）贴出来看看
<h1 style="color:# 麦片财富积分
fengshenone 发表于
能不能把你的结果（图）贴出来看看
d=0:0.001:0.11;
for j=1:length(d)&&上面有个小错误，应该将n改成d;
[t,y]=ode45('shiyan',[0:0.001:10],[-0.001,-0.001,0,0],[],d(j));
& & data=y(:,1);
& & m=length(data);
& & for k=1:m-1
& && &&&if 1.604*y(k,3)+6.183*10^(-11)*y(k,4)&=0 && 1.604*y(k+1,3)+6.183*10^(-11)*y(k+1,4)&0;
& && &&&plot(d(j),1.604*y(k,1)+6.183*10^(-11)*y(k,2));
& && &&&hold on
& && &&&else
& && &&&end
论坛优秀回答者
<h1 style="color:#6 麦片财富积分
关注者： 11
你是想画出和你给的例子一模一样的图吗？？不考虑内容的话，你的程序本身应该没有问题，我也算了一下，和你画的图一样。 你看你给的例子在d=0的时候纵坐标是0.28，而你画出来的纵坐标是0，说明很有可能是你本身程序的参数或者你的数学模型写错了，如果你参考某篇文献的，建议你还是给出这个数学模型的具体说明，不然的话也没什么人能帮你吧。
<h1 style="color:# 麦片财富积分
fengshenone 发表于
你是想画出和你给的例子一模一样的图吗？？不考虑内容的话，你的程序本身应该没有问题，我也算了一下，和你 ...
对& & 我就是要画出例子一样的图，我的模型是简谐激励振动条件下直管道的模型，模型为两个2阶非线性常微分方程组，模型是文献上的，应该没错误的，参数是根据文献上得来的，我主要是想通过文献验证来验证编写的程序。模型的求解是通过采用数值模拟的方法对系统运动形态进行仿真研究，主要考察流体平均流速、流速脉动振幅和基础简谐运动激励振幅对系统动态响应的影响。应用四阶Runge-Kutta方法直接迭代求解方程。迭代步长为0.001，而触发条件是：管道中点的速度从正向趋于0的时候，记下管道中点此时的位移。即式子：速度：1.604*y(:,3)+6.183*10^(-11)*y(:,4)从正向趋于0；位移：记下1.604*y(:,1)+6.183*10^(-11)*y(:,2)的值并描在图上（即纵坐标的值）。即例子上面的图。
<h1 style="color:# 麦片财富积分
没有人了么？帮忙分析一下，不胜感激。。
<h1 style="color:# 麦片财富积分
y(500:100:1400,2)中的2是什么意思？
<h1 style="color:# 麦片财富积分
这个要运行好久才出结果？
<h1 style="color:# 麦片财富积分
d=0:0.001:0.11;
for j=1:length(d)&&上面有个小错误，应该将n改成d;
[t,y]=ode45('shiyan',[0:0.001:10] ...
这个代码运行得到的就是这个结果！
站长推荐 /3
Powered byMatlab——混沌分岔图_逝水流年_新浪博客
Matlab——混沌分岔图
经过近期的研究发现，目前对于系统单参数分岔图的计算共有以下的几种方法：
1）最大值法
即对系统微分方程（组）进行求解，对求解的结果用getmax函数进行取点，并绘图。
2）Poincare截面法
对系统参数的每一次取值，绘制其Poincare截面，进而得到其分岔图。
这种方法需要注意的是，自治系统的Poincare截面是选取一超平面，平面上点的分布即构成一Poincare截面，非自治系统的Poincare截面则是根据系统激励的频率进行取点并绘图。
本帖将以Lorenz系统为例，对这两种方法进行比较
首先对第二种方法进行阐述。
编程如下（matlab）
Lorenz系统：
function dy = Lorenz(t,y)
% Lorenz系统
% 系统微分方程：
dx/dt = -a(x-y)
dy/dt = x(r-z)-y
dz/dt = xy-bz
% & &r=y(5)
dy=zeros(6,1);
dy(1)=-y(4)*(y(1)-y(2));
dy(2)=y(1)*(y(5)-y(3))-y(2);
dy(3)=y(1)*y(2)-y(6)*y(3);
dy(6)=0;&#8203;
随r的分岔图求解程序：——按照x=y平面取截面
function Lorenz_bifur_r
for r=linspace(1,500,1000);&%
舍弃前面迭带的结果，用后面的结果画图
[T,Y]=ode45('Lorenz',[0,1],[1;1;1;16;r;4]);&
&[T,Y]=ode45('Lorenz',[0,50],Y(length(Y),:));&
Y(:,1)=Y(:,2)-Y(:,1);&%
对计算结果进行判断，如果点满足x=y，则取点
for k=2:length(Y)&
Y(k,1)&0&&
&if Y(f,1)&0
y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));
Z=[Z r+abs(y)*i];
&else&&&&&if
Y(f,1)&0&&y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));&&&
Z=[Z r+abs(y)*i];&&
plot(Z,'.','markersize',1)title('Lorenz映射分岔图')xlabel('r'),ylabel('|y|
where x=y')
getmax法——取最大值法
function [Xmax] = getmax(y)
a=length(y);
for i=(a-1)/2:a
b=(y(i,1)-y(i-2,1))/2;
c=(y(i,1)+y(i-2,1))/2-y(i-1,1);
y(i-2,1)&=y(i-1,1)&y(i-1,1)&=y(i,1)&c==0
Xmax(j)=y(i-1,1);
y(i-2,1)&=y(i-1,1)&y(i-1,1)&=y(i,1)Xmax(j)=y(i-1,1)-b^2/(4*c);
function Lorenz_bifur_r_getmax% 最大值法求解分岔图
t0=[0 100];%积分时间
%bifurcation
for r=linspace(1,500,1000);& %r的变化精度
[t,y]=ode45('Lorenz',t0,[1;1;1;16;r;4]);
[Xmax]=getmax(y(:,1));
plot(r,Xmax,'b','markersize',1)
clear Xmax
&#8203;转载：
博客等级：
博客积分：0
博客访问：3,324
关注人气：0
荣誉徽章：查看: 7729|回复: 18|关注: 0
非线性常微分方程画分岔图问题
<h1 style="color:# 麦片财富积分
新手, 积分 5, 距离下一级还需 45 积分
应用四阶Runge-Kutta方法直接迭代求解非线性微分方程组,迭代步长为0.001.
方程的M文件是：
function dx=shiyan(t,x,flag,d)
a1=4.;a2=7.;a3=0.005;a4=0.2;a5=10;a6=50;a7=0.05;u=7.6;
n=20;p=20;
b12=0.438;b21=-b12;g11=-12.989;
g22=-46.714;h11=-6.995;h12=-8.050;
h21=-6.174;h22=-23.357;e=-7.612;
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.604*y(:,3)+6.183*10^(-11)*y(:,4)从正向趋于0的时候，在图中记录下参数n激励时1.604*y(:,1)+6.183*10^(-11)*y(:,2)的值（即以d为横坐标，以1.604*y(:,1)+6.183*10^(-11)*y(:,2)的值为纵坐标画图）；
解方程及其画分岔图的程序为：
for j=1:length(n)
[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.604*y(k,3)+6.183*10^(-11)*y(k,4)&=0 && 1.604*y(k+1,3)+6.183*10^(-11)*y(k+1,4)&0;
& && &&&plot(d(j),1.604*y(k,1)+6.183*10^(-11)*y(k,2));
& && &&&hold on
& && &&&else
& && &&&end
这样做出来的完全对不上，是不是我的分岔程序有问题，还是其他的地方出了问题，希望各位如果有时间帮忙看一下，到底哪里出了问题，我需要改进哪些部分。以后如果画图需要注意什么，因为这些真的只能自己琢磨，没有师兄的指点，有时候真的很头疼，先提前谢谢大家了~~
<h1 style="color:# 麦片财富积分
有能指点一下的吗？
<h1 style="color:# 麦片财富积分
:'(:'(:'(:'(:'(:'(:'(:'(:'(:'(:'(
<h1 style="color:# 麦片财富积分
╮(╯▽╰)╭& & 有人能可以帮忙下嘛？有偿~~~&&求大神
MATLAB 数学、统计与优化版块优秀回答者
关注者： 17
dairyqzy 发表于
╮(╯▽╰)╭& & 有人能可以帮忙下嘛？有偿~~~&&求大神
[t,y]=ode45('shiyan',[0:0.001:100],[-0.001,-0.001,0,0],[],d(j))
ode45是自适应步长，不可以你自己手动设置，如果你要设置定步长，那你得自己编写一个ode45的程序，百度上随便下载一个就有了。。
<h1 style="color:# 麦片财富积分
本帖最后由 dairyqzy 于
08:42 编辑
zhanglu 发表于
[t,y]=ode45('shiyan',[0:0.001:100],[-0.001,-0.001,0,0],[],d(j))
ode45是自适应步长，不可以你自己手 ...
应该是可以的啊，我看过有程序这么做过：动力学方程为
M 函数文件：dbfun.m
function ydot=dbfun(t,y,flag,d)&&
w=1;c=0.3;
ydot=[y(2);&&
- y(1)^3 - c*y(2)+d*cos(w*t)];
相应程序：
d=20:0.01:40;&&
for j=1:length(d)&&
[t,y]=ode23 ('dbfun',[0:T/100:70*T],[4,4],[],d(j));&&
plot(d(j),y(500:100:1400,2),'linewidth',5)&&
title('分岔图二')&&
xlabel(' 参数&&F0')&&
ylabel(' 位移')&&
这里也是自定义步长的。
是否还有其他原因呢？？？
MATLAB 数学、统计与优化版块优秀回答者
关注者： 17
dairyqzy 发表于
应该是可以的啊，我看过有程序这么做过：动力学方程为
M 函数文件：dbfun.m
function ydot=dbfun(t,y,fla ...
其他的原因就是你的程序有些语法压根就不对。自己查看ode23是怎么用的。。改成这样。
D=20:1:40;&&
for j=1:length(D)&&
& & d=D(j);
[t,y]=ode23 ('dbfun',[0:T/100:70*T],[4,4]);
plot(D(j),y(500:100:1400,2),'linewidth',5)&&
title('分岔图二')&&
xlabel(' 参数&&F0')&&
ylabel(' 位移')&&
function ydot=dbfun(t,y)&&
w=1;c=0.3;
ydot=[y(2);&&
(- y(1)^3 - c*y(2)+d*cos(w*t))/m];
<h1 style="color:# 麦片财富积分
zhanglu 发表于
其他的原因就是你的程序有些语法压根就不对。自己查看ode23是怎么用的。。改成这样。
D=20:1:40;&&
嗯，谢谢！！&&可是，上面那个用ODE45的程序怎么就出不来图呢，结果老是不对啊。而且ode23能自定义步长，那ode45应该也可以啊。能不能帮帮忙看下。我自己现在是实在看不出来了。。。& &麻烦了，谢谢！！
MATLAB 数学、统计与优化版块优秀回答者
关注者： 17
dairyqzy 发表于
嗯，谢谢！！&&可是，上面那个用ODE45的程序怎么就出不来图呢，结果老是不对啊。而且ode23能自定义步长， ...
我这边直接把23改为45是可以的，不知道你那边报错是什么？
<h1 style="color:# 麦片财富积分
zhanglu 发表于
我这边直接把23改为45是可以的，不知道你那边报错是什么？
应该是可以的啊，我看过有程序这么做过：动力学方程为
M 函数文件：dbfun.m
function ydot=dbfun(t,y,flag,d)
w=1;c=0.3;
ydot=[y(2);
- y(1)^3 - c*y(2)+d*cos(w*t)];
相应程序：
d=20:0.01:40;
for j=1:length(d)
[t,y]=ode23 ('dbfun',[0:T/100:70*T],[4,4],[],d(j));
plot(d(j),y(500:100:1400,2),'linewidth',5)
title('分岔图二')
xlabel(' 参数
ylabel(' 位移')
这里也是自定义步长的。
不好意思，刚才应该是我没表述清楚。上面这个程序我只是举个例子说明一下ode45是可以用自定义步长的。而我现在需要解决的问题是1楼的那4个非线性常微分方程组的问题。程序没有报错，但是做出来的分岔图与书上的不一致(就是我后面贴出来的那个图)。我试了很多次，但是都找不出原因，所以来求助您，希望您能帮忙看看。当然d的步长可以加大，t的范围可以缩小，这样仿真时间才不至于过长。麻烦了，谢谢！！
站长推荐 /3
Powered by