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

greatchina

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

主程序：
clear all;
clc
t0=0;tN=50;
x0=[0;0;0;0;0;0];
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=[x(1);x(2)];
x2=[x(3);x(4)];
x3=[x(5);x(6)];
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

有结果,但不对是不是?

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

greatchina

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

sigma665

原帖由 greatchina 于 2008-3-3 16:38 发表

                               
登录/注册后可看大图

鍔ㄥ姏瀛︽柟绋?

这是些什么字啊?

greatchina

请大家指点一下：

eight

原帖由 greatchina 于 2008-3-3 16:35 发表

                               
登录/注册后可看大图

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

主程 ...

请注意你的标题！发帖前，先阅读本版所有所有置顶帖！

eight

原帖由 greatchina 于 2008-3-3 16:42 发表

                               
登录/注册后可看大图

请大家指点一下：

问题的所有资料请尽量在一楼发完，不要分开三个楼层发。特别是“别人问一个你就提供一个”这种（虽然你不是这种情况），新手建议先看看版规