请教ode45步长问题

julian

我有一个方程组B*ddq+C*q=-II*ddx(:,2)  
变化为：
dy=A*y+FII,
FII=-II*ddx(N,2);
其中ddq是q的二阶导数，dy是y的一阶导数，A为8*8矩阵，y,f为八阶向量，ddx为地震加速度,2688行，2列，第一列为时间，第二列为加速度值。

已知ddx在0~53.74秒了取2688个值,时间步长为0.02.我用ode45时也取时间步长为0.02，在每一个加速度值时用oed45积分一次。我下面的处理好像不对，大侠帮忙啊！

我是这么处理的，但每次都说下面的N不是整数？我该怎么办？
y0=[0;0;0;0;0;0;0;0];%初始值
tspan=linspace(0,53.74,2688)';%0~53.74取2688个点刚好是步长为0.02
[t,y]=ode45(@funX,tspan,y0);

function dy=funX(t,y)
global B ddx C II;

N=t/0.02;%%%%%%%%%%关键就此处,t不能依次0，0.02,0.04,0.06....53.74

if N==0
N=1;
else
N=N+1;
end

FI=-II*ddx(N,2);
ZERO=zeros(4);I=eye(4);
A=[ZERO,I;-B\C,ZERO];
FII=[0;0;0;0;B\FI];
dy=A*y+FII;%y=[q1,q2,q3,q4,dq1,dq2,dq3,dq4]

好像oed45是变步长的，我用ode45是不是不行？若行我该怎么处理啊？不行该用什么方法呢？

我自己编写的ode45也不行。

y0=[0;0;0;0;0;0;0;0];%初始值
[t,y]=rk4(@funX,0,53.74,2688,y0);

function [t,y]=rk4(f,aa,bb,M,y0)
N=length(y0);
h=(bb-aa)/M;%步长
y=zeros(N,M+1);
T=zeros(1,M+1);
t=zeros(1,M+1);
Y=zeros(N,M+1);
T=aa:h:bb;
t(1)=aa;t(M+1)=bb;
Y(:,1)=y0;
y(:,1)=Y(:,1);

for j=1:M
    K1=h*feval(f,T(j),Y(:,j));
    K2=h*feval(f,T(j)+h/2,Y(:,j)+K1/2);
    K3=h*feval(f,T(j)+h/2,Y(:,j)+K2/2);
    K4=h*feval(f,T(j)+h,Y(:,j)+K3);
    Y(:,j+1)=Y(:,j)+(K1+2*K2+2*K3+K4)/6;
    y(:,j+1)=Y(:,j+1);
    t(:,j+1)=T(:,j+1);
end

t=t';
y=y';

大侠帮帮忙哈！

[ 本帖最后由 julian 于 2008-11-1 15:59 编辑 ]

ch_j1985

步长应该是可以自己设定的，如t=linspace(0,30,50)
doc ode45

julian

首先谢谢你的建议，这个我知道，但我在设定以后，我在被调用函数内要使用积分时间点，但为什么被调用函数内不是按步长增加的呢？
即：
tspan=linspace(0,53.74,2688)';%0~53.74取2688个点刚好是步长为0.02，即时间增长为0，0.02，0.04，0.08......
[t,y]=ode45(@funX,tspan,y0);

%%%调用函数
function dy=funX(t,y)
global B ddx C II;

N=t/0.02;%%%%%%%%%%关键就此处,t不能依次0，0.02,0.04,0.06....53.74

if N==0
N=1;
else
N=N+1;
end

FI=-II*ddx(N,2);
ZERO=zeros(4);I=eye(4);
A=[ZERO,I;-B\C,ZERO];
FII=[0;0;0;0;B\FI];
dy=A*y+FII;%y=[q1,q2,q3,q4,dq1,dq2,dq3,dq4]