Summary: The paper presents a randomized algorithm which evaluates the partition function of an arbitrary ferromagnetic Ising system to any specified degree of accuracy. The running time of the algorithm increases only polynomially with the size of the system (i.e., the number of sites) and a parameter which controls the accuracy of the result. Further approximation algorithms are presented for the mean energy and the mean magnetic moment of ferromagnetic Ising systems. The algorithms are based on Monte Carlo simulation of a suitably defined ergodic Markov chain. The states of the chain are not, as is customary, Ising spin configurations, but spanning subgraphs of the interaction graph of the system. It is shown that the expectations of simple operators on these configurations give numerical information about the partition function and related quantities. The performance guarantees for the algorithms are rigorously derived and rest on the fact that the Markov chain in question is rapidly mixing, i.e., converges to its equilibrium distribution in a polynomial number of steps. This is apparently the first time that rapid mixing has been demonstrated at all temperatures for a Markov chain related to the Ising model.