回复 15楼 的帖子
呵呵！难怪做不出来。那是不是得按非自治系统分叉图的方法做呀？我按照论坛上面的帖子重新编了主程序：clc;
clear;
global M w
m=2;
C=0.0001;
eta=0.018;
R=0.025;
L=0.01;
w=M*eta*R*L*(R/C)^2*(L/(2*R))^2/(m*C);
range=10:300;
T=2*pi/w;
k=0;
yy1=[];
step=pi/100;
for M=range
x0=;
k=k+1;
tspan=;
=ode45('youmofun',tspan,x0);
x0=y(end,:);
j=1;
for i=1000:1200
tspan=;
=ode45('youmofun',tspan,x0);
yy1(k j)=y(end,1);
j=j+1;
x0=y(end,:);
end
end
bifdata=yy1(:,1000:end);
plot(M,bifdata,'k.','markersize',1);
运行出现如下错误：
??? Error using ==> mrdivide
Matrix dimensions must agree.
请问这样做合适吗？ ode45('youmofun',tspan,x0);
这句话调用的不对哦，呵！参考论坛里面的分岔图程序吧！
回复 17楼 的帖子
function dx=youmofun(t,x)global M w
m=2;
C=0.0001;
eta=0.018;
R=0.025;
L=0.01;
theta=0;
rho=0.1;
g=9.8;
w=M*eta*R*L*(R/C)^2*(L/(2*R))^2/(m*C);
a=atan((x(3)+2*x(2))/(x(1)-2*x(4)))-pi/2*sign((x(3)+2*x(2))/(x(1)-2*x(4)))-pi/2*sign(x(3)+2*x(2));
G=2/sqrt(1-x(1)^2-x(3)^2)*(pi/2+atan((x(3)*cos(a)-x(1)*sin(a))/sqrt(1-x(1)^2-x(3)^2)));
S=(x(1)*cos(a)+x(3)*sin(a))/(1-(x(1)*cos(a)+x(3)*sin(a))^2);
V=(2+(x(3)*cos(a)-x(1)*sin(a))*G)/(1-x(1)^2-x(3)^2);
fx=sqrt((x(1)-2*x(4))^2+(x(3)+2*x(2))^2)/(1-x(1)^2-x(3)^2)*(3*x(1)*V-sin(a)*G-2*cos(a)*S)+8*pi*R*C/((L^2)*sqrt(x(1)^2+x(3)^2))*(1-1/sqrt(1-x(1)^2-x(3)^2))*cos(a);
fy=sqrt((x(1)-2*x(4))^2+(x(3)+2*x(2))^2)/(1-x(1)^2-x(3)^2)*(3*x(3)*V+cos(a)*G-2*sin(a)*S)+8*pi*R*C/((L^2)*sqrt(x(1)^2+x(3)^2))*(1-1/sqrt(1-x(1)^2-x(3)^2))*sin(a);
w*t+theta=x(5);
dx=zeros(5,1);
dx(1)=x(2);
dx(2)=fx/M+rho*cos(w*t+theta);
dx(3)=x(4);
dx(4)=fx/M+rho*cos(w*t+theta)-g/C*w^2;
dx(5)=w;
我把方程的程序改成上面的那样，运行的时候还是提示矩阵维数不一致的错误，:@L 请问还是方程程序写的有问题吗？主程序和论坛上面的主程序一样。 w*t+theta=x(5);
这是什么意思 我是按照论坛上面的程序写的:@)
它是轮盘的偏心力在x，y上面的分量
具体的无量纲方程如下：
回复 15楼 的帖子
院长开始求神了，呵呵。开个玩笑的。
页:
1
[2]