greatchina 发表于 2008-3-3 16:35

R-K方法求解3自由度系统振动方程时遇到不解

刚接触MATLAB时间不长,在编写下面的R-K方法求解3自由度系统振动方程时发现,结果总和所给标准答案不符,多选择了几组初值,也是不一样,不知道是什么原因,请各位帮忙给看看,是不是程序本身就有问题,谢谢

主程序:
clear all;
clc
t0=0;tN=50;
x0=;
h=0.01;
t=t0:h:tN;N =length(t);j=1;
for i=1:N
t1=t0+h;
K1=m_chap2_ex1_1_sub(t0,x0);
K2=m_chap2_ex1_1_sub(t0+h/2,x0+h*K1'/2);
K3=m_chap2_ex1_1_sub(t0+h/2,x0+h*K2'/2);
K4=m_chap2_ex1_1_sub(t0+h,x0+h*K3');
x=x0+(h/6)*(K1'+2*K2'+2*K3'+K4');
yy1(j)=x(1);yy2(j)=x(2);yy3(j)=x(3);yy4(j)=x(4);yy5(j)=x(5);yy6(j)=x(6);
t0=t1;x0=x;j=j+1;
end
t=0:h:tN;
plot (t,yy2); grid on

子程序:
function ydot=vdpol(t,x)
x1=;
x2=;
x3=;
m1=1;
m2=1;
m3=1;
KK1=2;KK2=2;KK3=1;KK4=2;
p0=1;
w=3;
ydot(1)=(-KK1*x1(2)-KK2*(x1(2)-x2(2))+p0*sin(w*t))/m1;
ydot(2)=x1(1);
ydot(3)=(KK2*(x1(2)-x2(2))-KK3*(x2(2)-x3(2)))/m2;
ydot(4)=x2(1);
ydot(5)=(KK3*(x2(2)-x3(2))-KK4*x3(2))/m3;
ydot(6)=x3(1);

[ 本帖最后由 eight 于 2008-3-3 17:16 编辑 ]

sigma665 发表于 2008-3-3 16:37

回复 楼主 的帖子

有结果,但不对是不是?

建议断点调试,一行行运行,看那边开始值不对

greatchina 发表于 2008-3-3 16:38

方程



[ 本帖最后由 eight 于 2008-3-3 17:15 编辑 ]

sigma665 发表于 2008-3-3 16:39

原帖由 greatchina 于 2008-3-3 16:38 发表 http://www.chinavib.com/forum/images/common/back.gif
鍔ㄥ姏瀛︽柟绋?
这是些什么字啊?

greatchina 发表于 2008-3-3 16:42

请大家指点一下:

eight 发表于 2008-3-3 17:16

原帖由 greatchina 于 2008-3-3 16:35 发表 http://www.chinavib.com/forum/images/common/back.gif
刚接触MATLAB时间不长,在编写下面的R-K方法求解3自由度系统振动方程时发现,结果总和所给标准答案不符,多选择了几组初值,也是不一样,不知道是什么原因,请各位帮忙给看看,是不是程序本身就有问题,谢谢

主程 ...
请注意你的标题!发帖前,先阅读本版所有所有置顶帖!

eight 发表于 2008-3-3 17:19

原帖由 greatchina 于 2008-3-3 16:42 发表 http://www.chinavib.com/forum/images/common/back.gif
请大家指点一下:
问题的所有资料请尽量在一楼发完,不要分开三个楼层发。特别是“别人问一个你就提供一个”这种(虽然你不是这种情况),新手建议先看看版规
页: [1]
查看完整版本: R-K方法求解3自由度系统振动方程时遇到不解