Array indices for differential equations

조회 수: 2 (최근 30일)
Sam Butler
Sam Butler 2023년 3월 19일
댓글: Torsten 2023년 3월 20일
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
Sam Butler
Sam Butler 2023년 3월 20일
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
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개)

카테고리

Help CenterFile Exchange에서 Programming에 대해 자세히 알아보기

제품


릴리스

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by