ode45 error-time-dependent parameter

Dear all,

I would appreciate it if you could help me fix the error. 9 time-dependent parameter are existed (A_x, ..., D_t).

The error is:

Index exceeds matrix dimensions.

Error in STdof/MDOF3st3dof (line 56)

dg13=g(14);

Error in odearguments (line 87)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to

yp0.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan,

y0, options, varargin);

Error in STdof (line 34)

[t,g] =

ode45(@MDOF3st3dof,tspan,g0,options,tspan,A_x,V_x,D_x,A_z,V_z,D_z,A_t,V_t,D_t,m1x,m2x,m3x,m1z,m2z,m3z,I1,I2,I3,c1x,c2x,c3x,c1z,c2z,c3z,c1t,c2t,c3t,k1x,k2x,k3x,k1z,k2z,k3z,k1t,k2t,k3t);

function STdof
close all, clear,  clc
 
t_start=0;
t_end=32.82;
dt=0.02;
 
A_x=load ('D:\tax.out');  
V_x=load ('D:\tvx.out');  
D_x=load ('D:\tdx.out');  
A_z=load ('D:\taz.out');  
V_z=load ('D:\tvz.out');  
D_z=load ('D:\tdz.out');  
A_t=load ('D:\tat.out');  
V_t=load ('D:\tvt.out');  
D_t=load ('D:\tdt.out');  

L1=3;   L2=6; 
m1x=100;    m2x=150;    m3x=200;
m1z=400;    m2z=300;    m3z=210;
I1=130;     I2=160;     I3=190;
c1x=22000;  c2x=18000;  c3x=15000;
c1z=13000;  c2z=15000;  c3z=14000;
c1t=23000;  c2t=18000;  c3t=16000;
k1x=5000000;   k2x=4000000;    k3x=7600000;
k1z=6000000;   k2z=4400000;    k3z=7400000;
k1t=32000;     k2t=42000;      k3t=6500000;
%%========================================================== 
g0=[0,0,0,0,0,0,0,0,0,0,0,0];
Nstep=(t_end-t_start)/dt+1;
tspan=t_start:dt:t_end;
options = odeset('RelTol',1e-9,'AbsTol',1e-9);    
 
[t,g] = ode45(@MDOF3st3dof,tspan,g0,options,tspan,A_x,V_x,D_x,A_z,V_z,D_z,A_t,V_t,D_t,m1x,m2x,m3x,m1z,m2z,m3z,I1,I2,I3,c1x,c2x,c3x,c1z,c2z,c3z,c1t,c2t,c3t,k1x,k2x,k3x,k1z,k2z,k3z,k1t,k2t,k3t);
 
 function dg=MDOF3st3dof(t,g,t1,A_x,V_x,D_x,A_z,V_z,D_z,A_t,V_t,D_t,m1x,m2x,m3x,m1z,m2z,m3z,I1,I2,I3,c1x,c2x,c3x,c1z,c2z,c3z,c1t,c2t,c3t,k1x,k2x,k3x,k1z,k2z,k3z,k1t,k2t,k3t)
   
 A_x = interp1(t1,A_x,t)      
 V_x = interp1(t1,V_x,t) 
 D_x = interp1(t1,D_x,t) 
 A_z = interp1(t1,A_z,t) 
 V_z = interp1(t1,V_z,t) 
 D_z = interp1(t1,D_z,t) 
 A_t = interp1(t1,A_t,t) 
 V_t = interp1(t1,V_t,t) 
 D_t = interp1(t1,D_t,t) 
 
 dg1=g(2);
 dg2=(-1/m1x)*(c1x*g(2)-c1x*V_x+k1x*g(1)-k1x*D_x);
 dg5=g(6);
 dg6=(-1/m3x)*(m2x*A_x-c1x*g(2)+c1x*V_x-k1x*g(1)+k1x*D_x+c3x*g(6)+k3x*g(5));
 dg7=g(8);
 dg8=(-1/m1z)*(2*c1z*g(8)-2*c1z*V_z+2*k1z*g(7)-2*k1z*D_z);
 dg11=g(12);
 dg12=(-1/m3z)*(m2z*A_z-2*c1z*V_z-2*k1z*g(7)+2*k1z*D_z+c3z*g(12)+k3z*g(11));
 dg13=g(14);
 dg14=(-L1^2/(2*I1))*(c1z*g(14)-c1z*V_t+k1z*g(13)-k1z*D_t);
 dg17=g(18);
 dg18=(-1/I3)*(I2*A_t-L1^2/2*c1z*g(14)+L1^2/2*c1z*V_t-L1^2/2*k1z*g(13)+L1^2/2*k1z*D_t+c3t*g(18)+k3t*g(17));
 
