A stiffness matrix H′ is designed as PiTCPi or a linear combination form thereof, where a spatial first-order difference operator, which operates on a state vector, of a displacement vector at each discrete grid point is Pi and the matrix of elastic constants is C. The mass matrix T′ is derived using finite-difference approximations at the discrete grid points so that the result of the matrix δTnum operating on a state vector corresponding to rotations or translations is 0, where T is a positive definite diagonal mass matrix and δTnum=T′−T.