# Array indices for differential equations

1 visualizzazione (ultimi 30 giorni)
Sam Butler il 19 Mar 2023
Commentato: Torsten il 20 Mar 2023
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 CommentiMostra 5 commenti meno recentiNascondi 5 commenti meno recenti
Sam Butler il 20 Mar 2023
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 il 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.

Accedi per commentare.

### Categorie

Scopri di più su Programming in Help Center e File Exchange

R2022b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by