求助:高手帮忙看一下程序

handb

下面是一般力学和振动理论上的一个程序
用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)];
我按照例子中的exemple去试程序,总是出错,麻烦高人帮忙看看是什么原因?谢谢

xjzuo

请将你的调用方式及出错信息先给出.

eight

原帖由 xjzuo 于 2007-4-2 12:03 发表
请将你的调用方式及出错信息先给出.

xjzuo说得对，我也已经在置顶贴中做了说明，看来需要规定与版规或者置顶帖不符合的帖子一律删除才行

handb

不好意思,先谢谢xjzhou和老八,我的调用方式如下:
%       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.
这个程序怎么总出错?提示如下
??? Error: File: E:\MATLAB6p5p1\work\vbr_rk.m Line: 23 Column: 16
Missing variable or function.

xjzuo

请将 vtb1_3 函数给出.

虽然我修改后已经可以画出z(1,:), z(2,:)曲线.

handb

谢谢xjzuo,我也不知道vtb1_3是个什么子函数，原来程序就是我在1楼贴出来的那个，程序中并没有提及vtb1_3，但从调用的函数看，有可能就是调用微分方程的表达式。能把你改的程序贴出来吗？谢谢了

xjzuo

很难想象你都不知道vtb1_3是什么...?
程序我稍微修改了一下, 直接存为myrk.m, 再调用之，即可．
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myrk
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.
z=vbrrk('vbrrkex',u,t,x0);      
plot(t,z(1,:));        % Plots displacement versus time.
hold on
plot(t,z(2,:),'r-');   % Plots velocity versus time.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      
function z = vbrrk(rkf,u,t,x0)

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function zd=vbrrkex(z,u,t)
zd=[z(2);
    -.1*z(2)-z(1)+u(2)];

handb

谢谢xjzuo了，其实我的问题我用了ode45()算了，看这个R-K法程序是求二阶微分方程的，和我的问题相似，我就想看这个程序的效果如何，可程序就是不通，我手比较低，想看一下高手怎么调的，只好求助了，以后还会向你讨教，再次感谢xjzuo!至于vtb1_3是否是原创者把他的程序存为这个名字、抑或笔误就不得而知了，重要的是这个程序能用R-K法求解就成。