How can I represent the sparse mass, damping and stiffness matrices in state space?

3 views (last 30 days)
How can I represent the sparse mass, damping and stiffness matrices in state space?

Answers (1)

David Goodmanson
David Goodmanson on 30 Nov 2021
HI 'M',
Suppose you have n positions and n corresponding velocities. Let the state space (column) vector be in the order
y = [x1; x2; ....xn; x1'; x2'; ...xn'].
In matrix notation with mass, damping and stiffness being M, C and K (all nxn matrices), then the equation of motion is
M*x'' = -C*x'-K*x
Let I be the nxn identity matrix and 0 indicate an nxn matrix of zeros. Then
[I 0] * [x' ] = [ 0 I] * [x ]
[0 M] [x''] [-C -K] [x']
[I 0] * [y]' = [ 0 I] * [y] (a)
[0 M] [ ] [-C -K] [ ]
where the notation above represents 2nx2n matrices consisting of four nxn blocks. For the signs, C and K are intended to be positive in the sense that they have positive eigenvalues, so the solution oscillates and decays.
Taking M to the other side, he first order ode looks like
[y]' = [ 0 I ] * [y] (b)
[ ] [-M\C -M\K] [ ]
except that if M is sparse, you would not really want to be dividing by M and getting a full matrix as a result. Matlab odes have a feature that lets you keep M on the left using the odeset function. The first case below implements that, using sparse matrices everywhere. (In this example the matrices are anything but sparse, but it's an illustration that the sparse class works). The second case uses full matrices and verifies that dividing by M as in (b) agrees with the first case.
n = 5 % number of position variables
M = poseig(n);
Ms = sparse(M);
C = .01*poseig(n);
Cs = sparse(C);
K = poseig(n);
Ks = sparse(K);
M1s = [speye(n,n) sparse(n,n);sparse(n,n) Ms]; % as in (a)
opt = odeset('Mass',M1s);
% need 2n intitial conditions, arbitrarily choose 1:2n
[t y] = ode45(@(t,x) funs(t,x,n,Cs,Ks),[0 100],[1:2*n],opt);
plot(t,y); grid on
% case 2
[t y] = ode45(@(t,x) funM(t,x,n,M,C,K),[0 100], [1:2*n]);
plot(t,y); grid on
function dydt = funs(t,y,n,Cs,Ks)
dydt = [sparse(n,n) speye(n,n); -Cs -Ks]*y;
function dydt = funM(t,y,n,M,C,K)
dydt = [zeros(n,n) eye(n,n); (-M\C) (-M\K)]*y;
function m = poseig(n)
% create a symmetric matrix with positive eigenvalues
a = rand(n,n);
a = a+a';
[v lam] = eig(a);
m = (v*abs(lam))/v;

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by