ODE45 Multiple Degrees of Freedom Problem

Hello,

I am currently trying to model the displacement of a 3DOF system using matlabs ODE45 solver. However as much I try to fix the errors I am getting I do not seem to be making any progress. At the moment I am just trying to get the solver to run without error and will then plot this at a later date. Any advice / help would be greatly appreciated.

These are the errors that I am getting.

Thanks in advance.

[Stephen23]

Original question by Maximilian retrieved from Google Cache:

https://web.archive.org/web/20221110062251/https://webcache.googleusercontent.com/search?q=cache:RCWUlJ79uDcJ:https://www.mathworks.com/matlabcentral/answers/1847148-ode45-multiple-degrees-of-freedom-problem

"ODE45 Multiple Degrees of Freedom Problem"

Hello,

I am currently trying to model the displacement of a 3DOF system using matlabs ODE45 solver. However as much I try to fix the errors I am getting I do not seem to be making any progress. At the moment I am just trying to get the solver to run without error and will then plot this at a later date. Any advice / help would be greatly appreciated.

% This script produces a transient response of a 3DOF dynamic model of a
% jacket platform
clear all
clc
% Define global variables which will be used in the other functions
% global M k F I 
%% Input Parameters
E=204e+9;                       % Youngs modulus (GN/m^2)
R=1.25;                         % Outer radius of steel jacket (m)
r=1.23;                         % Inner radius of steel jacket (m)
M=9000e+3;                      % Deck mass (kg)
I=1.35e+9;                      % Inertia (kgm^2)
a= pi*(R^2-r^2);                % Cross-sectional area (m^2)
I2nd= (pi/4)*(R^4-r^4);         % Second moment of area (m^4)
l= [16 14 17 13] ;              % Length to centre of mass (m)
L = 20 ;                        % Leg Length (m)
Nrun = length(L);
omega = zeros(3, Nrun);         % Natural Frequencies
%% Defining forcing
% Wave specs : Height=2m ; Period=4s
fxa= 0.00;
fxb= 0.00;
fxc= 0.00;
fxd= 0.00;
fya= 371.56;
fyb= 378.96;
fyc= 371.56;
fyd= 351.44;
K=(3*E*I2nd) / (L^3) % Jacket stiffness (N/m)
% Define the mass, stiffness, and damping matirces m, k and c 
% and amplitude of force F.  We are using proportional damping
m=[  M 0 0 ; 0 M 0 ; 0 0 I ];
k= [ 4.*K 0 (2.*K*(l(1)-l(2))) ;  0 4.*K (2.*K*(l(4)-l(3))) ; 2.*K*(l(1)-l(2)) 2.*K*(l(4)-l(3)) 2.*K*(l(1)^2+l(2)^2+l(3)^2+l(4)^2) ];
c= 0.2*k;
F=[  fxa+fxb+fxc+fxd ; fya+fyb+fyc+fyd ; (fxa+fxb)*l(1)-(fxc+fxd)*l(2)-(fya+fyd)*l(3)+(fyb+fyc)*l(4) ];
% ODE Solution
% Initial Conditions
tspan = [0 1000] ;
z0 = [0 0 0 0 0 0] ;
% Solution
[t, z] = ode45(@(t, z) doffunc(t, z, k, M, I), tspan, z0) ;
%% ODE Function
function zdot=doffunc(t, z, k, M, I)
zdot = zeros(6, 1) ;
zdot(1)= z(4);
zdot(2)= z(5);
zdot(3)= z(6);
zdot(4)= ( - 0.2*((k(1)*z(4)) - (k(3)*z(6))) - (k(1)*z(1)) - (k(3)*z(3)) + F(1))/M;
zdot(5)= ( - 0.2*((k(5)*z(5)) - (k(6)*z(6))) - (k(5)*z(2)) - (k(6)*z(3)) + F(2))/M;
zdot(6)= ( - 0.2*((k(7)*z(4)) - (k(8)*z(5)) - (k(9)*z(6))) - (k(7)*z(1)) - (k(8)*z(2)) - (k(9)*z(3)) + F(3))/I; 
end

These are the errors that I am getting.

[Rena Berman]

