How can I model second order ODE with matrices and external forcing?
1 view (last 30 days)
Show older comments
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.
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 Comments
Answers (1)
Alan Stevens
on 15 Jul 2021
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 Comments
Alan Stevens
on 19 Jul 2021
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
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!