Index in position 1 exceeds array bounds. Index must not exceed 1.

Hello everyone,

I hope you are all doing well. I am using MATLAB to solve the differential equations. But I found the problem always happens when I increased the variables.

My codes for m.file and function please see below.

Many thanks for your help and supporting.

Best wishes,

Yu

clear; clc;
tic
% Parameters
g = 9.81;
ms = 17839000;
cs = 1.1029*10^6;
ct = 6.9515*10^7;
cp = 3.4079*10^6;
ch = 1.3*10^5;
ks = 6.6026*10^4;
kp = -3.0746*10^8;
kh = 4.470*10^6;
kt = 1.5944*10^10;
mt = 2254000; % Tower+RNA
m = 20093000;
Ip = 1.251*10^10; % platform inertia moment of pitch
It = 2.9359*10^9;
ma = 9.64*10^6; % Added mass for platform surge
mh = 2.480*10^7; % Added mass for platform heave
Ia = 1.16*10^10; % Added mass for platform pitch
hs = 14.94;
ht = 56.50;
Iac = -1.01*10^8;
z = 14;
height_t = 129.13;
options = odeset('reltol',1e-13);
for x = 0:0.1:10  % Initial displacement
    input_p = x * pi /180; % Platform initial angle
    h = figure ;
    input_t = 0; % TTD initial angle
    d_t = .125 ;
    t_span = [0,0.1,100]; % Time span
    y0 = [0 0 0 0 0 input_t 0 input_p];
    % Perform integration to find displacements
    [t_out,y_out] = ode45(@(t,y) Reduced_Degree_model(t,y,ma,ht,ms,m,hs,Iac,ks,cs,mh,kh,ch,It,mt,g,kt,ct,Ip,Ia,z,cp), t_span, y0, options);
    PtfmPitch_deg = ( y_out (: ,7) *180/ pi ) ; % PtfmPitch_deg time domain output
    StDev_PtfmPitch_deg = std ( PtfmPitch_deg ) ; % PtfmPitch_deg standard deviation
    TDspFA_deg = ( y_out (: ,5) *180/ pi ) ; % TTDspFA_deg time domain output
    StDev_TTDspFA_deg = std ( TTDspFA_deg ) ; % TTDspFA_deg standard deviation
    TTDspFA_m = height_t * sind (( y_out (: ,5) *180/ pi ) ) ; % TTDspFA_m time domain output
    StDev_TTDspFA_m = std ( TTDspFA_m ) ; % TTDspFA_m standard deviation
    subplot (2 ,2 ,1) ; plot ( t_out , PtfmPitch_deg ) ; hold on ; grid on ; xlabel ('time (s)') ;
    ylabel ('Pitch (deg)') ; title ('Pitch (deg)') 
    subplot (2 ,2 ,2) ; plot ( t_out , TTDspFA_m ) ; hold on ; grid on ; xlabel ('time (s)') ;
    ylabel ('TTDsp (m)') ; title ('TTDsp (m)') ; 
    subplot (2 ,2 ,3) ; plot ( t_out , TTDspFA_deg ) ; hold on ; grid on ; xlabel ('time (s)') ;
    ylabel ('TTDsp (deg )') ; title ('TTDsp (deg )') ;
    % Create Figure 1 , Figure 2 , Figure 3 ,...
    base_file_name = sprintf (' Ini_plat_angle % .1f.fig ',x ) ;
    fullfileName = fullfile ( data_folder , base_file_name ) ;
    saveas ( gcf , fullfileName ) ; % Saves the generated figure
end

function dy = Reduced_Degree_model(y,ma,ht,ms,m,hs,Iac,ks,cs,mh,kh,ch,It,mt,g,kt,ct,Ip,Ia,z,cp,...
    t_span, y0)

dy = zeros(8,2); 
% surge:

    dy(1,1) = y(2,1);
    dy(2,1) = 1/(m+ma)*(-mt*ht*dy(6,1)+ms*hs*dy(8,1)-Iac*dy(8,1)-ks*y(1)+ks*z*y(7)-cs*y(2));
  

