The third vertical derivative of gravity anomalies is applied in gravity prospecting for estimating the parameters of anomalous geological bodies and in physical geodesy for computing the third correction term within the framework of solving the M. S. Molodensky problem using the method of analytic continuation. The primary objective of this study was to develop and analyze algorithms for calculating this derivative by two distinct methods: (1) numerical differentiation based on a five-point finite difference scheme, wherein gravity anomalies are specified along the vertical Oz axis relative to a reference point with a constant height step, and (2) an analytical approach based on expansion into a truncated series of spherical harmonics utilizing coefficients from the global geopotential model EIGEN-6C4 up to degree N = 2190. Implementation of the developed algorithms and corresponding software yielded two global digital models of the third vertical derivative of gravity anomalies. The analysis demonstrated that at solving the M. S. Molodensky problem with a third-order analytic approximation, computation of this derivative is more economically efficient when performed analytically, i.e. via expansion into a truncated series of spherical functions compared to numerical differentiation techniques