Loop not returning all values

joe brady

2022 11 25

0 个回答

6 次查看（30 天）

0 个投票

Below is the code for my matlab file that is only returning u(1) and u(2):

%% Clear MATLAB
clear;
clc;
%% Parameters
tlength=16.74; % time (s)
h=0.02; % time step (s)
N=1:ceil(tlength/h);
t=0.02:h:tlength; % time steps
w0=31.416; % Natural circular freq
zeta0=0.05;
eta=2;
tau=0.12;
m=7; 
R=eta*m*w0^2;
v=zeros(1,837);
y3=zeros(1,837);
%% Initial Conditions
v(1)=1; % Velocity at t=0.02s
y3(1)=1; % Displacement of EDD at t=0s
ag(1)=-0.1; % Ground Acceleration at t=0s
a(1)=-1; % Acceleration of frame
u(1)=1; % Displacement of frame
%% Define Equations
f1=@(t,u,v) v;
f2=@(t,v,y3) v-y3/tau;
%% RK4 Loop
for i=1:N
    t(i+1)=t(i)+h;
    k1=h*f1(t(i),u(i),v(i));
    l1=h*f2(t(i),v(i),y3(i));
    k2=h*f1(t(i)+h/2,u(i)+k1*h/2,v(i)+l1*h/2);
    l2=h*f1(t(i)+h/2,v(i)+k1*h/2,y3(i)+l1*h/2);
    k3=h*f1(t(i)+h/2,u(i)+k2*h/2,v(i)+l2*h/2);
    l3=h*f1(t(i)+h/2,v(i)+k2*h/2,y3(i)+l2*h/2);
    k4=h*f1(t(i)+h,u(i)+k3*h,v(i)+l3*h);
    l4=h*f1(t(i)+h,v(i)+k3*h,y3(i)+l3*h);
    v(i+1)=v(i)+h/6*(k1+2*k2+2*k3+k4);
    y3(i+1)=y3(i)+h/6*(l1+2*l2+2*l3+l4);
    
end
plot(t,v)

5 个评论
显示 3更早的评论隐藏 3更早的评论

joe brady 2022-11-25

编辑：Jan 2022-11-26

在 MATLAB Online 中打开

Thanks for your reply, i have now added a definition for u and i am still struggling - the same error is happening:

%% Clear MATLAB
clear;
clc;
%% Parameters
tlength=16.74; % time (s)
h=0.02; % time step (s)
N=1:ceil(tlength/h);
t=0.02:h:tlength; % time steps
w0=31.416; % Natural circular freq
zeta0=0.05;
eta=2;
tau=0.12;
m=7; 
R=eta*m*w0^2;
%% Initial Conditions
t(1)=0.02;
v(1)=0; % Velocity at t=0.02s
y3(1)=0; % Displacement of EDD at t=0s
ag(1)=-0.1; % Ground Acceleration at t=0s
a(1)=0; % Acceleration of frame
u(1)=1; % Displacement of frame
vedd(1)=0; % Velocity of EDD.
%% Define Equations
f1=@(t,u,v,y3) v;
f2=@(t,u,v,y3) -u*w0^2-2*zeta0*w0*v-eta*y3*w0^2-ag;
f3=@(t,u,v,y3) v-y3/tau;
%% RK4 Loop
for i=1:N
    t(i+1)=t(i)+h;
    
    F1_1=h*f1(t(i),u(i),v(i),y3(i));
    F2_1=h*f2(t(i),u(i),v(i),y3(i));
    F3_1=h*f3(t(i),u(i),v(i),y3(i));
    F1_2=h*f1(t(i)+h/2,u(i)+h/2*F1_1,v(i)+h/2*F1_1,y3(i)+h/2*F1_1);
    F2_2=h*f2(t(i)+h/2,u(i)+h/2*F2_1,v(i)+h/2*F2_1,y3(i)+h/2*F2_1);
    F3_2=h*f3(t(i)+h/2,u(i)+h/2*F3_1,v(i)+h/2*F3_1,y3(i)+h/2*F3_1);
    F1_3=h*f1(t(i)+h/2,u(i)+h/2*F1_2,v(i)+h/2*F1_2,y3(i)+h/2*F1_2);
    F2_3=h*f2(t(i)+h/2,u(i)+h/2*F2_2,v(i)+h/2*F2_2,y3(i)+h/2*F2_2);
    F3_3=h*f3(t(i)+h/2,u(i)+h/2*F3_2,v(i)+h/2*F3_2,y3(i)+h/2*F3_2);
    F1_4=h*f1(t(i)+h,u(i)+h*F1_3,v(i)+h*F1_3,y3(i)+h*F1_3);
    F2_4=h*f2(t(i)+h,u(i)+h*F2_3,v(i)+h*F2_3,y3(i)+h*F2_3);
    F3_4=h*f3(t(i)+h,u(i)+h*F3_3,v(i)+h*F3_3,y3(i)+h*F3_3);
    v(i+1)= v(i)+1/6*(F1_1+2*F1_2+2*F1_3+F1_4);
    a(i+1)= a(i)+1/6*(F2_1+2*F2_2+2*F2_3+F2_4);
    vedd(i+1)= vedd(i)+1/6*(F3_1+2*F3_2+2*F3_3+F3_4);
end