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

10 次查看（过去 30 天）
Yu
Yu 2024-1-29
编辑： Torsten 2024-3-4
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
Sam Chak 2024-1-29
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.
  4 个评论
显示 2更早的评论
Yu
Yu 2024-3-4
Ran in:
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
Torsten 2024-3-4
编辑：Torsten 2024-3-4
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);

请先登录，再进行评论。

更多回答（2 个）

Jon
Jon 2024-1-29
编辑：Jon 2024-1-29
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
  1 个评论
Yu
Yu 2024-1-30
Sorry for my late response and many thanks for your comments. I am going to check my codes with your suggestions.

请先登录，再进行评论。

Dyuman Joshi
Dyuman Joshi 2024-1-29
  • 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.

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by