如何求出外力为零时其对应的速度序列x'（i）

陆永杰

问题：如何求出（或者找出）外力为零时其对应的速度序列x'（i）
1
ode求出的常微分方程是数值解，现在可以求出系统的位移和速度的一系列数值解，可以画出相图（x-x’）、位移的时间历程（tau-x）、油膜力的时间历程（tau-fx、fy）、外力（即油膜力）和位移的相图（x-fx）。
2
问题是：
如何求出（或者找出）外力（也就是油膜力fx、fy，它们也是数值解，不一定恰好为零）为零时其对应的速度序列x'（i），求大家给指点一下。（是用枚举法遍历吗,matlab求解）
3
以下是常微分方程的.m文件和求解fx或fy对应速度序列的ode.m文件
a
常微分方程的.m文件
function
dxdt=force(t,x,flag,w)
r=0.005;p=0.004;L=0.008; c=0.001;m=0.4055;u=0.01;g=9.8;
a1=c/r;M=m*c*w*a1^2/(L*r*u);G1=g/c/w^2;
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);% x(1)=x,x(2)=x',x(3)=y,x(4)=y'
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);
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);
dxdt=zeros(4,1);
dxdt=[ x(2);

-fx/M+p*sin(t);

x(4);

-fy/M+p*cos(t)+G1];
b
求解fx或fy对应的速度序列的ode.m文件
n=1000;w=pi*n/30;x0=[0.42;0;0.21;0];T=2*pi/w;
options=odeset('RelTol',1.e-7,'AbsTol',[1.e-8 1.e-8 1.e-8 1.e-8]);
[tm,x]=ode15s('force',[0:T/10:1*T],x0,[],w);
% x1=x(:,1);x2=x(:,2);x3=x(:,3);x4=x(:,4);
m1=x(:,1);m2=x(:,2);m3=x(:,3);m4=x(:,4);%用于计算动能差序列
x1=m1;x2=m2;x3=m3;x4=m4;
r=0.005;L=0.004; c=0.001;m=0.4055;p=0.001;u=0.0178;g=9.8;
a1=c/r;M=m*c*w*a1^2/(L*r*u);G1=g/c/w^2;
a=atan((x3+2*x2)./(x1-2*x4))-pi/2*sign((x3+2*x2)./(x1-2*x4))-pi/2*sign(x3+2*x2);
G=2./sqrt(1-x1.^2-x3.^2).*(pi/2+atan((x3.*cos(a)-x1.*sin(a))./sqrt(1-x1.^2-x3.^2)));
S=(x1.*cos(a)+x3.*sin(a))./(1-(x1.*cos(a)+x3.*sin(a)).^2);
% x(1)=x,x(2)=x',x(3)=y,x(4)=y'
V=(2+(x3.*cos(a)-x1.*sin(a)).*G)./(1-x1.^2-x3.^2);
fx=sqrt((x1-2*x4).^2+(x3+2*x2).^2)./(1-x1.^2-x3.^2).*(3*x1.*V-sin(a).*G-2.*cos(a).*S);
fy=sqrt((x1-2*x4).^2+(x3+2*x2).^2)./(1-x1.^2-x3.^2).*(3*x3.*V+cos(a).*G-2.*sin(a).*S);

y1=[]; y2=[];
for i=1:length(x1)

for j=1:length(x2)

for k=1:length(x3)

for l=1:length(x4)

if ( a=atan((x3(k)+2*x2(j))./(x1(i)-2*x4(l)))-pi/2*sign((x3(k)+2*x2(j))./(x1(i)-2*x4(l)))-pi/2*sign(x3(k)+2*x2(j)) ;

and
G=2./sqrt(1-x1(i).^2-x3(k).^2).*(pi/2+atan((x3(k).*cos(a)-x1.*sin(a))./sqrt(1-x1(i).^2-x3(k).^2)));


and S=(x1(i).*cos(a)+x3(k).*sin(a))./(1-(x1(i).*cos(a)+x3(k).*sin(a)).^2);
% x(1)=x,x(2)=x',x(3)=y,x(4)=y'

and V=(2+(x3(k).*cos(a)-x1(i).*sin(a)).*G)./(1-x1(i).^2-x3(k).^2);

sqrt((x1(i)-2*x4(l)).^2+(x3(k)+2*x2(j)).^2)./(1-x1(i).^2-x3(k).^2).*(3*x1(i).*V-sin(a).*G-2.*cos(a).*S)<1e-3)

y1=[y1;i,j,k,l];

y2=[y2;x1(i),x2(j),x3(k),x4(l)];

else y1=y1; y2=y2;

end

end

end

end

end
save('results.mat','y1','y2')
红色部分如何处理才可找出fx为0时其对应的时刻（tau），从而进一步求出其对应的x'（x2序列）

sogooda

问题描述好复杂啊，看了半天才好想看明白了。
就是说对应一个序列x=0时，另一个序列y的值吧
可以先求x中最接近0的两个值，求出对应于x的y，再进行插值。