Array indices for differential equations

Hello, I have a function file (rk4.m) that contains equations that desine a system of motion as shown below.

function xv=rk4(x,t,dt,Input)
  
   dth = dt/2;
   t2 = t + dth;
   
   dummy_x = ax(x,t,Input);
   xA = x + dth*dummy_x;
   
%     Input.thm1 = Input.thm1_h;
%     Input.dthm1 = Input.dthm1_h;
   
   dummy_x = ax(xA,t2,Input);
   xB = dth*dummy_x;
   xB2 = x + 2*xB;
   xB = x + xB;
   
   dummy_x = ax(xB,t2,Input);
   xC = x + dt*dummy_x;
   
   dummy_x = ax(xC,t2,Input);
   xD = xC + dth*dummy_x;
   
   xv = (xA + xB2 + xD)/3;
   
end

function dx = ax(x,t,In)
global J1 J2 Jn c1 c2 k1 k2 cn kn kn3 T beta vin times
%
% Below are the differential equations to be solved, written in 1st order
% format
%
UJvibs=x(6,1)^2*(cos(beta)*(sin(beta))^2*sin(2*x(5,1)))/(1-(sin(beta))^2*(cos(x(5,1)))^2)^2;
%
dx(1,1)=x(2,1); % shaft 1
dx(2,1)=(-c1*(x(2,1)-x(8,1))-k1*(x(1,1)-x(7,1))+cn*(x(4,1)-x(2,1))+kn*(x(3,1)-x(1,1))+kn3*(x(3,1)-x(1,1))^5+c2*(x(10,1)-x(2,1))+k2*(x(7,1)-x(1,1)))/(J1);%
dx(3,1)=x(4,1); % NES
dx(4,1)=(-cn*(x(4,1)-x(2,1))-kn*(x(3,1)-x(1,1))-kn3*(x(3,1)-x(1,1))^5)/Jn;
dx(5,1)=x(6,1); % UJ input
x(6,1) = vin+0.02*vin*cos(36*vin*times(i))+0.005*vin*cos(72*vin*times(i));
dx(7,1)=x(8,1); % UJ output
dx(8,1)=x(6,1)*cos(beta)/(1-(sin(beta))^2*(cos(x(5,1)))^2) - UJvibs;
dx(9,1)=x(10,1); % shaft 2
dx(10,1)=(-c2*(x(10,1)-x(2,1))-k2*(x(9,1)-x(1,1)))/J2;

When i run the script file, I get the following error:

Array indices must be positive integers or logical values.

Error in rk4>ax (line 40)

x(6,1) = vin+0.02*vin*cos(36*vin*times(i))+0.005*vin*cos(72*vin*times(i));

I can see that the array sizes are different due to the time values introduced but I don't know how to fix this. The script file is below for reference. There is addition function file (Main.m) however this does not effect this issue.

disp('NES status?');
disp('type L for Locked or A for Active');
Clf=input('','s');
comp1=strcmp(Clf,'L');
comp2=strcmp(Clf,'l');
if comp1==1
    NESstate='l';
    NEScase='locked';
disp('Locked case running'); 
elseif comp2==1
    NESstate='l';
    NEScase='locked';
disp('Locked case running'); 
else
    NESstate='a';
    NEScase='active';
disp('Active case running');   
end
%% Make parameters available to other functions
%
global times dt Fs J1 J2 Jn J_UJ J_UJ_out J_UJ_in c1 c2 k1 k2 cn kn kn3 beta T cM kM UJacc vin

%
%% Definition of parameters and preparation for solving the system of ODEs
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
dt=2e-4; % taken from sampling rate used in experiments 2e-4
Fs=1/dt;

% %
 times=0:dt:5;
 times=times(1):dt:times(end);

%
 fstart=0/60;      % initial shaft speed (Hz)
 fend=100/60;       % end shaft speed
%  Input.w0 = 2*pi*fstart; % convert initial shaft speed to rad/s
%  Input.wrate = 2*pi*(fend-fstart)/(times(end)-times(1)); % calculate the rate with which the shaft speed increases
%
    beta = 20*pi/180; %confirmed was 19.8 
    J_UJ_in=0.06543;%%confirmed
    J_UJ_out=9.0737e-5; %confirmed
    J_UJ=J_UJ_in+J_UJ_out;
%     T = Input.wrate; % rate of mean speed
    vin = 60;
    c1=0.15; % damping for shaft was 2
    k1=540; %confirmed
    J2=0.002931; % PTO inertia
    c2=0.2; % PTO damping was 2
    k2=8500; % PTO stiffness (confirmed: shaft + coupling in series)
    Jn=4.0085e-4; %+1.25e-4(RING Inertia for test vehicle) confirmed WAS 4.0085E-4
    UJacc=T;
    Jmock=3.31e-4; % inertia of replacement disk for locked tests


if NESstate=='l'
    J1=1.13e-4+Jmock;% note shaft inertia without the impactor;
    cn=0;
    kn=0;
    kn3=0;
    M = [J1,0;0,J2];
    K = [k1+k2,-k2;-k2,k2];
    [V,D]=eig(M\K);
    w=D.^0.5/2/pi;

    
