求助:高手帮忙看一下程序
下面是一般力学和振动理论上的一个程序用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=; % Creates force matrix
% x0=; % 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 =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 于 2007-4-2 12:03 发表
请将你的调用方式及出错信息先给出.
xjzuo说得对，我也已经在置顶贴中做了说明，看来需要规定与版规或者置顶帖不符合的帖子一律删除才行 不好意思,先谢谢xjzhou和老八,我的调用方式如下:
% Example
t=0:.5:20; % Creates time vector
u=; % Creates force matrix
x0=; % 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. 请将 vtb1_3 函数给出.
虽然我修改后已经可以画出z(1,:), z(2,:)曲线. 谢谢xjzuo,我也不知道vtb1_3是个什么子函数，原来程序就是我在1楼贴出来的那个，程序中并没有提及vtb1_3，但从调用的函数看，有可能就是调用微分方程的表达式。能把你改的程序贴出来吗？谢谢了
回复
很难想象你都不知道vtb1_3是什么...?程序我稍微修改了一下, 直接存为myrk.m, 再调用之，即可．
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myrk
t=0:.5:20; % Creates time vector
u=; % Creates force matrix
x0=; % 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)]; 谢谢xjzuo了，其实我的问题我用了ode45()算了，看这个R-K法程序是求二阶微分方程的，和我的问题相似，我就想看这个程序的效果如何，可程序就是不通，我手比较低，想看一下高手怎么调的，只好求助了，以后还会向你讨教，再次感谢xjzuo!至于vtb1_3是否是原创者把他的程序存为这个名字、抑或笔误就不得而知了，重要的是这个程序能用R-K法求解就成。
页:
[1]