huazi071783 发表于 2011-3-7 11:06

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

本帖最后由 huazi071783 于 2011-3-7 16:23 编辑

用荣格-库塔发解结构振动微分方程，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=;               %初始值
tspan=;
=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=;

huazi071783 发表于 2011-3-7 13:56

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

hsfy919 发表于 2011-3-7 23:00

回复 2 # huazi071783 的帖子

应该是初值的问题吧，你用试试

huazi071783 发表于 2011-3-8 10:20

回复 3 # hsfy919 的帖子

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

huazi071783 发表于 2011-3-8 16:38

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

hsfy919 发表于 2011-3-9 18:45

本帖最后由 hsfy919 于 2011-3-9 18:45 编辑

回复 5 # huazi071783 的帖子

是你的方程输错了“dy=;”，再检查一下

huazi071783 发表于 2011-3-10 08:57

回复 6 # hsfy919 的帖子

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

hsfy919 发表于 2011-3-10 10:39

回复 7 # huazi071783 的帖子

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

huazi071783 发表于 2011-3-10 13:29

回复 8 # hsfy919 的帖子

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

hsfy919 发表于 2011-3-10 20:07

回复 9 # huazi071783 的帖子

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

huazi071783 发表于 2011-3-10 21:49

回复 10 # hsfy919 的帖子

主任，非常感谢，终于弄懂了
页: [1]
查看完整版本: 请教用荣格库塔法解结构振动微分方程，程序问题，结果不正确