Physics-informed neural networks (PINNs) use automatic differentiation to solve partial differential equations (PDEs) by penalizing the PDE in the loss function at a random set of points in the domain of interest. Here, we develop a Petrov-Galerkin version of PINNs based on the nonlinear approximation of deep neural networks (DNNs) by selecting the *trial space* to be the space of neural networks and the *test space* to be the space of Legendre polynomials. We formulate the *variational residual* of the PDE using the DNN approximation by incorporating the variational form of the problem into the loss function of the network and construct a *variational physics-informed neural network* (VPINN). By integrating by parts the integrand in the variational form, we lower the order of the differential operators represented by the neural networks, hence effectively reducing the training cost in VPINNs while increasing their accuracy compared to PINNs that essentially employ delta test functions. For shallow networks with one hidden layer, we analytically obtain explicit forms of the *variational residual*. We demonstrate the performance of the new formulation for several examples that show clear advantages of VPINNs over PINNs in terms of both accuracy and speed.