用荣格-库塔法，如何实现外部力为环境激励问题

huazi071783 · 发表于 2011-5-12 22:09

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

本帖最后由 huazi071783 于 2011-5-12 22:11 编辑

用荣格-库塔发解结构振动微分方程，求系统相应x，Mx''+cx'+kx=F(t),我设激振力F（t）为环境激励，就是用randn生成的随机数组。不知道能不能以我下面编的程序实现，就是逐个的解，虽然能计算，但计算速度特慢。如果激励为函数形势的如正弦函数就特别快。还有一下红字部分问题。下面是我的程序
clear;clc
global m k c f
m=0.5;
c=0.02;
k=80;
f=5;
x0=[0;0];
tspan=[0:0.01:100];
[t,y]=ode45(@randvibration,tspan,x0);
figure(1)
set(gca,'FontSize',30)
plot(y(2000:10000,1),y(2000:10000,2),'r')
调用函数
function dy=randvibration(t,y)
global m k c f
load('randsignal.mat');
F=randsignal;
n=floor(t./0.01);
%n=t./0.01;这样会出错，提示n到一定大的时候就不是整数了，按理说不会出现小数的。为什么？
% dy=[y(2);-(c./m)*y(2)-(k./m)*y(1)+f*cos(omeg*t); 这个计算特别块
dy(2)=-(c./m)*y(2)-(k./m)*y(1)+f*F(n+1);
dy(1)=y(2);
dy=dy';

结果不知道是否正确，请高手指点

huazi071783 · 发表于 2011-5-12 22:46

这是算出的结果画出的相轨迹图

meiyongyuandeze · 发表于 2011-5-13 08:23

本帖最后由 meiyongyuandeze 于 2011-5-13 08:25 编辑

不知道楼主程序里面为什么都是用的点除呢？程序运行慢是因为没此迭代都要读取数据的原因，可以尝试在主程序中先.mat数据后用参数传递的方法进入状态方程，这样可能会节省一定的计算时间。没有你的.mat文件，所以也没法运行，随便说说看法！

huazi071783 · 发表于 2011-5-13 10:15

回复 3 # meiyongyuandeze 的帖子

请问怎么用参数传递呢？没有明白你的意思？指教！

huazi071783 · 发表于 2011-5-13 14:16

回复 3 # meiyongyuandeze 的帖子

谢谢，解决了，果然是你说的那样，改完后计算非常快了

meiyongyuandeze · 发表于 2011-5-13 14:30

回复 5 # huazi071783 的帖子

呵呵，客气，一起学习！

账号		自动登录	找回密码
密码			我要加入

[分形与混沌] 用荣格-库塔法，如何实现外部力为环境激励问题

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。