请教ode45步长问题
我有一个方程组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=;%初始值
tspan=linspace(0,53.74,2688)';%0~53.74取2688个点刚好是步长为0.02
=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=;
FII=;
dy=A*y+FII;%y=
好像oed45是变步长的，我用ode45是不是不行？若行我该怎么处理啊？不行该用什么方法呢？
我自己编写的ode45也不行。
y0=;%初始值
=rk4(@funX,0,53.74,2688,y0);
function =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 编辑 ]
回复 楼主 julian 的帖子
步长应该是可以自己设定的，如t=linspace(0,30,50)doc ode45
回复 沙发 ch_j1985 的帖子
首先谢谢你的建议，这个我知道，但我在设定以后，我在被调用函数内要使用积分时间点，但为什么被调用函数内不是按步长增加的呢？即：
tspan=linspace(0,53.74,2688)';%0~53.74取2688个点刚好是步长为0.02，即时间增长为0，0.02，0.04，0.08......
=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=;
FII=;
dy=A*y+FII;%y=
页:
[1]