Array indices for differential equations
Show older comments
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!
7 Comments
Torsten
on 20 Mar 2023
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
on 20 Mar 2023
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
on 20 Mar 2023
Torsten
on 20 Mar 2023
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
on 20 Mar 2023
Torsten
on 20 Mar 2023
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.
Answers (0)
Categories
Find more on Mathematics and Optimization in Help Center and File Exchange
on 19 Mar 2023
on 20 Mar 2023
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