% heave
    dy(3,1) = y(4,1);
    dy(4,1) = -1/(m+mh)*(kh*y(3)+ch*y(4));

% tower
    dy(5,1) = y(6,1);
    dy(6,1) = 1/(It+mt*(ht)^2)*(-mt*ht*dy(2,1)-(kt+mt*g*ht)*y(5)+kt*y(7)-ct*y(6)+ct*y(8));

% platform
    dy(7,1) = y(8,1);
    dy(8,1) = 1/(Ip+ms*hs^2+Ia)*((ms*hs-Iac)*dy(2,1)+ks*z*y(1)+kt*y(5)-(kp+ks*(z)^2+kt)*y(7)+cy*y(6)+(ct+cp)*y(8));

Answer 1

[Sam Chak]

Hi @Yu

In addition to the technical points raised by @Jon and @Dyuman Joshi, there is a major issue with the way you have modeled the dynamics. If you observe the three state equations, you'll notice that they are coupled together. Specifically, depends on and , but both and also depend on , creating interdependent nested loops. Unfortunately, this is not allowed in the MATLAB ODE solver, at least not in the presented form.

dy(2,1) = 1/(m + ma)*(- mt*ht*dy(6,1) + ms*hs*dy(8,1) - Iac*dy(8,1) - ks*y(1) + ks*z*y(7) - cs*y(2));
%                             ^^^^^^^         ^^^^^^^
dy(6,1) = 1/(It + mt*(ht)^2)*(- mt*ht*dy(2,1) - (kt + mt*g*ht)*y(5) + kt*y(7) - ct*y(6) + ct*y(8));
%                                     ^^^^^^^
dy(8,1) = 1/(Ip + ms*hs^2 + Ia)*((ms*hs - Iac)*dy(2,1) + ks*z*y(1) + kt*y(5) - (kp + ks*(z)^2 + kt)*y(7) + cy*y(6) + (ct + cp)*y(8));
%                                              ^^^^^^^

To address this issue, you can rearrange the equations by moving all the first-order derivatives to the left-hand side and create a Mass matrix. This special matrix can be specified using the odeset() command. Alternatively, you can utilize the odeToVectorField() command to convert the improper ODEs into a standardized system of first-order ODEs. The latter approach involves some symbolic work and requires the Symbolic Math Toolbox.

[Yu]

Hello Sam,

Thank you for your clarifications. The motion equations are associated with a 4 DOFs system. Therefore, there are some coupling terms. Actually, it is also what I was thinking. I am not sure whether ODE sovler could be feasible for the coupling equations.

I have seen that some people solved the matrixes by using ODE sovler. Can I rearrange the equations into matrix form and use ODE solver to solve them? There are mass matrix, damping matrix, and stiffness matrix.

[Sam Chak]

Hi @Yu

I wanted to remind you that there are some examples provided in the links I shared earlier. If you haven't had a chance to do so yet, you can take the time to review them and work on the code tomorrow. Once you have reviewed the examples, please feel free to post any updated code or let me know if you have any questions. We are here to assist you.

[Yu]

Hi @Sam Chak,

Thank you for your guildance. Now it seems that the codes work after I transfer the equation to martix form and use the 'odefn'. But I am not sure wether it is correct (I run it successfully however it seems not very accurate.). I am not sure if I am supposed to set the calculation accuracy. Please see the codes below:

z = 14; 
ht = 56.50;
hp = 14.94;
height_t = 129.13;
g = 9.81;
m = 20093000; % total mass
mp = 17839000; % platform mass
mt = 2254000; % tower+RNA mass
ma = 9.64*10^6; % Added mass for platform surge
mh = 2.480*10^7; % Added mass for platform heave
Ia = 1.16*10^10; % Added mass for platform pitch
Iac = -1.01*10^8;
cs = 1.1029*10^6;
ct = 6.9515*10^7;
cp = 3.4079*10^6;
ch = 1.3*10^5;
ks = 6.6026*10^4;
kp = -3.0746*10^8;
kh = 4.470*10^6;
kt = 1.5944*10^10;
tspan = 0:200;

