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

Question

Justin Burzachiello 2021-7-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/878893-how-can-i-model-second-order-ode-with-matrices-and-external-forcing

评论： Justin Burzachiello 2021-7-19

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Alan Stevens 2021-7-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/878893-how-can-i-model-second-order-ode-with-matrices-and-external-forcing#answer_747468

在 MATLAB Online 中打开

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

6 个评论
显示 4更早的评论隐藏 4更早的评论

Alan Stevens 2021-7-19

在 MATLAB Online 中打开

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