Second order differential equation with large matrices

Greetings,

I am tryig to solve this differential equation M * q'' + K * q= - K_d * q' - K_p * (q - q_d) where M is a 15x15 matrix, q is 15x1 vector, K is a 15x15 matrix, K_d and K_p are 15x15 known matrices and q_d is a 15x1 vector which is also known. q'' is a second time derivative. Which solver is the most applicable for this equation? and is there way to solve this in a matrix form and not expand the whole equation. Thanks in advance.

[Torsten]

q_d must be 15x1, not 8x1.

Rewrite your system as a system of first-order differential equations and use ode45, e.g., to solve it.

[David Togonidze]

typo.. my bad

[Jan]

As Torsten has explained alreayd, rewrite to 2nd order ODE to a system of order 1 equations. You can use matrices to evaluate this system. This is one of the fundamental design ideas of Matlab. Simply write the system as matrix equation.

[David Togonidze]

Thanks for the tips guys 👍

[David Togonidze]

One question guys. Let's say matrix M is not constant and is dependent on q, so at every time step M changes accordingly to the values of q calculated at previous time step (starting from initial condition q0). How can i incorporate that into the odesystem function for ode45 solver? Thank you in advance

[Torsten]

One question guys. Let's say matrix M is not constant and is dependent on q, so at every time step M changes accordingly to the values of q calculated at previous time step (starting from initial condition q0). How can i incorporate that into the odesystem function for ode45 solver? Thank you in advance

By dividing through M (assuming M is non-singular):

q'' = M \ (- K * q - K_d * q' - K_p * (q - q_d) )

If M is singular, define M as mass-matrix for the ODE solver in the options-structure and use ode15s or ode23t instead of ode45.

Answer 1

[Sam Chak]

Hi @David Togonidze

SInce this is a linear system, if you have the Control System Toolbox, then you can manipulate the matrix differential equation and convert it to a State-Space Model, that is similar to the system of 1st-order differential equations as advised by @Torsten, and then run the simulation using the lsim() command.

The following example demonstrates a 2-DOF Mass-Spring-Damper System.

M  = diag([1 1])

M = 2×2

1 0 0 1

Kd = 2*diag([1 1])

Kd = 2×2

2 0 0 2

Kp = 0.5*diag([1 1])

Kp = 2×2

0.5000 0 0 0.5000

K  = 0.5*diag([1 1])

K = 2×2

0.5000 0 0 0.5000

Convert the system of ODEs to State-space Model:

where the state vector is defined as

% State matrix
A  = [zeros(2) eye(2);
      M\(-Kp-K) M\(-Kd)]

A = 4×4

0 0 1 0 0 0 0 1 -1 0 -2 0 0 -1 0 -2

% Input matrix
% B  = [zeros(2); M\Kp]         % based on your original Kp*qd
B  = [zeros(2); M\(Kp+K)]       % modified to track desired qd

B = 4×2

0 0 0 0 1 0 0 1

% Output matrix
C  = [eye(2) zeros(2)]

C = 2×4

1 0 0 0 0 1 0 0

% Feedforward matrix
D  = 0;
sys = ss(A, B, C, D)

sys =
 
  A = 
       x1  x2  x3  x4
   x1   0   0   1   0
   x2   0   0   0   1
   x3  -1   0  -2   0
   x4   0  -1   0  -2
 
  B = 
       u1  u2
   x1   0   0
   x2   0   0
   x3   1   0
   x4   0   1
 
  C = 
       x1  x2  x3  x4
   y1   1   0   0   0
   y2   0   1   0   0
 
  D = 
       u1  u2
   y1   0   0
   y2   0   0
 
Continuous-time state-space model.

t  = 0:0.02:20;
u1 = max(0,min(t-1,1));     % qd1
u2 = max(0,min(t-1,1));     % qd2
U  = [u1; u2];              % qd array
lsim(sys, U, t)
% lsim(sys, U, t, x0)       % Use this if there is a vector x0 of initial state values
grid on

[David Togonidze]

Thank you a lot for such a detailed answer.. definitely will try this