贴出我的转子运动程序,在用poincare方法求它的稳定点时,这个方程的jacobian是什么呢?

胡晓宇

function xprime=odetest_sub(t,s)

%m1 m2转子1,2的质量;e1 e2为转子的偏心量;c1x为轴承1在x方向阻尼;g重力加速度
%w1,w2为转子的转速 Nb1为轴承1的转子数量 gama为游隙
%s(1)=s1     x1 转子1在x方向的位移 垂直方向
%s(2)=s2     dx1 转子1在x方向速度
%s(3)=s3     y1 转子1在y方向位移    水平方向
%s(4)=s4     dy1 转子1在y方向速度


global m e c1x c1y c2x c2y g F w Nb1 Nb2 gama1 gama2...
     r1 R1 r2 R2 kn1 kn2

xprime=zeros(4,1);

%*******************************计算轴承力*********************************
%轴承1深沟球轴承的轴承力fq1(2:1) qx1 qy1 转子1在x，y方向位移s(1),s(3)
%轴承1中心在x，y方向的位移分别为s(1),s(3) 其轴承力x，y方向分别为fq1
fq1=zeros(2,1);
wc1=(w*r1)/(r1+R1);%滚子公转转速 轴承1的内圈转速为w1，外圈为转子1的转速0
for No_ball1=1:Nb1
sita(No_ball1)=2*pi/Nb1*(No_ball1-1)+wc1*t;% 第No_ball个滚珠的位置角
Clearance1(No_ball1,1)=s(3)*sin(sita(No_ball1) ) ...
+ s(1)*cos( sita(No_ball1) )-gama1;% 滚珠与内滚道的间隙变化。

if Clearance1(No_ball1)<=0;% 判断哪几个滚动体受到接触力
Clearance1(No_ball1) = 0;
end
fs1 = abs((Clearance1(No_ball1))^1.5 );        
fq1(1,1) = fq1(1,1)+kn1*fs1*cos(sita(No_ball1));% x方向轴承力
fq1(2,1) = fq1(2,1)+kn1*fs1*sin(sita(No_ball1));% y方向轴承力
end

%轴承2深沟球轴承的轴承力fq2(2:1) qx2 qy2 转子1在x，y方向位移s(1),s(3)
%轴承2中心在x，y方向的位移分别为s(1),s(3) 其轴承力x，y方向分别为fq2
fq2=zeros(2,1);
wc2=(w*r2)/(r2+R2);%滚子公转转速，内圈转速为w1，外圈为转子1的转速0
for No_ball2=1:Nb2
sita(No_ball2)=2*pi/Nb2*(No_ball2-1)+wc2*t;% 第No_ball个滚珠的位置角
Clearance2(No_ball2,1)=s(3)*sin( sita(No_ball2) ) ...
+ s(1)*cos(sita(No_ball2))-gama2;% 滚珠与内滚道的间隙变化。

if Clearance2(No_ball2)<=0;% 判断哪几个滚动体受到接触力
Clearance2(No_ball2) = 0;
end
fs2 = abs((Clearance2(No_ball2))^1.5 );        %??????????????????????
fq2(1,1) = fq2(1,1)+kn2*fs2*cos(sita(No_ball2));% x方向轴承力
fq2(2,1) = fq2(2,1)+kn2*fs2*sin(sita(No_ball2));% y方向轴承力
end

%*******************************转子运动的微分方程*********************************
xprime(1)=s(2);
%ds1=s2
xprime(2)=e*w^2.*cos(w*t)+g+F/m-c1x*s(2)/m-fq1(1,1)/m-fq2(1,1)/m;
%ds2
xprime(3)=s(4);
%ds3=s4
xprime(4)=e*w^2.*sin(w*t)-c1y*s(4)/m-fq1(2,1)/m-fq2(2,1)/m;
%ds4
-------------------------------------------------------------------------------------------
方程如上,转子运动由四个方程组成,fq1(1,1),fq2(1,1),fq1(2,1),fq2(2,1)分别是轴承力,由程序看,都是分段函数,在用poincare映射法求不动点时,需要求poincare算子的jacobian,我的方程jacobian怎么求呢?分别对变量s(i)求偏导？可里面还有fq1(1,1),fq2(1,1),fq1(2,1),fq2(2,1)分段函数，郁闷中。盼学长们给予解答。

胡晓宇

怎么没有人回答呢，盼大家浏览后给点意见呀

cam_1980

那你的雅克比矩阵也应该是分段的，分段来求解雅克比矩阵，当自变量属于哪一段的时候，这是就将对应的函数表达式代入微分方程右边，然后再求解雅克比矩阵

胡晓宇

多谢cam_1980
还有个问题,我求jacobian矩阵,是求poincare算子的. 说白了,就是在用newton_raphson和poincare映射求系统稳定点,即求P(x)-x的零点.P为poincare算子.newton_raphson需要求P(x)-x的jacobian矩阵，你知道怎么求么？盼告知。再次感谢！

无水1324

newton_raphson需要求P(x)-x的jacobian矩阵，你知道怎么求

P(x)-x当这个一个多元函数，
对每一个元素求导组成一个矩阵，
比如说：
f1=x1^2+x2^3+x2*x1
f2=x1+x2*x1^2
它的jacobian矩阵就是
[2x1+x2, 3x2^2+x1;
1+2x1x2,x1^2]
这个就是一个jacobian矩阵了

胡晓宇

先对无水表示衷心的感谢！
以您列出的系统为例
f1=x1^2+x2^3+x2*x1
f2=x1+x2*x1^2
要求其稳定点，通常求的是P(x)-x的零点。则需要求P(x)的jacobian，P是庞加莱映射算子。而您给出的矩阵是向量场f的jacobian，是不是针对您的系统，可以用f代替P(x)呢？盼告知！针对您的系统，它的庞加莱算子是什么呢？

无水1324

应该使用f代替P(x)-x，这里面的x不会就是一个量吧？

胡晓宇

x就是参数组成的向量，这个例子中应该是x=[x1，x2]
用f代替P(x)-x,有理论依据么？P(x)是x的庞加莱映射

sizhiyuan2006

学习了，谢谢

mxlzhenzhu

不是什么系统都有导数的。