请教用荣格库塔法解结构振动微分方程，程序问题，结果不正确

huazi071783

用荣格-库塔发解结构振动微分方程，Mx''+cx'+kx=F(t),我设激振力为正弦周期力F（t）=f*cos(omeg*t),变化后：x'=y;dy=-(c./m)*y-(k./m)*x+F(t);求解后画出结果的相轨迹图，但是得到的结果矩阵y的第一列全为零，第二列是周期值看似正确，为什么第一列为零呢？这结果不对啊！不知道什么原因，请高手解答，谢谢！
clear;clc
global m omeg k c f
m=0.5;                 %质量
c=0.2;                   %阻尼
k=1;                      %刚度
omeg=6.3;
f=5;
x0=[0;0];               %初始值
tspan=[0:0.01:100];
[t,y]=ode45('vibration',tspan,x0);      %y为结果矩阵
plotmatrix (y);  figure(gcf)
figure(2)
plot(y(:,1),y(:,2))                %相轨迹图

调用结构振动函数
function dy=vibration(t,y)
global m omeg k c f
dy=[y(1);-(c./m)*y(2)-(k./m)*y(1)+f*cos(omeg*t)];

huazi071783

本来得到的结果y（:,1）应该是周期震荡值，但是结果全是0，怎么解释呢？

hsfy919

回复 2 # huazi071783 的帖子

应该是初值的问题吧，你用[0.1,0]试试

huazi071783

回复 3 # hsfy919 的帖子

主任啊，我已经试过了，y(:,1)是单调递增的，也不是震荡值，不对啊，晕了

huazi071783

还没有解决，高手帮我看看是什么问题

hsfy919

回复 5 # huazi071783 的帖子

是你的方程输错了“dy=[y(2);-(c./m)*y(2)-(k./m)*y(1)+f*cos(omeg*t)];”，再检查一下

huazi071783

回复 6 # hsfy919 的帖子

谢谢主任，这样写的确结果不为0了，我这个是故意这样写的，在方程中y(2)代表的是x',y(1)代表的是原方程的x，如果结果写成dy=[y(2);-(c./m)*y(2)-(k./m)*y(1)+f*cos(omeg*t)]时得到的结果第一列应该是x’，第二列应该是dx'也就是dy。如果结果写成dy=[y(1);-(c./m)*y(2)-(k./m)*y(1)+f*cos(omeg*t)]，结果第一列应该是x，第二列应该是dy。不知道是不是这样？我是想输出解方程后x的值，然后再求一次解出y值，这样可作出相轨迹图来。呵呵，指教

hsfy919

回复 7 # huazi071783 的帖子

不是的，利用龙格-库塔方法求出的数值解是对y的迭代，其结果是[y(1);y(2)]的排列，直接就能画相轨迹图

huazi071783

回复 8 # hsfy919 的帖子

还是没有明白，结构振动响应应该是x的时间序列，相轨迹也应该是基于x时间序列得出。用龙格-库塔发求解得到的结果，像你说的那样按[y(1),y(2)]排列，那么这结果各代表什么？第一列y(1)是y，也就是x'？y(2)是dy？还是第二列y(2)是基于y(1)重构后的时间延迟后的y值?我看了matlab的help也没明白，呵呵，谢主任

hsfy919

回复 9 # huazi071783 的帖子

我习惯把未知量用y表示，第一列是y(1)(也就是x),第二列是y(2)(也就是x')，没有dy。

huazi071783

回复 10 # hsfy919 的帖子

主任，非常感谢，终于弄懂了