X0 = [0;10;10;0;0;0;0;0]; % [x10;x20;x30;x40; dx1dt0;dx2dt0;dx3dt0;dx4dt0]

[t, X] = ode45(@(t,X) odefn(t,X,tspan,g,m,ms,mt,It,Ip,ma,Ia,mh,Iac,ks,kt,kp,kh,cs,ct,cp,ch,z,ht,hp), tspan, X0);

Unrecognized function or variable 'ms'.

Error in solution>@(t,X)odefn(t,X,tspan,g,m,ms,mt,It,Ip,ma,Ia,mh,Iac,ks,kt,kp,kh,cs,ct,cp,ch,z,ht,hp) (line 25)
[t, X] = ode45(@(t,X) odefn(t,X,tspan,g,m,ms,mt,It,Ip,ma,Ia,mh,Iac,ks,kt,kp,kh,cs,ct,cp,ch,z,ht,hp), tspan, X0);

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

subplot(3,1,1)
plot(t,X(:,1)),grid
xlabel('t'),ylabel('surge')
subplot(3,1,2)
plot(t,X(:,2)),grid
xlabel('t'),ylabel('Tower Pitch')
subplot(3,1,3)
plot(t,X(:,3)),grid
xlabel('t'),ylabel('Platform Pitch')

function dXdt = odefn(t,X,tspan,g,m,mp,mt,ma,Ia,It,Ip,mh,Iac,ks,kt,kp,kh,cs,ct,cp,ch,z,ht,hp)

x = X(1:4);
xdot = X(5:8);
M = [m mt*ht Iac 0;...
    mt*ht It 0 0;...
    Iac 0 Ip 0;...
    0 0 0 m];
A = [ma 0 -mp*hp 0;...
    0 mt*ht^2 0 0;...
    -mp*hp 0 Ia 0;...
    0 0 0 mh];
C = [cs 0 0 0;...
    0 ct -ct 0;...
    0 -ct ct+cp 0;...
    0 0 0 ch];
K = [ks 0 -ks*z 0;...
    0 kt-mt*ht -kt 0;...
   -ks*z -kt kp+ks*z^2+kt+mp*g*hp 0;...
    0 0 0 kh];

xddot = (M+A)\(-K*x-C*xdot);

dXdt = [xdot; xddot];

end

Thank you for your support!

Best wishes,

Yu

[Torsten]

The list of variables for the ode function is inconsistent (see above).

For higher accuracy, use e.g.

options = odeset('RelTol',1e-10,'AbsTol',1e-10)
[t, X] = ode45(@(t,X) odefn(t,X,tspan,g,m,ms,mt,It,Ip,ma,Ia,mh,Iac,ks,kt,kp,kh,cs,ct,cp,ch,z,ht,hp), tspan, X0, options);

Answer 2

[Jon]

It looks like your argument list for your function Reduced_Degree_model, is not consistent with the call you make to ode45 with the anonymous function. I immediately see that the first argument in the call is t, as it should be, but the first argument in the function is y. So time, a scalar is getting passed in as y, and then you try to access the second element of a scalar and you get the error message you report.

You should check that the rest of the arguments for Reduced_Degree_model are consistent between the call to ode45 with the anonymous function and the function definition

[Yu]

Sorry for my late response and many thanks for your comments. I am going to check my codes with your suggestions.

Answer 3

[Dyuman Joshi]

• You have not defined "t" and many other variables as input to the function "Reduced_Degree_model".
• You have defined "tspan" and "y0" as inputs to "Reduced_Degree_model", which is incorrect.
• The variable "cy" in the last line of code (where dy(8,1) is defined) seems to be a typo.
• dy is to be pre-allocated as 8x1 not 8x2.
• There are typos in names of variables used in the for loop.

There might be other mistakes in your code. I suggest you go through your code line by line and check thoroughly.

[Yu]

Thank you for your detailed suggestions. I will modify my codes following your suggestions.