Array indices for differential equations
조회 수: 1 (최근 30일)
이전 댓글 표시
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
Torsten
2023년 3월 20일
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.
답변 (0개)
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!