dg=[dg1;dg2;dg5;dg6;dg7;dg8;dg11;dg12;dg13;dg14;dg17;dg18];
 
end
save('result.txt','g','-ASCII')
display 'Done!'
end

Answer 1

[Sam Chak]

Hi @Elahe S

You may not fully understand some of the MATLAB error messages if you are not a programmer or if you haven't encountered differential equations in matrix form. The code has been fixed now, and you can learn by example. However, you need to inject, feed, or load the time-dependent parameters into the Workspace.

%% §1: Time-dependent parameters
t_start = 0;
t_end   = 6;    % 32.82;
dt      = 0.02;

tAVD    = t_start:dt:t_end;
tx      = tAVD;     % should get time data from tax, tvx, tdx
tz      = tAVD;     % should get time data from taz, tvz, tdz
tt      = tAVD;     % should get time data from tat, tvt, tdt
 
A_x     = tanh(tx); % load('D:\tax.out');  
V_x     = tanh(tx); % load('D:\tvx.out');  
D_x     = tanh(tx); % load('D:\tdx.out');  

A_z     = tanh(tz); % load('D:\taz.out');  
V_z     = tanh(tz); % load('D:\tvz.out');  
D_z     = tanh(tz); % load('D:\tdz.out');  

A_t     = tanh(tt); % load('D:\tat.out');  
V_t     = tanh(tt); % load('D:\tvt.out');  
D_t     = tanh(tt); % load('D:\tdt.out');  

%% §2: Call ode solver
tspan   = t_start:dt:t_end;
% Nstep   = (t_end-t_start)/dt + 1;
g0      = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
options = odeset('RelTol', 1e-9, 'AbsTol', 1e-9);    
[t, g]  = ode45(@(t, g) MDOF3st3dof(t, g, tAVD, A_x, V_x, D_x, A_z, V_z, D_z, A_t, V_t, D_t), tspan, g0, options);

%% §3: Plot results
tiledlayout(4, 3, 'TileSpacing', 'Compact');
nexttile, plot(t, g(:,1)),  grid on, title('g1')
nexttile, plot(t, g(:,2)),  grid on, title('g2')
nexttile, plot(t, g(:,3)),  grid on, title('g3')
nexttile, plot(t, g(:,4)),  grid on, title('g4')
nexttile, plot(t, g(:,5)),  grid on, title('g5')
nexttile, plot(t, g(:,6)),  grid on, title('g6')
nexttile, plot(t, g(:,7)),  grid on, title('g7')
nexttile, plot(t, g(:,8)),  grid on, title('g8')
nexttile, plot(t, g(:,9)),  grid on, title('g9')
nexttile, plot(t, g(:,10)), grid on, title('g10')
nexttile, plot(t, g(:,11)), grid on, title('g11')
nexttile, plot(t, g(:,12)), grid on, title('g12')