elseif NESstate=='a'
    J1=2.08e-4; % shaft inertia with impactor
    cn=0.002; % NES damping
    kn=0.9401; % NES linear confirmed
    kn3=21.185; % NES nonlinear quintic confirmed was 38870
    M = [J1,0,0;0,Jn,0;0,0,J2];
    K = [k1+kn+k2,-kn,-k2;-kn,kn,0;-k2,0,k2];
    [V,D]=eig(M\K);
    w=D.^0.5/2/pi;
end
    %
tic % this command initialises a time counter. The program will count how 
%     much time is spent until it reaches the command "toc" 
%
%% Call function Main to perform the calculations
%
x= Main([]);

Thak you for the help in advance!

[Torsten]

A variable "i" is nowhere set in the function "ax". Thus you cannot reference it.

Further it seems strange that you set

x(6,1) = vin+0.02*vin*cos(36*vin*times(i))+0.005*vin*cos(72*vin*times(i));

and not

dx(6,1) = vin+0.02*vin*cos(36*vin*times(i))+0.005*vin*cos(72*vin*times(i));

[Sam Butler]

By setting the equation equal to dx(6,1) that would give the derivative of x(6,1). The equation shown is for x(6,1), the derivative of x(6,1) is not required.

Could you provide an example in my code of how the i variable can be used? Sorry I am new to MATLAB.

[Torsten]

By setting the equation equal to dx(6,1) that would give the derivative of x(6,1). The equation shown is for x(6,1), the derivative of x(6,1) is not required.

If you know

x(6,1) = vin+0.02*vin*cos(36*vin*times(i))+0.005*vin*cos(72*vin*times(i));

you must remove it from the variables you solve for. You should set it at the beginning of function "ax" as

x6 = vin+0.02*vin*cos(36*vin*t)+0.005*vin*cos(72*vin*t);

This also solves the problem with the unknown "i".

[Sam Butler]

Thank you for your response!

Like so?

function xv=rk4(x,t,dt,Input)
  
   dth = dt/2;
   t2 = t + dth;
   
   dummy_x = ax(x,t,Input);
   xA = x + dth*dummy_x;
   
%     Input.thm1 = Input.thm1_h;
%     Input.dthm1 = Input.dthm1_h;
   
   dummy_x = ax(xA,t2,Input);
   xB = dth*dummy_x;
   xB2 = x + 2*xB;
   xB = x + xB;
   
   dummy_x = ax(xB,t2,Input);
   xC = x + dt*dummy_x;
   
   dummy_x = ax(xC,t2,Input);
   xD = xC + dth*dummy_x;
   
   xv = (xA + xB2 + xD)/3;
   
end

function dx = ax(x,t,In)


global J1 J2 Jn c1 c2 k1 k2 cn kn kn3 T beta vin times

x6 = vin+0.02*vin*cos(36*vin*t)+0.005*vin*cos(72*vin*t);
%
UJvibs=x(6,1)^2*(cos(beta)*(sin(beta))^2*sin(2*x(5,1)))/(1-(sin(beta))^2*(cos(x(5,1)))^2)^2;
%
dx(1,1)=x(2,1); % shaft 1
dx(2,1)=(-c1*(x(2,1)-x(8,1))-k1*(x(1,1)-x(7,1))+cn*(x(4,1)-x(2,1))+kn*(x(3,1)-x(1,1))+kn3*(x(3,1)-x(1,1))^5+c2*(x(10,1)-x(2,1))+k2*(x(7,1)-x(1,1)))/(J1);%
dx(3,1)=x(4,1); % NES
dx(4,1)=(-cn*(x(4,1)-x(2,1))-kn*(x(3,1)-x(1,1))-kn3*(x(3,1)-x(1,1))^5)/Jn;
dx(5,1)=x6; % UJ input
% x(6,1) = vin+0.02*vin*cos(36*vin*times(i))+0.005*vin*cos(72*vin*times(i));
dx(7,1)=x(8,1); % UJ output
dx(8,1)=x(6,1)*cos(beta)/(1-(sin(beta))^2*(cos(x(5,1)))^2) - UJvibs;
dx(9,1)=x(10,1); % shaft 2
dx(10,1)=(-c2*(x(10,1)-x(2,1))-k2*(x(9,1)-x(1,1)))/J2;
end

Also, as there was no time domain in the equations before, and just constants, the graph produced straight lines. Now this x6 is the input velocity and is also time dependant, do the other dx equations need to be transferred to be in the time domain?

I'm new to MATLAB so appreciate all the help I can get.

[Torsten]

Now this x6 is the input velocity and is also time dependant, do the other dx equations need to be transferred to be in the time domain?

I don't know. I only see that you lack an unknown now (x(6)) and the dy are not counted continuously from 1 to 9, but from 1 to 10 with a dx(6,1) missing. Further, you still use x(6,1) in the computation of dx(8,1).

You should check your equations for validity.

[Sam Butler]

The input velocity is time dependent yes, so I believe all the other equations will need to be time dependent. The x6 is the velocity from a motor which 'drives' the system. What is the best way to change all of the equations to become time depndent?

It makes more sense for me to change it from 1 to 9 also i believe.

[Torsten]

If one of the inputs to the equations is time-dependent, the other equations will be time-dependent as well.

So if your original equations were correct, they will remain correct.