『转贴』Runge-Kutta实现多自由度系统响应的MATLAB程序(自振动论坛）

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)];