%% §4: ODE function
function dg = MDOF3st3dof(t, g, t1, A_x, V_x, D_x, A_z, V_z, D_z, A_t, V_t, D_t)
    
    %% §4.1: Interpolate time-dependent parameters
    A_x = interp1(t1, A_x, t);
    V_x = interp1(t1, V_x, t);
    D_x = interp1(t1, D_x, t);
    A_z = interp1(t1, A_z, t);
    V_z = interp1(t1, V_z, t);
    D_z = interp1(t1, D_z, t);
    A_t = interp1(t1, A_t, t);
    V_t = interp1(t1, V_t, t);
    D_t = interp1(t1, D_t, t);
    
    %% §4.2: Constants
    L1  = 3;        L2  = 6; 
    m1x = 100;      m2x = 150;      m3x = 200;
    m1z = 400;      m2z = 300;      m3z = 210;
    I1  = 130;      I2  = 160;      I3  = 190;
    c1x = 22000;    c2x = 18000;    c3x = 15000;
    c1z = 13000;    c2z = 15000;    c3z = 14000;
    c1t = 23000;    c2t = 18000;    c3t = 16000;
    k1x = 5000000;  k2x = 4000000;  k3x = 7600000;
    k1z = 6000000;  k2z = 4400000;  k3z = 7400000;
    k1t = 32000;    k2t = 42000;    k3t = 6500000;
    
    %% §4.3: Differential equations
    dg1 = g(2);
    dg2 = (-1/m1x)*(c1x*g(2) - c1x*V_x + k1x*g(1) - k1x*D_x);
    dg3 = g(4);
    dg4 = (-1/m3x)*(m2x*A_x - c1x*g(2) + c1x*V_x - k1x*g(1) + k1x*D_x + c3x*g(4) + k3x*g(3));
    dg5 = g(6);
    dg6 = (-1/m1z)*(2*c1z*g(5) - 2*c1z*V_z + 2*k1z*g(5) - 2*k1z*D_z);
    dg7 = g(8);
    dg8 = (-1/m3z)*(m2z*A_z-2*c1z*V_z - 2*k1z*g(5) + 2*k1z*D_z + c3z*g(8) + k3z*g(7));
    dg9 = g(10);
    dg10= (-L1^2/(2*I1))*(c1z*g(10) - c1z*V_t + k1z*g(9) - k1z*D_t);
    dg11= g(12);
    dg12= (-1/I3)*(I2*A_t - L1^2/2*c1z*g(10) + L1^2/2*c1z*V_t - L1^2/2*k1z*g(9) + L1^2/2*k1z*D_t + c3t*g(12) + k3t*g(11));
    
    dg  = [dg1; dg2; dg3; dg4; dg5; dg6; dg7; dg8; dg9; dg10; dg11; dg12];

end
% save('result.txt','g','-ASCII')
% display 'Done!'

[Elahe S]

thank you so much for your constructive response. I have another question.

How dg1 to dg12 could be saved and exported at each step?

using

A3x(:,1)=dg4;

save('Ac2x.txt','A3x','-ASCII')

saved A3x just at the last step.

[Torsten]

[t, g]  = ode45(@(t, g) MDOF3st3dof(t, g, tAVD, A_x, V_x, D_x, A_z, V_z, D_z, A_t, V_t, D_t), tspan, g0, options);
DG_SAVE = zeros(numel(t),12);
for i = 1:numel(t)
  dg = MDOF3st3dof(t(i), g(i,:), tAVD, A_x, V_x, D_x, A_z, V_z, D_z, A_t, V_t, D_t);
  DG_SAVE(i,:) = dg(1:12);
end

[Elahe S]

Thank you.

Answer 2

[Rishi]

Hi Elahe,

I understand that you want to know why you are getting ‘Index exceeds matrix dimensions’ as an error.

This error message indicates that you are trying to access an element of a matrix using an index that does not exist. In your code, the error occurs at the following line:

d13 = g(14);

The function 'MDOF3st3dof' is called by the 'ode45' solver, which expects the function to return the derivatives of the state vector 'g'. However, the state vector 'g' is initialized with 12 elements (g0=[0,0,0,0,0,0,0,0,0,0,0,0];), and the code is trying to access the 14th element, which doesn't exist. This is the source of the error.

To resolve the error, you need to ensure that you are only accessing elements within the bounds of the g vector. Since g has 12 elements, you should only access indices 1 through 12. The differential equations need to be corrected to match the actual dimension of the system.

Hope this helps!