(Answers Dev) Restored edit

Answer 1

[Jan]

The error message tells you, that Matlab fails in the function F. Unfortunately you post the other code, but not the failing function. Please edit the question and add this function also to allow to help you.

Or do you mean this vector?

F=[fxa+fxb+fxc+fxd; ...
   fya+fyb+fyc+fyd; ...
   (fxa+fxb)*l(1)-(fxc+fxd)*l(2)-(fya+fyd)*l(3)+(fyb+fyc)*l(4) ];

If so, include in in the list of provided parameters:

[t, z] = ode45(@(t, z) doffunc(t, z, k, M, I, F), tspan, z0) ;
...
function zdot=doffunc(t, z, k, M, I, F)

Then inside doffunc the values of the vector F is used, not the function stored in the File F.m .

[Maximilian]

Hello Jan,

Thank you so much I can't believe I missed this.

Thank you for your help

Answer 2

[Torsten]

% This script produces a transient response of a 3DOF dynamic model of a
% jacket platform
clear all
clc
% Define global variables which will be used in the other functions
% global M k F I 
%% Input Parameters
E=204e+9;                       % Youngs modulus (GN/m^2)
R=1.25;                         % Outer radius of steel jacket (m)
r=1.23;                         % Inner radius of steel jacket (m)
M=9000e+3;                      % Deck mass (kg)
I=1.35e+9;                      % Inertia (kgm^2)
a= pi*(R^2-r^2);                % Cross-sectional area (m^2)
I2nd= (pi/4)*(R^4-r^4);         % Second moment of area (m^4)
l= [16 14 17 13] ;              % Length to centre of mass (m)
L = 20 ;                        % Leg Length (m)
Nrun = length(L);
omega = zeros(3, Nrun);         % Natural Frequencies
%% Defining forcing
% Wave specs : Height=2m ; Period=4s
fxa= 0.00;
fxb= 0.00;
fxc= 0.00;
fxd= 0.00;
fya= 371.56;
fyb= 378.96;
fyc= 371.56;
fyd= 351.44;
K=(3*E*I2nd) / (L^3) % Jacket stiffness (N/m)

K = 9.1650e+06

% Define the mass, stiffness, and damping matirces m, k and c 
% and amplitude of force F.  We are using proportional damping
m=[  M 0 0 ; 0 M 0 ; 0 0 I ];
k= [ 4.*K 0 (2.*K*(l(1)-l(2))) ;  0 4.*K (2.*K*(l(4)-l(3))) ; 2.*K*(l(1)-l(2)) 2.*K*(l(4)-l(3)) 2.*K*(l(1)^2+l(2)^2+l(3)^2+l(4)^2) ];
c= 0.2*k;
F=[  fxa+fxb+fxc+fxd ; fya+fyb+fyc+fyd ; (fxa+fxb)*l(1)-(fxc+fxd)*l(2)-(fya+fyd)*l(3)+(fyb+fyc)*l(4) ];
% ODE Solution
    % Initial Conditions
    tspan = [0 10] ;
    z0 = [0 0 0 0 0 0] ;
    % Solution
    [t, z] = ode45(@(t, z) doffunc(t, z, k, M, I, F), tspan, z0) ;
    plot(t,z)

%% ODE Function
function zdot=doffunc(t, z, k, M, I, F)
zdot = zeros(6, 1) ;
zdot(1)= z(4);
zdot(2)= z(5);
zdot(3)= z(6);
zdot(4)= ( - 0.2*((k(1)*z(4)) - (k(3)*z(6))) - (k(1)*z(1)) - (k(3)*z(3)) + F(1))/M;
zdot(5)= ( - 0.2*((k(5)*z(5)) - (k(6)*z(6))) - (k(5)*z(2)) - (k(6)*z(3)) + F(2))/M;
zdot(6)= ( - 0.2*((k(7)*z(4)) - (k(8)*z(5)) - (k(9)*z(6))) - (k(7)*z(1)) - (k(8)*z(2)) - (k(9)*z(3)) + F(3))/I; 
end

[Maximilian]

Hello Torsten,

Thank you so much I can't believe I missed this.

Thank you for your help