非线性振动程序有错，三天未果，大家来看看怎么回事？？
%利用逐步积分法求非线性系统的响应% 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();
velocity=zeros();
accelaration=zeros();
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) 原来我估计是由于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 编辑 ] 同时k的表达式是我从参考文献上引用的，别人都做出来了，这是为什么啊？
请大家见附件的参考文献
[ 本帖最后由 mcx8305 于 2007-7-24 11:13 编辑 ] 程序算出k(1)=NaN,但我从第一步直接得出的k=???e6(好像是这个数量级，具体的值我忘了)
大家根据程序、图、参考文献看看这到底是什么问题。
也许问题很简单，但却有代表性，欢迎讨论，也帮我出出主意~~论文用。 花大侠，你帮我看看吧~~ 出现响应突然无限增大可麻烦了，你这个程序按照直接积分得来的，因为一时还很难看出来问题，可能是你积分中推导有问题，或者这个问题不适合直接积分，为什么不用ode45求解呢，应该更简单
回复 #5 mcx8305 的帖子
谢谢你的信任，不过很惭愧，图像方面我还凑合。数学上就不行，基础不好。以后不要指名哦，知道的都会回答你的:handshake ，多亏女侠不计较还帮你出了注意。 这个问题最好还是用ode23tb求解，具体使用可以看一下help！ 由于所求问题的需要，后续还将分析非线性微分方程组，故选择了直接积分法。 谢谢各位大侠相助，我再看看，有情况会传上来。欢迎大家继续交流，集思广益！！向花大侠致歉！
页:
[1]