The stochastic dynamics of networks of biochemical reactions in living cells are typically modelled using chemical master equations (CMEs). The stationary distributions of CMEs are seldom solvable analytically, and few methods exist that yield numerical estimates with computable error bounds. Here, we present two such methods based on mathematical programming techniques. First, we use semidefinite programming to obtain increasingly tighter upper and lower bounds on the moments of the stationary distribution for networks with rational propensities. Second, we employ linear programming to compute convergent upper and lower bounds on the stationary distributions themselves. The bounds obtained provide a computational test for the uniqueness of the stationary distribution. In the unique case, the bounds collectively form an approximation of the stationary distribution accompanied with a computable $\ell^1$-error bound. In the non-unique case, we explain how to adapt our approach so that it yields approximations of the ergodic distributions, also accompanied with computable error bounds. We illustrate our methodology through two biological examples: Schl\"ogl's model and a toggle switch model.