Runge-Kutta实现多自由度系统响应的MATLAB程序

yanyongju

function z = vbr_rk(rkf,u,t,x0)
%vbr_rk   vbr_rk(rkf,u,t,x0)
%       Runge-Kutta fourth order solution to a first order DE
%       t is a row vector from the initial time to the final time
%       step h.  
%       x0 is the initial value of the function.
%       The force matrix u is ordered such that the nth column of u is the force vector u evaluated at time n*dt.
%       rkf is a variable containing the name of the function file.  
%       The equations in the function file must be written in first order state space form.
%       See example vbr_rk_ex.m.
%       Example
%       t=0:.5:20;                             % Creates time vector
%       u=[zeros(1,length(t));sin(t*1.1)];     % Creates force matrix
%       x0=[1 ;0];                             % Creates initial state vector.
%       x=vtb1_3('vbr_rk_ex',u,t,x0);           % Runs analysis.
%       plot(t,x(1,:));                        % Plots displacement versus time.
%       plot(t,x(2,:));                        % Plots velocity versus time.

n=length(t);
z=zeros(length(x0),length(t));
z(:,1)=x0;
h=t(2)-t(1);

for l1=1:(n-1);
   z1=z(:,l1);
   u1=u(:,l1);
   u2=u(:,l1+1);

   k1=h*feval(rkf,z1,u1,t(l1));
   k2=h*feval(rkf,z1+.5*k1,u1,t(l1)+.5*h);
   k3=h*feval(rkf,z1+.5*k2,u1,t(l1)+.5*h);
   k4=h*feval(rkf,z1+k3,u1,t(l1)+h);
   z(:,l1+1)=z(:,l1)+1/6*(k1+2*k2+2*k3+k4);
end
其中rkf为动力微分方程，形如
function [zd]=vbr_rk_ex(z,u,t)
%    function for    2
%                  dx      dx
%                m --  + c -- +k x = f(t)
%                    2     dt
%                  dt                     dx
%    where m=1,k=1,c=.1, and z(1)=x, z(2)=--
%                                         dt
zd=[z(2);
    -.1*z(2)-z(1)+u(2)];

wei124

这个程序我没有看懂，是不是太特殊了。怎么没有地震波的输入

yanyongju

u表示激励输入
对于多自由度,关键在于rkf的形式

yanyongju

不客气！大家互相学习！

handb

这个程序怎么总出错?提示如下
??? Error: File: E:\MATLAB6p5p1\work\vbr_rk.m Line: 23 Column: 16
Missing variable or function.
怎么才能让程序计算呢?请解答,谢谢!

多情清秋

原帖由 handb 于 2007-3-30 15:54 发表
这个程序怎么总出错?提示如下
??? Error: File: E:\MATLAB6p5p1\work\vbr_rk.m Line: 23 Column: 16
Missing variable or function.
怎么才能让程序计算呢?请解答,谢谢!

应该是参数没给对，仔细研究一下程序的注释

handb

先谢谢多情清秋,但我用的是程序中的exemple啊???

handb

也就是这一段
%       Example
       t=0:.5:20;                             % Creates time vector
      u=[zeros(1,length(t));sin(t*1.1)];     % Creates force matrix
       x0=[1 ;0];                             % Creates initial state vector.
       x=vtb1_3('vbr_rk_ex',u,t,x0);           % Runs analysis.
%上一个语句无法调用子程序,我改为z=vbr_rk('vbr_rk_ex',u,t,x0); 否则提示找不到子程序         
       plot(t,z(1,:));                        % Plots displacement versus time.
      plot(t,z(2,:));                        % Plots velocity versus time.
结果出现我描述的错误，烦请高人指点，谢谢

WSYcxl

Runge-Kutta法是自启动的。
记得硕士论文答辩时，一个老师也这样问过：你这没有初始条件，难道系统就自己振动起来了？
当时一下子把我问蒙了，确实嘛，离开了激励，也没有初始条件，自己就振动起来了，没道理。呵呵。

yejet

正确的代码

1. t=0:.5:20;                             % Creates time vector
   
2. u=[zeros(1,length(t));sin(t*1.1)];     % Creates force matrix
   
3. x0=[1 ;0];                             % Creates initial state vector.
   
4. x=vbr_rk('vbr_rk_ex',u,t,x0);           % Runs analysis.
   
5. plot(t,z(1,:));                        % Plots displacement versus time.
   
6. plot(t,z(2,:));                        % Plots velocity versus time.

复制代码

建议楼主简单学习一下matlab，如果你要用matlab程序的话

[ 本帖最后由 无水1324 于 2007-7-23 13:29 编辑 ]

utuu

问下 ，如果是二自由度的动力学方程，那rkf该怎么构造呢？

博大广阔

我在找这个程序呢

S11021023

谢谢楼主！！！！！！！！！！

yghit08

自己编也不困难吧，还能自己设定很多参数，自己还明白怎么设定参数

0810064

每天来这里的目的就是为了学习，期待着有一天也能为这个论坛做一点贡献。