There is a mathematical formula, but it requires knowing the pressure and velocity distribution around the surface of the object:
$$D=\int_{S_{upper}}\left[-pcos(\theta)+\tau_wsin(\theta)\right]dA +\int_{S_{lower}}\left[psin(\theta)+\tau_wcos(\theta)\right]dA$$
where $S_{lower}$ and $S_{upper}$ refer to the lower and upper surfaces, respectively, and $\theta$ is the angle between the local unit normal vector to the surface and the vertical direction.
To get the pressure and velocity under specific conditions, you need one of three things: experiment data, a computational simulation, or an analytical solution to the navier stokesNavier Stokes equations for your particular case of interest. Analytical solutions are few and hard to come by for the general case with arbitrary flow domain shapes. Experiment data can be costly to obtain and would require using instruments that necessarily disrupt the flow, making it more difficult to accurately measure what you want (although, well designed experiments minimize disruption as much as possible).
Computational simulation, in my opinion, is the best compromise between the other two. It provides a "virtual" experiment using analytical equations that need to be solved iteratively. Nonetheless, as others have stated, it is non-trivial to obtain useful solutions from computational simulation. Experienced CFD analysts spend a lot of time preparing cases carefully, creating appropriate meshes, using/implementing proper gradient schemes, numerical fluxes, time stepping schemes, turbulence models, transition models, shock tailoring, boundary layer treatment, etc... While CFD may be more ubiquitously accessible, it still takes a lot more effort to produce useful solutions.
When computers were still relatively low power and were not yet capable of solving navier stokesNavier Stokes equations, there were other methods for estimating the flow field parameters. One such method is to first solve the "Euler Equations", which are essentially an inviscid form of the navier stokes equations. Of course, if the fluid is only treated as inviscid, the drag would always be zero. To estimate the viscous solution from the inviscid one, we look at the individual streamlines of the inviscid solution on the body of interest and integrate from start to finish to sort of "build" a boundary layer up iteratively until a convergence criterion is met. Once you've built up enough of a boundary layer thickness, you can use various correlations to obtain the resulting pressure $p$ and wall shear stress $\tau_w$ that you need to compute drag. This method is often referred to in literature as the "Momentum Integral Equation".