Let m,n≥2 be positive integers. Denote by Mm the set of m×m complex matrices and by w(X) the numerical radius of a square matrix X. Motivated by the study of operations on bipartite systems of quantum states, we show that a linear map ϕ:Mmn→Mmn satisfies w(ϕ(A⊗B))=w(A⊗B)for all A∈Mm and B∈Mn if and only if there is a unitary matrix U∈Mmn and a complex unit ξ such that ϕ(A⊗B)=ξU(φ1(A)⊗φ2(B))U∗for all A∈Mm and B∈Mn, where φk is the identity map or the transposition map X↦Xt for k=1,2, and the maps φ1 and φ2 will be of the same type if m,n≥3. In particular, if m,n≥3, the map corresponds to an evolution of a closed quantum system (under a fixed unitary operator), possibly followed by a transposition. The results are extended to multipartite systems.