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

15 views (last 30 days)
Yu on 29 Jan 2024
Edited: Torsten on 4 Mar 2024
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));

Sam Chak on 29 Jan 2024
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 on 4 Mar 2024
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
Best wishes,
Yu
Torsten on 4 Mar 2024
Edited: Torsten on 4 Mar 2024
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);

Jon on 29 Jan 2024
Edited: Jon on 29 Jan 2024
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 on 30 Jan 2024
Sorry for my late response and many thanks for your comments. I am going to check my codes with your suggestions.

Dyuman Joshi on 29 Jan 2024
• 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 on 30 Jan 2024
Thank you for your detailed suggestions. I will modify my codes following your suggestions.