四阶龙骨库塔法求解动力方程报错求高手指点

软件宝贝

求解如此方程D2y+3.14*Dy+0.002*y=25*sin(25.12*t)
function ydot=vdpol(t,y)
ydot(1)=25*sin(25.12*t)-3.14*y(1)-985.96*y(2);
ydot(2)=y(1);
clear
t0=0;
tN=20;
y0=[0:0.25];
h=0.001;
t=t0:h:tN;
N=length(t);
j=1;
for i=1:N
t1=t0+h;
K1=vdpol(t0,y0);
K2=vdpol(t0+h/2,y0+h*K1'/2);
K3=vdpol(t0+h/2,y0+h*K2'/2);
K4=vdpol(t0+h,y0+h*K3');
y1=y0+(h/6)*(K1'+2*K2'+2*K3'+K4');
yy1(j)=y1(1);
yy2(j)=y1(2);
t0=t1;
y0=y1;
j=j+1;
end
??? Attempted to access y(2); index out of bounds because numel(y)=1.

Error in ==> vdpol at 2
ydot(1)=25*sin(25.12*t)-3.14*y(1)-985.96*y(2);

博大广阔

y应该是列向量

博大广阔

ydot(1)=25*sin(25.12*t)-3.14*y(1)-985.96*y(2);
ydot(2)=y(1);
改为
ydot(1，：)=25*sin(25.12*t)-3.14*y(1)-985.96*y(2);
ydot(2，：)=y(1);

软件宝贝

回复 3 # 博大广阔 的帖子

非常感谢:@)

软件宝贝

回复 3 # 博大广阔 的帖子

博大广阔楼主你好，我这样修改以后
ydot(1，：)=25*sin(25.12*t)-3.14*y(1)-985.96*y(2);
ydot(2，：)=y(1);
还是会出现同样的问题，还望高手进一步赐教

博大广阔

function aa=aa()
clc
clear all
t0=0;
tN=20;
y0=[0;0.25];
h=0.001;
t=t0:h:tN;
N=length(t);
j=1;
for i=1:N
t1=t0+h;
K1=vdpol_1(t0,y0);
K2=vdpol_1(t0+h/2,y0+h*K1/2);
K3=vdpol_1(t0+h/2,y0+h*K2/2);
K4=vdpol_1(t0+h,y0+h*K3);
y1=y0+(h/6)*(K1+2*K2+2*K3+K4);
yy1(j)=y1(1);
yy2(j)=y1(2);
t0=t1;
y0=y1;
j=j+1;
end
end
function ydot=vdpol_1(t,y)

ydot(1,:)=25*sin(25.12*t)-3.14*y(1)-985.96*y(2);
ydot(2,:)=y(1);
end

博大广阔

新起的函数名不能已有的函数名重合、、再者主函数应该放在前面。调用函数的格式不对。我给你修正过来了。