How to solve 2nd order coupled differential equations?

I'd like to solve the 2 equation system below. I usually reduce the problem to a system of first order differencial equations and then use, for instance, ode113. However, here I don't know what to do because in each equation there are both the second derivatives of the variables.

Could anyone please help?

thanks a lot in advance

Answer 1

[Torsten]

M*x'' + A*x' + B*x = v
->
x'' = M^(-1) * (v - A*x' - B*x)
-> 
u1' = u2;
u2' = M^(-1) * (v - A*u2 - B*u1)

[Gloria]

Hi Torsten, thank you for your answer. Unfortunately I don't think this solution works because if

u1=x

u2=x'

u3=y

u4=y'

in the equation for u2' there would be also a (u4)'

u2' = M^(-1) * (v - A*u2 - B*u1-Cu4-Du4')

and for this reason I don't know how I could do

[Torsten]

u1 = [theta;phi];
u2 = [dtheta/dt;dphi/dt];
u = [u1;u2]
du/dt = [du1/dt;du2/dt] = [u2;M^(-1)*(v-A*u2-B*u1)]

4 differential equations in the unknown vector of functions u = [theta;phi;dtheta/dt;dphi/dt]

Answer 2

[Sam Chak]

Hi @Gloria,

If there are time-dependent elements in the mass matrix, the idea is to define the generalized 'Mass' matrix in odeset(). The system

can be rearranged to

,

where and can be defined as a state vector :

.

In the 1st-order ODE form, the system can be expressed in terms of state variables as

.

and rewritten compactly in the state-space representation:

,

where is the identity matrix and ℳ is the generalized mass matrix of the state-space that contains time-dependent elements.

Here is an example since we don't know the parameters of the system.

% Parameters
lr     = 1;
h      = 1;
md     = 1;
rd     = 1;
ms     = 1;
rs     = 1;
w      = 1;

% Generalized Mass matrix
M      = @(t) [eye(2), zeros(2); zeros(2), [lr+2*md*h^2+(2*md*rd^2+ms*rs^2)*(cos(w*t))^2, 1/2*(2*md*rd^2+ms*rs^2)*(sin(w*t))^2; 1/2*(2*md*rd^2+ms*rs^2)*(sin(w*t))^2, lr+2*md*h^2+(2*md*rd^2+ms*rs^2)*(sin(w*t))^2]];
opt    = odeset('Mass', M, 'MStateDependence', 'none');

% ode113 solver
tspan  = linspace(0, 30, 3001);
x0     = [4 -3 -2 1];
[t, x] = ode113(@odefcn, tspan, x0, opt);

% Solution Plot
plot(t, x, 'linewidth', 1.5), 
grid on, xlabel('t'), ylabel('\bf{x}')
legend({'$\theta$', '$\phi$', '$\dot{\theta}$', '$\dot{\phi}$'}, 'interpreter', 'latex', 'location', 'best', 'fontsize', 14)

function xdot = odefcn(t, x)    % right-hand side of state-space
    xdot      = zeros(4, 1);
    
    % Parameters
    c         = 1;
    dc        = 1;
    w         = 1;
    md        = 1;
    rd        = 1;
    ms        = 1;
    rs        = 1;
    Iz        = 1;
    k         = 1;
    dk        = 1;
    h         = 1;
    
    % Matrices
    C11       =   c*dc^2 - w*(2*md*rd^2 + ms*rs^2)*(sin(w*t))^2;
    C12       =   Iz*w + 2*w*(2*md*rd^2 + ms*rs^2)*(cos(w*t))^2;
    C21       = - Iz*w - 2*w*(2*md*rd^2 + ms*rs^2)*(sin(w*t))^2;
    C22       =   c*dc^2 + w*(2*md*rd^2 + ms*rs^2)*(sin(w*t))^2;
    C         = [C11, C12; C21, C22];               % damping matrix
    K         = diag([k*dk^2, k*dk^2]);             % stiffness matrix
    F         = 2*md*rd*h*w^2*[cos(w*t); sin(w*t)]; % force matrix
    
    % System
    xdot(1:2) = x(3:4);
    xdot(3:4) = F - C*x(3:4) - K*x(1:2);
end

[Gloria]

Hi @Sam Chak

Thanks a lot for your answer!

I'd have a question: why you don't define also C as a time-dependent matrix?

[Gloria]

and also in this equation:

xdot(3:4) = F - C*x(3:4) - K*x(1:2);

why there is not the matrix M?

[Sam Chak]

Hi @Gloria,

Q: I'd have a question: why you don't define also C as a time-dependent matrix?

The time-dependent matrix C is defined in odefcn() because it appears on the right-hand side of the rearranged equation.

Q: why there is not the matrix M?

Because the time-variant mass matrix M appears on the left-hand side of the rearranged equation, and I haven't resolved the coupling components, it is defined at the solver level. Now, after verifying that M is invertible at time , you can actually define M inside odefcn() and proceed as advised by @Torsten:

xdot(3:4) = M\(F - C*x(3:4) - K*x(1:2));

[Gloria]

Thank you very much @Sam Chak. I still have a doubt. In the equation for xdot(3) there should be xdot(4) and in the equation for xdot(4) there should be xdot(3), do you agree? But you wrote xdot(3:4) = M\(F - C*x(3:4) - K*x(1:2))

[Torsten]

x(1) = theta
x(2) = phi
x(3) = dtheta/dt
x(4) = dphi/dt

Thus

xdot(3:4) = d/dt ([dtheta/dt;dphi/dt]) = [d^2theta/dt^2;d^2phi/dt^2]
x(3:4) = [dtheta/dt;dphi/dt]
x(1:2) = [theta;phi]

so

[d^2theta/dt^2;d^2phi/dt^2] = M^(-1) * (F - C*[dtheta/dt;dphi/dt] - K*[theta;phi]) 

or

M * [d^2theta/dt^2;d^2phi/dt^2] + C*[dtheta/dt;dphi/dt] + K*[theta;phi] = F

[Sam Chak]

Hi @Gloria

If matrix remains on the left-hand side, then yes, and are still coupled.

However, when matrix is moved to the right-hand side of Eq. (1), it is essentially the same as solving a system of two algebraic equations:

,

,

where

.

In other words, and are considered to be fully decoupled in this form

,

which can be expressed in a relatively lengthy decoupled format:

So, the expression xdot(3:4) = M\(F - C*x(3:4) - K*x(1:2)) serves as a concise representation of the extensive decoupled form in MATLAB code, and it corresponds to Equation (2) mathematically.