非线性振动程序有错，三天未果，大家来看看怎么回事？？

mcx8305

%利用逐步积分法求非线性系统的响应
% 9300*y''+5*y'+2*73397*(1-((y+sign(y')*85e-6)/(2*190e-6))^2)*((y+sign(y')*85e-6)/2)-(73397*(1-(y/190e-6)^2)*y)*sign(y')=5*cos(30*t)
%clc;    %Clear Command Window
clear all;
close all;
deltat=0.002;
t=0:deltat:6;
p=5*cos(30*t);
% initial the displacement and velocity
y=0.1;v=1;m=9300;c=5;
len=length(t)-1;
displacement=zeros([1,len]);
velocity=zeros([1,len]);
accelaration=zeros([1,len]);
for i=1:len

if y>=0&y<0.001
y=y+0.001;
else if y<0&y>-0.001
y=y-0.001;
else
y=y;
end
end

k=(2*73397*(1-((y+sign(v)*85e-4)/(2*190e-4))^2)*((y+sign(v)*85e-4)/2)-(73397*(1-(y/190e-4)^2)*y)*sign(v))/y;
a=(p(i)-c*v-2*73397*(1-((y+sign(v)*85e-4)/(2*190e-4))^2)*((y+sign(v)*85e-4)/2)-(73397*(1-(y/190e-4)^2)*y)*sign(v))/m;
equivelentk=k+6*m/(deltat^2)+3/deltat*c;
deltap=(p(i+1)-p(i))+m*(6/deltat*v+3*a)+c*(3*v+deltat/2*a);
deltay=deltap/equivelentk;
deltav=3/deltat*deltay-3*v-deltat/2*a;
y=y+deltay;
v=v+deltav;
displacement(i)=y;
velocity(i)=v;
accelaration(i)=a;

k(i)=k;

end
figure(1)
plot(t(1:len),displacement)
figure(2)
plot(t(1:len),velocity)
figure(3)
plot(t(1:len),accelaration)

figure(4)
plot(t(1:len),k)

k(1)

mcx8305

原来我估计是由于y太小，导致函数k出了问题，因此把下列程序加上了：
if y>=0&y<0.001
y=y+0.001;
else if y<0&y>-0.001
y=y-0.001;
else
y=y;
end
end
但是问题依然存在。
显示的k(1)=NaN,而其他的k(i)都是0，如图4所示。

[ 本帖最后由 mcx8305 于 2007-7-24 11:08 编辑 ]

mcx8305

同时k的表达式是我从参考文献上引用的，别人都做出来了，这是为什么啊？
请大家见附件的参考文献

[ 本帖最后由 mcx8305 于 2007-7-24 11:13 编辑 ]

mcx8305

程序算出k(1)=NaN,  但我从第一步直接得出的k=???e6(好像是这个数量级，具体的值我忘了)
大家根据程序、图、参考文献看看这到底是什么问题。

也许问题很简单，但却有代表性，欢迎讨论，也帮我出出主意~~论文用。

mcx8305

花大侠，你帮我看看吧~~

咕噜噜

出现响应突然无限增大可麻烦了，你这个程序按照直接积分得来的，因为一时还很难看出来问题，可能是你积分中推导有问题，或者这个问题不适合直接积分，为什么不用ode45求解呢，应该更简单

花如月

谢谢你的信任，不过很惭愧，图像方面我还凑合。数学上就不行，基础不好。以后不要指名哦，知道的都会回答你的:handshake ，多亏女侠不计较还帮你出了注意。

octopussheng

这个问题最好还是用ode23tb求解，具体使用可以看一下help！

mcx8305

由于所求问题的需要，后续还将分析非线性微分方程组，故选择了直接积分法。

mcx8305

谢谢各位大侠相助，我再看看，有情况会传上来。欢迎大家继续交流，集思广益！！
向花大侠致歉！