In this thesis we consider the stochastic p-Laplace and the symmetric p-Stokes system. The overall aim is to solve the systems, by which we mean that we find numerical values for the solution variable. Both models are related to non-Newtonian fluids, more specifically, power-law fluids. We discuss the derivation of fluid models for non-Newtonian fluids based on physical principles. In particular, we highlight the natural occurrence of stochasticity in the models, when the underlying physical quantities lack regularity. Mathematically, the results of the thesis are two folded. First, we address regularity results for stochastic p-Laplace and p-Stokes systems. Second, we use the regularity for the error analysis of a new numerical scheme for the p-Laplace system. The p-Laplace system arises naturally as the gradient flow of the p-Dirichlet energy. Exactly the gradient flow structure enables improved energy estimates. Exploiting the energy estimates, we show advanced temporal regularity results for strong solutions in optimal function spaces. By optimal we mean that the solution process has the same temporal regularity as the driving Wiener process. The p-Stokes system can also be interpreted as a gradient flow on the space of divergence free vector fields. We extend, at least partially, the results for the p-Laplace system to the p-Stokes system and prove optimal temporal regularity of the solution process. Regularity is the key ingredient in the a priori error analysis. Based on the established regularity we successfully bound the error of a new numerical scheme to the solution of the p-Laplace system with optimal rates in time and space. The new numerical scheme approximates time-averaged values of the solution variable. This leads to relaxed regularity requirements and improved robustness of the algorithm. In fact, our algorithm is stable even for possibly discontinuous solutions. Due to the approximation of time-averaged values we need to sample averaged Wiener increments. We supplement the algorithm by a simple sampling strategy that allows for an efficient implementation of the algorithm. We enrich our theoretical findings by suitable numerical simulations. In particular, we discuss different notions of error measures.