帮忙看看我的这个程序行吗？楼主和各位高手们，怎么运行下来就要几个小时啊？
function dx=sdzz1(t,x)dx=zeros(14,1)
%参数设置
m1=4
m2=32.1
m3=50
R=25
L=12
c=0.11
miu=0.018
c1=1050
c2=2100
omiga=25
k=2.5*10^7
dalta1=0.01
ks1=7.5*10^7
ks2=2.5*10^9
cs1=350
cs2=500
if x(10)>dalta1
ks=ks1
cs=cs1
elseif (x(10)>=0)&(x(10)<=dalta1 )
ks=0
cs=0
else
ks=ks2
cs=cs2
end
b=0.02
P=m2/2
delta=((miu*omiga*R*L)/P)*(R/c)^2*((L/2*R)^2)
pi=3.14159
%无量纲化
epus1=c1/(m1*omiga)
epus3=c1/(m1*omiga)
epus2=c2/(m2*omiga)
epus4=cs/(m3*omiga)
eta1=k/(m1*omiga^2)
eta3=k/(m1*omiga^2)
eta2=k/(m2*omiga^2)
eta4=ks/(m3*omiga^2)
M1=(m1*c*omiga^2)/(delta*P)
M3=(m1*c*omiga^2)/(delta*P)
M4=(m3*c*omiga^2)/(delta*P)
g=9.8;
G=g/(omiga^2)
%微分方程组程序
dx(1)=x(4)
dx(2)=x(5)
dx(3)=x(6)
%dy(1)=x(11)%y(5)
dx(7)=x(11)
%dy(2)=x(12)%y(6)
dx(8)=x(12)
%dy(3)=x(13)%y(7)
dx(9)=x(13)
%dy(4)=x(14)%y(8)
dx(10)=x(14)
%油膜力x,y方向的分量表达式
=youmoli(x(1),x(7),dx(1),dx(7))
=youmoli(x(3),x(9)-x(10),dx(3),dx(9)-dx(10))
%继续
dx(4)=-epus1*x(4)-eta1*(x(1)-x(2))+1/M1*fx1
dx(11)=-epus1*x(11)-eta1*(x(7)-x(8))+1/M1*fy1
dx(5)=-epus2*x(5)-eta2*(2*x(2)-x(1)-x(3))+b*cos(omiga*t)
dx(12)=-epus2*x(12)-eta2*(x(8)-x(7)-x(9))+b*sin(omiga*t)-G
dx(6)=-epus3*x(6)-eta3*(x(3)-x(2))+1/M3*fx2
dx(13)=-epus3*x(13)-eta3*(x(9)-x(8))+1/M3*fy2-G
dx(14)=-epus4*x(14)-eta4*x(10)-1/M4*fy2-G 请各位高手能不能帮我看下这个程序，还有个youmoli函数我也发上去，急死了都调好几天了，一运行就要几个小时不知道怎么回事？
function =youmoli(x,y,dx,dy)
alpha=atan((y+2*dx)/(x-2*dy))-pi/2*sign((y+2*dx)/(x-2*dy))-pi/2*sign(y+2*dx);
S=(x*cos(alpha)+y*sin(alpha))/(1-(x*cos(alpha)+y*sin(alpha)^2));
L=(2/(sqrt(1-x^2-y^2)))*(pi/2+atan(((y*cos(alpha)-x*sin(alpha))/sqrt(1-x^2-y^2))))
V=(2+(y*cos(alpha)-x*sin(alpha))*L)/(1-x^2-y^2)
fx=(sqrt((x-2*dy)+(y+2*dx)^2)/(1-x^2-y^2))*(3*x*V-sin(alpha)*L-2*cos(alpha)*S)
fy=(sqrt((x-2*dy)+(y+2*dx)^2)/(1-x^2-y^2))*(3*y*V-cos(alpha)*L-2*sin(alpha)*S)
回复 #56 无水1324 的帖子
ode45就是变步长的RK法。你的这种调用方式只不过是在输出时按照时间等间隔地输出计算结果，但是其中的计算过程是定步长的。
可以 help ode45，查看说明
回复 #22 无水1324 的帖子
请问你的频率响应曲线是怎么得到的，是用的什么方法，我现在遇到个方程和你的差不多，不知怎么求解，麻烦给指导下，谢谢！回复 #57 xiaojin831108 的帖子
见附件，在Mathworks公司主页上找到的。下载压缩文件后将扩展名中的txt去掉解压即可
回复 #63 shenyongjun 的帖子
谢谢！原来是这样的，所以我见到很多书上都把这种叫做定步长了！
还有我i用这种方法时总是出现周期解不收敛（就是说Poincare 图上出现的点相差很小10-5左右），是不是因为我数出时间取等间隔的缘故？
[ 本帖最后由 无水1324 于 2007-3-30 08:32 编辑 ]
回复 #64 xujy 的帖子
可以用谐波平衡，多尺度法求，并得到其频率响应方程。然后就可以画出来了！ 你们帮我看下那个解方程的程序好吗？我怎么就是调试不出结果，每选一个初值就要运行很久，好急的不知道是方法的问题还是程序的问题，我是在matlab中直接调用的ode45函数回复 #66 无水1324 的帖子
跟等间隔取数据点没有关系。出现你说的这种情况，一个原因就是你的系统阻尼较小，收敛速度太慢；另一个原因就是系统处于混沌状态。回复 #68 xiaojin831108 的帖子
你的程序是不是有问题？运行时出现x未定义的提示：??? Input argument "x" is undefined.
Error in ==> sdzz1 at 20
if x(10)>dalta1
回复shenyongjun老师！！
要是系统处于混沌状态的话就可以直接做分叉图或poincare图了吗？就是不解响应可以吗？申老师，我是在matlab中运行ode45时给x赋的初值，这样应该就不会说x没定义了，还有就是这个系统是不是和初值关系很大啊？
我是做的书上的松动转子系统的仿真，参数也是按他们给定的参数，这个好像东北大学的闻邦春老师和他的学生做的很多，他们说是用4阶runge-kutta解微分方程组.
现在就是调试不出来，不知道是不是程序哪有问题还是我本来理解的就不对！谢谢你们帮我解决问题！！ 解响应之后，才能够画出分岔图和poincare图 原帖由 xiaojin831108 于 2007-3-30 12:09 发表
要是系统处于混沌状态的话就可以直接做分叉图或poincare图了吗？就是不解响应可以吗？
申老师，我是在matlab中运行ode45时给x赋的初值，这样应该就不会说x没定义了，还有就是这个系统是不是和初值关系很大 ...
rk法迭代需要一定的收敛过程，所以前面迭代出来的结果是没有收敛的结果
对于振动问题而言，大家更多的是关心稳态解，因此迭代后舍弃前面若干步的解，从收敛后取值