[讨论]参数激励非线性系统分析

xiaojin831108

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方向的分量表达式

[fx1,fy1]=youmoli(x(1),x(7),dx(1),dx(7))
[fx2,fy2]=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

xiaojin831108

请各位高手能不能帮我看下这个程序，还有个youmoli函数我也发上去，急死了都调好几天了，一运行就要几个小时不知道怎么回事？
function [fx,fy]=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)

shenyongjun

ode45就是变步长的RK法。
你的这种调用方式只不过是在输出时按照时间等间隔地输出计算结果，但是其中的计算过程是定步长的。
可以 help ode45，查看说明

xujy

请问你的频率响应曲线是怎么得到的，是用的什么方法，我现在遇到个方程和你的差不多，不知怎么求解，麻烦给指导下，谢谢！

shenyongjun

见附件，在Mathworks公司主页上找到的。
下载压缩文件后将扩展名中的txt去掉解压即可

无水1324

谢谢！
原来是这样的，所以我见到很多书上都把这种叫做定步长了！

还有我i用这种方法时总是出现周期解不收敛（就是说Poincare 图上出现的点相差很小10-5左右），是不是因为我数出时间取等间隔的缘故？

[ 本帖最后由 无水1324 于 2007-3-30 08:32 编辑 ]

无水1324

可以用谐波平衡，多尺度法求，并得到其频率响应方程。然后就可以画出来了！

xiaojin831108

你们帮我看下那个解方程的程序好吗？我怎么就是调试不出结果，每选一个初值就要运行很久，好急的不知道是方法的问题还是程序的问题，我是在matlab中直接调用的ode45函数

shenyongjun

跟等间隔取数据点没有关系。出现你说的这种情况，一个原因就是你的系统阻尼较小，收敛速度太慢；另一个原因就是系统处于混沌状态。

shenyongjun

你的程序是不是有问题？运行时出现x未定义的提示：

??? Input argument "x" is undefined.

Error in ==> sdzz1 at 20
if x(10)>dalta1

xiaojin831108

要是系统处于混沌状态的话就可以直接做分叉图或poincare图了吗？就是不解响应可以吗？
   
申老师，我是在matlab中运行ode45时给x赋的初值，这样应该就不会说x没定义了，还有就是这个系统是不是和初值关系很大啊？
  我是做的书上的松动转子系统的仿真，参数也是按他们给定的参数，这个好像东北大学的闻邦春老师和他的学生做的很多，他们说是用4阶runge-kutta解微分方程组.
  现在就是调试不出来，不知道是不是程序哪有问题还是我本来理解的就不对！谢谢你们帮我解决问题！！

无水1324

解响应之后，才能够画出分岔图和poincare图

yejet

原帖由 xiaojin831108 于 2007-3-30 12:09 发表
要是系统处于混沌状态的话就可以直接做分叉图或poincare图了吗？就是不解响应可以吗？
   
申老师，我是在matlab中运行ode45时给x赋的初值，这样应该就不会说x没定义了，还有就是这个系统是不是和初值关系很大 ...

rk法迭代需要一定的收敛过程，所以前面迭代出来的结果是没有收敛的结果

对于振动问题而言，大家更多的是关心稳态解，因此迭代后舍弃前面若干步的解，从收敛后取值

xiaojin831108

我不知道上面我发的程序哪里出了问题？？？每次都要运行好久，参数我是按照书上面的给的，不知道是不是解微分方程的程序有问题，麻烦高手们有空帮我看下好不？特别感谢，很急！！！！

shenyongjun

这样吧，你把要仿真的方程给出来，其中的参数也给出来，我看能不能帮你解决这个问题。因为每个人的编程习惯都不一样，看你的这个程序太费劲了！