How can I model second order ODE with matrices and external forcing?

Hello. I am new to MATLAB-based modeling of ODEs, and I was wondering if someone here can help me with simulating the dynamics of the following ODE. Links to resources, code, and/or vocabulary would be appreciated.

M is the mass matrix, E is the damping matrix, K is the stiffness matrix, B describes inputs, and C describes outputs.

https://morwiki.mpi-magdeburg.mpg.de/morwiki/index.php/Linear_1D_Beam

My goal is the eventually model an impulse response or simply the evolution of the system with a collection of initial values.

For background, I don't really have sufficient background in FEA or numerical analysis. I am an intern learning the math and the MATLAB on the fly. I read into the documentation of ODE solvers, but I was unable to find a simple way to incorporate the equation above.

• Thank you

Answer 1

[Alan Stevens]

Like this perhaps (I've made up arbitrary data; you will obviously have to replace it with your own)

x0 = [1; -1];
v0 = [0; 0];
X0 = [x0; v0];

tspan = [0 1];

[t,X] = ode45(@fn, tspan, X0);
x = X(:,1:2);
v = X(:,3:4);

C = [2, 0; 0, 2];
y = C*x';

plot(t,x,'k',t,y,'r'),grid
xlabel('t'), ylabel('x (black) & y (red)')

function dXdt = fn(t,X)
    M = [1, 0.1; 0.5, 0.5];
    E = [2, 0; 0, 1];
    K = [1, 1; 0.1, 0.2];
    B = [1, 0; 0, 1];
    
    u = [1; 0];
    x = X(1:2);
    v = X(3:4);
    dxdt = v;
    dvdt = M\(B*u - K*x - E*v);
    
    dXdt = [dxdt; dvdt];

end

[Justin Burzachiello]

Thank you Alan. I have not yet confirmed if this works for my exact matrices, but this certainly helps me develop the code.

[Justin Burzachiello]

Also, I actually just wrote a follow-up question to this one in https://www.mathworks.com/matlabcentral/answers/879468-how-read-matrix-market-data-into-matlab since I am struggling with figuring out how to import Matrix Market data. I just wanted to mention this in case you are familiar with the topic. Thanks for your help!

[Alan Stevens]

See https://uk.mathworks.com/matlabcentral/answers/879468-how-read-matrix-market-data-into-matlab now.

[Justin Burzachiello]

Hi Alan. I tried altering your code above to work with the actual matrices from my data, but I am having trouble. How would you alter the code above if the matrices have the following dimensions?

E is 18x18

M is 18x18

K is 18x18

B is 18x1

C is 1x18

I made u a 1x18 vector of ones, while x and xDot are both 18x1 vectors of zeros. I get the error from the attached picture. My current code:

x0 = zeros(18,1);
xd0 = zeros(18,1);
X0 = [x0; xd0];

tspan = [0 10];

%error occurs when calling ode45
[t, X] = ode45(@fn, tspan, X0);

function dXdt = fn(t, X)
    u = ones(1,18);
    x = X(1:18, 1);
    xd = X(19:36, 1);
    
    [M, ~, ~, ~] = mmread('LF10_M.m');
    [E, ~, ~, ~] = mmread('LF10_E.m');
    [K, ~, ~, ~] = mmread('LF10_K.m');
    [B, ~, ~, ~] = mmread('LF10_B.m');
    
    M = full(M);
    E = full(E);
    K = full(K);
    B = full(B);
    
    xdd = M\(B*u - K*x - E*xd);
    
    dXdt = [xd; xdd];
end

[Alan Stevens]

Like this (no point in calling the files every time fn is called by the ode function):

[M, ~, ~, ~] = mmread('LF10_M.m');
[E, ~, ~, ~] = mmread('LF10_E.m');
[K, ~, ~, ~] = mmread('LF10_K.m');
[B, ~, ~, ~] = mmread('LF10_B.m');

M = full(M);
E = full(E);
K = full(K);
B = full(B);
  
x0 = zeros(18,1);
xd0 = zeros(18,1);
X0 = [x0; xd0];

tspan = [0 1];

%error occurs when calling ode45
[t, X] = ode15s(@(t,X) fn(t,X,M,E,K,B), tspan, X0);   %%%%% ode15s works best

plot(t,X(:,1:18))

function dXdt = fn(t, X,M,E,K,B)
    u = ones(18,1);                  %%%%%%%%%%%%
    x = X(1:18, 1);
    xd = X(19:36, 1);
    
    xdd = M\(B.*u - K*x - E*xd);     %%%%%%%%%%%%%%  B.*u
    
    dXdt = [xd; xdd];
end

[Justin Burzachiello]

Thank you, Alan. Also, thanks for the help with inputting the data matrices as arguments.