The Feynman-Kac formula enables one to represent the solutions to certain parabolic partial differential equations (PDEs) as integrals over the Wiener space; as a consequence, when using a Monte-Carlo approach to approximate the solution of a parabolic PDE, one is led to approximate an integral over the Wiener space. The present paper proposes a new approach to the approximation of such integrals, that is well adapted to situations where the coefficients of the underlying differential operator fail to satisfy Hörmander's condition. This new approach is based on cubature formulae, which can be best understood first in a finite-dimensional context: considering a positive measure \(\mu\) on \(\mathbb{R}^d\) such that \(\int | P(x)| d\mu (x)<\infty\) whenever \(P\) varies in the algebra \(\mathbb{R}_m[X_1,\ldots,X_d]\) of all polynomials in \(d\) variables having degree at most \(m\), Tchakaloff's theorem asserts that one may fix \(n\) points \(x_1,\ldots,x_n\in \text{supp}(\mu)\) and \(n\) positive weights \(\lambda_1,\ldots,\lambda_n\) in such a way that \[ \int_{\mathbb{R}^d} P(x)\,d\mu (x)=\sum_{i=1}^n\lambda_iP(x_i),\quad\forall P\in\mathbb{R}_m[X_1,\ldots,X_d]. \] For an arbitrary regular function \(f:\mathbb{R}^d\rightarrow\mathbb{R}\), one may then approximate the integral \(\int f(x)\,d\mu (x)\) by the sum \(\sum_{i=1}^n\lambda_if(x_i)\): the quality of this approximation depends on how ``close'' \(f\) is to a polynomial of degree at most \(m\), and this in turn may be investigated via the Taylor formula. Analogously, one may give a Taylor expansion for some regular functionals of the diffusion process corresponding to the parabolic PDE under consideration, and in this context, iterated Stratonovich integrals are playing the role earlier assigned to polynomials. Having recalled such stochastic Taylor expansions, the authors first state an analog of Tchakaloff's theorem for the Wiener measure, enabling one to replace the mean values of some iterated Stratonovich integrals under the Wiener measure by mean values taken with respect to a linear combination of the form \(\sum_{i=1}^n\lambda_i\delta_{\omega_i}\), the \(\omega_i\)'s being some fixed paths of bounded variation. Based on this theorem, an algorithm is then proposed to approximate the expected value of certain functionals of a diffusion, by ``solving many carefully chosen ODE's and taking a weighted average of these solutions''. An understanding of the algebraic relations connecting iterated Stratonovich integrals is necessary in order to give explicit constructions of cubature formulae on the Wiener space: such algebraic relations are exposed in a separate section, and two cubature formulae on the Wiener space are then given explicitly.