How can I represent the sparse mass, damping and stiffness matrices in state space?
3 次查看(过去 30 天)
显示 更早的评论
How can I represent the sparse mass, damping and stiffness matrices in state space?
0 个评论
回答(1 个)
David Goodmanson
2021-11-30
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']
or
[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);
figure(1)
plot(t,y); grid on
% case 2
[t y] = ode45(@(t,x) funM(t,x,n,M,C,K),[0 100], [1:2*n]);
figure(2)
plot(t,y); grid on
function dydt = funs(t,y,n,Cs,Ks)
dydt = [sparse(n,n) speye(n,n); -Cs -Ks]*y;
end
function dydt = funM(t,y,n,M,C,K)
dydt = [zeros(n,n) eye(n,n); (-M\C) (-M\K)]*y;
end
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;
end
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!