Mathematical Analysis of Some Typical Problems in Geodesy by Means of Computer Algebra

Abstract

There are many complicated and fussy mathematical analysis processes in geodesy, such as the power series expansions of the ellipsoid’s eccentricity, high order derivation of complex and implicit functions, operation of trigonometric function, expansions of special functions and integral transformation. Taking some typical mathematical analysis processes in geodesy as research objects, the computer algebra analysis are systematically carried out to bread, deep and detailed extent with the help of computer algebra analysis method and the powerful ability of mathematical analysis of computer algebra system. The forward and inverse expansions of the meridian arc in geometric geodesy, the nonsingular expressions of singular integration in physical geodesy and the series expansions of direct transformations between three anomalies in satellite geodesy are established, which have more concise form, stricter theory basis and higher accuracy compared to traditional ones. The breakthrough and innovation of some mathematical analysis problems in the special field of geodesy are realized, which will further enrich and perfect the theoretical system of geodesy.

1. Introduction

Geodesy is the science of accurately measuring and understanding three fundamental properties of the Earth: its geometric shape, its orientation in space, and its gravity field, as well as the changes of these properties with time. There are many fussy symbolic problems to be dealt with in geodesy, such as the power series expansions of the ellipsoid’s eccentricity, high order derivation of complex and implicit functions, expansions of special functions and integral transformation. Many geodesists have made great efforts to solve these problems, see [1, 2, 3, 4, 5, 6, 7, 8]. Due to historical condition limitation, they mainly disposed of these problems by hand, which were not perfectly solved yet. Traditional algorithms derived by hand mainly have the following problems: (1) The expressions are complex and lengthy, which makes the computation process very complicated and time-consuming. (2) Some approximate disposal is adopted, which influences the computation accuracy. (3) Some formulae are numerical and only apply to a specific reference ellipsoid, which are not convenient to be generalized.

In computational mathematics, computer algebra, also called symbolic computation, is a scientific area that refers to the study and development of algorithms and software for manipulating mathematical expressions and other mathematical objects. Software applications that perform symbolic calculations are called computer algebra systems, which are more popular today. Computer algebra systems, like Mathematica, Maple and Mathcad, are powerful tools to solve some mathematical derivation in geodesy, see [9, 10, 11]. By means of computer algebra analysis method and computer algebra system Mathematica, we have already solved many complicated mathematical problems in special fields of geodesy in the past few years; see [12, 13, 14, 15].

The main contents and research results presented in this chapter are organized as follows: In Section 2, we discuss the forward and inverse expansions of the meridian arc often used in geometric geodesy. In Section 3, the nonsingular expressions of singular integration in physical geodesy are derived. In Section 4, we discuss series expansions of direct transformations between three anomalies in satellite geodesy. Finally in Section 5, we make a brief summary.

2. The forward and inverse expansions of the meridian arc in geometric geodesy

The forward and inverse problem of the meridian arc is one of the fundamental problems in geometric geodesy, which seems to be a solved one. Briefly reviewing the existing methods, however, one will find that the inverse problem was not perfectly and ideally solved yet. This situation is due to the complexity of the problem itself and the lack of advanced computer algebra systems. Yang had given the direct expansions of the inverse transformation by means of the Lagrange series method, but their coefficients are expressed as polynomials of coefficients of the forward expansions, which are not convenient for practical use, see [6, 7]. Adams expressed the coefficients of inverse expansions as a power series of the eccentricity $e$ by hand, but expanded them up to eighth-order terms of $e$ at most, see [1]. Due to these reasons, the forward and inverse expansions of the meridian arc are discussed by means of Mathematica in the following sections. Their coefficients are uniformly expressed as a power series of the eccentricity and extended up to tenth-order terms of $e$.

2.1. The forward expansion of the meridian arc

The meridian arc from the equator where the latitude is from $B=0$ to $B$ is

$$X=a\left(1-e^2\right){\int }_0^B{\left(1-e^2{sin}^{2}B\right)}^{-3/2}\mathit{dB}$$E1

where $X$ is the meridian arc; $B$is the geodetic latitude; $a$ is the semi-major axis of the reference ellipsoid; $e$ is the first eccentricity of the reference ellipsoid.

Expanding the integrand in Eq. (1) and integrating it item by item using Mathematica as follows:

Then one arrives at

$$X=a\left(1-e^2\right)\left(K_{0}B+K_{2}sin2B+K_{4}sin4B+K_{6}sin6B+K_{8}sin8B+K_{10}sin10B\right)$$E2

where

$$\left\{\begin{array}{l}{K}_0=1+\frac{3}{4}{e}^2+\frac{45}{64}{e}^4+\frac{175}{256}{e}^6+\frac{11025}{16384}{e}^8+\frac{43659}{65536}{e}^{10}\\ K_2=-\frac{3}{8}{e}^2-\frac{15}{32}{e}^4-\frac{525}{1024}{e}^6-\frac{2205}{4096}{e}^8-\frac{72765}{131072}{e}^{10}\\ K_4=\frac{15}{256}{e}^4+\frac{105}{1024}{e}^6+\frac{2205}{16384}{e}^8+\frac{10395}{65536}{e}^{10}\\ K_6=-\frac{35}{3072}{e}^6-\frac{105}{4096}{e}^8-\frac{10395}{262144}{e}^{10}\\ K_8=\frac{315}{131072}{e}^8+\frac{3465}{524288}{e}^{10}\\ K_{10}=-\frac{693}{1310720}{e}^{10}\end{array}\right.$$E3

Eqs. (2) and (3) can also be derived using the binomial theorem by hand which consumes much more time, see [6, 7, 8, 9]. The denominator 693 of the last coefficient $K_{10}$ is mistaken as 639 in Ref. [9].

2.2. The inverse expansion of the meridian arc using the Hermite interpolation method

Differentiation to the both sides of Eq. (1) yields:

$$\frac{\mathit{dX}}{\mathit{dB}}=\frac{a\left(1-e^2\right)}{{\left(1-e^2{sin}^{2}B\right)}^{3/2}}$$E4

Define $\psi$ as

$$\psi =\frac{X}{a\left(1-e^2\right)K_0}$$E5

Substituting Eq. (5) into Eq. (4) yields:

$$\frac{\mathit{dB}}{K_0\mathit{d\psi }}={\left(1-e^2{sin}^{2}B\right)}^{3/2}$$E6

Suppose that the inverse solution of Eq. (6) has the following form:

$$B=\psi +a_{2}sin2\psi +a_{4}sin4\psi +a_{6}sin6\psi +a_{8}sin8\psi +a_{10}sin10\psi$$E7

Eq. (7) has five coefficients to be determined. Once these coefficients are known, the inverse problem can be solved.

We consider that the values of differentiation Eq. (6) at the beginning and end points can be treated as interpolation constraints. It is generally known that

$$B^{\prime }\left(0\right)=K_0$$E8

$$B^{\prime }\left(\frac{\pi }{2}\right)=K_0{\left(1-e^2\right)}^{3/2}$$E9

The further derivation of Eq. (6) with respect to $\psi$ yields $B^{″}\left(\psi \right)$. Unfortunately, $B^{″}\left(\psi \right)$ is equal to zero at $\psi =0$, $\psi =\frac{\pi }{2}$. Hence, one differentiates Eq. (6) twice and it yields $B^{‴}\left(\psi \right)$. Omitting the derivative procedure by means of Mathematica, one arrives at

$$B^{‴}\left(0\right)=-3e^2-\frac{27}{4}{e}^4-\frac{729}{64}{e}^6-\frac{4329}{256}{e}^8-\frac{381645}{16384}{e}^{10}$$E10

$$B^{‴}\left(\frac{\pi }{2}\right)=3e^2-\frac{15}{4}{e}^4+\frac{57}{64}{e}^6+\frac{3}{256}{e}^8-\frac{51}{16384}{e}^{10}$$E11

The further derivation of $B^{‴}\left(\psi \right)$ with respect to $\psi$ yields $B^{\left(4\right)}\left(\psi \right)$, but $B^{\left(4\right)}\left(\psi \right)$ is equal to zero at $\psi =0$, $\psi =\frac{\pi }{2}$. Hence, one differentiates $B^{‴}\left(\psi \right)$ twice and it yields$B^{\left(5\right)}\left(\psi \right)$. Then one arrives at

$$B^{\left(5\right)}\left(0\right)=12e^2+90e^4+\frac{4455}{16}{e}^6+\frac{20145}{32}{e}^8+\frac{4924935}{4096}{e}^{10}$$E12

Making use of the five interpolation constraints in Eqs. (8–12) and differentiating Eq. (7) correspondingly, one arrives at a set of linear equations for the unknown coefficients

$$\left(\begin{array}{lllll}2& 4& 6& 8& 10\\ -2& 4& -6& 8& -10\\ -8& -64& -216& -512& -1000\\ 8& -64& -216& -512& 1000\\ 32& 1024& 7776& 32768& 100000\end{array}\right)\left(\begin{array}{c}{a}_2\\ a_4\\ a_6\\ \begin{array}{l}{a}_8\\ a_{10}\end{array}\end{array}\right)=\left(\begin{array}{c}{B}^{\prime }\left(0\right)-1\\ B^{\prime }\left(\frac{\pi }{2}\right)-1\\ B^{‴}\left(0\right)\\ \begin{array}{l}{B}^{‴}\left(\frac{\pi }{2}\right)\\ B^{\left(5\right)}\left(0\right)\end{array}\end{array}\right)$$E13

The solution to Eq. (13) is

$$\left(\begin{array}{c}{a}_2\\ a_4\\ a_6\\ \begin{array}{l}{a}_8\\ a_{10}\end{array}\end{array}\right)={\left(\begin{array}{lllll}2& 4& 6& 8& 10\\ -2& 4& -6& 8& -10\\ -8& -64& -216& -512& -1000\\ 8& -64& -216& -512& 1000\\ 32& 1024& 7776& 32768& 100000\end{array}\right)}^{-1}\left(\begin{array}{c}{B}^{\prime }\left(0\right)-1\\ B^{\prime }\left(\frac{\pi }{2}\right)-1\\ B^{‴}\left(0\right)\\ \begin{array}{l}{B}^{‴}\left(\frac{\pi }{2}\right)\\ B^{\left(5\right)}\left(0\right)\end{array}\end{array}\right)$$E14

Omitting the main operations by means of Mathematica, one arrives at

$$\{\begin{array}{l}{a}_2=\frac{3}{8}{e}^2+\frac{3}{16}{e}^4+\frac{213}{2048}{e}^6+\frac{255}{4096}{e}^8+\frac{20861}{524288}{e}^{10}\hfill \\ a_4=\frac{21}{256}{e}^4+\frac{21}{256}{e}^6+\frac{533}{8192}{e}^8+\frac{197}{4096}{e}^{10}\hfill \\ a_6=\frac{151}{6144}{e}^6+\frac{151}{4096}{e}^8+\frac{5019}{131072}{e}^{10}\hfill \\ a_8=\frac{1097}{131072}{e}^8+\frac{1097}{65536}{e}^{10}\hfill \\ a_{10}=\frac{8011}{2621440}{e}^{10}\hfill \end{array}$$E15

2.3. The inverse expansion of the meridian arc using Lagrange’s theorem method

Suppose that

$$y=x+f\left(x\right)$$E16

with $\left|f\left(x\right)\right|<<\left|x\right|$ and $y\approx x$. Lagrange’s theorem states that in a suitable domain the solution of Eq. (16) is

$$x=y+\sum _{n=1}^{\infty }\frac{{\left(-1\right)}^n}{n!}\frac{d^{n-1}}{{\mathit{dy}}^{n-1}}{\left[f\left(y\right)\right]}^n$$E17

Suppose$f\left(x\right)$is defined by the following finite trigonometric series:

$$f\left(x\right)=\alpha sin2x+\beta sin4x+\gamma sin6x+\delta sin8x+\epsilon sin10x$$E18

where the coefficients $\alpha =O\left(e^2\right)$, $\beta =O\left(e^4\right)$, $\gamma =O\left(e^6\right)$, $\delta =O\left(e^8\right)$, $\epsilon =O\left(e^{10}\right)$ are small enough for the condition $\left|f\left(x\right)\right|<<\left|x\right|$. In deriving the inversion we shall truncate the infinite Lagrange expansion at eighth-order terms of $e$ and drop higher powers. Inserting Eq. (18) into Eq. (17), one arrives at

$$x=y-f\left(y\right)+\frac{1}{2!}\frac{d}{\mathit{dy}}{\left[f\left(y\right)\right]}^2-\frac{1}{3!}\frac{d^2}{{\mathit{dy}}^2}{\left[f\left(y\right)\right]}^3+\frac{1}{4!}\frac{d^3}{{\mathit{dy}}^3}{\left[f\left(y\right)\right]}^4-\frac{1}{5!}\frac{d^4}{{\mathit{dy}}^4}{\left[f\left(y\right)\right]}^5$$E19

One should make use of several trigonometric identities to calculate the derivatives, substituting them into Eq. (19) and grouping terms according to the trigonometric functions. It is a difficult and time-consuming work to do by hand, but could be easily realized by means of Mathematica. Omitting the main procedure, one arrives at

$$x=y+d_{2}sin2y+d_{4}sin4y+d_{6}sin6y+d_{8}sin8y+d_{10}sin10y$$E20

where

$$\left\{\begin{array}{l}{d}_2=-\alpha -\mathit{\alpha \beta }-\mathit{\beta \gamma }+\frac{1}{2}{\alpha }^3+{\mathit{\alpha \beta }}^2-\frac{1}{2}{\alpha }^2\gamma +\frac{1}{3}{\alpha }^3\beta -\frac{1}{12}{\alpha }^5\\ d_4=-\beta +{\alpha }^2-2\mathit{\alpha \gamma }+4{\alpha }^2\beta -\frac{4}{3}{\alpha }^4\\ d_6=-\gamma +3\mathit{\alpha \beta }-3\mathit{\alpha \delta }-\frac{3}{2}{\alpha }^3+\frac{9}{2}{\mathit{\alpha \beta }}^2+9{\alpha }^2\gamma -\frac{27}{2}{\alpha }^3\beta +\frac{27}{8}{\alpha }^5\\ d_8=-\delta +2{\beta }^2+4\mathit{\alpha \gamma }-8{\alpha }^2\beta +\frac{8}{3}{\alpha }^4\\ d_{10}=-\epsilon +5\mathit{\alpha \delta }+5\mathit{\beta \gamma }-\frac{25}{2}{\mathit{\alpha \beta }}^2-\frac{25}{2}{\alpha }^2\gamma +\frac{125}{6}{\alpha }^3\beta -\frac{125}{24}{\alpha }^5\end{array}\right.$$E21

Substituting $x$ for $B$and $y$ for $\psi$ in Eq. (16), the coefficients $\alpha$, $\beta$, $\gamma$, $\delta$, $\epsilon$ are consistent with ${\alpha }_2$, ${\alpha }_4$, ${\alpha }_6$, ${\alpha }_8$ in Eq. (3). According to Eq. (20) and denoting $a_2$, $a_4$, $a_6$, $a_8$, $a_{10}$ as the new coefficients, the inverse expansion of the meridian arc can be written as

$$B=\psi +a_{2}sin2\psi +a_{4}sin4\psi +a_{6}sin6\psi +a_{8}sin8\psi +a_{10}sin10\psi$$E22

where

$$\left\{\begin{array}{l}{a}_2=-{\alpha }_2-{\alpha }_2{\alpha }_4-{\alpha }_4{\alpha }_6+\frac{1}{2}{\alpha }_2^3+{\alpha }_2{\alpha }_4^2-\frac{1}{2}{\alpha }_2^2{\alpha }_6+\frac{1}{3}{\alpha }_2^3{\alpha }_4-\frac{1}{12}{\alpha }_2^5\\ a_4=-{\alpha }_4+{\alpha }_2^2-2{\alpha }_2{\alpha }_6+4{\alpha }_2^2{\alpha }_4-\frac{4}{3}{\alpha }_2^4\\ a_6=-{\alpha }_6+3{\alpha }_2{\alpha }_4-3{\alpha }_2{\alpha }_8-\frac{3}{2}{\alpha }_2^3+\frac{9}{2}{\alpha }_2{\alpha }_4^2+9{\alpha }_2^2{\alpha }_6-\frac{27}{2}{\alpha }_2^3{\alpha }_4+\frac{27}{8}{\alpha }_2^5\\ a_8=-{\alpha }_8+2{\alpha }_4^2+4{\alpha }_2{\alpha }_6-8{\alpha }_2^2{\alpha }_4+\frac{8}{3}{\alpha }_2^4\\ a_{10}=-{\alpha }_{10}+5{\alpha }_2{\alpha }_8+5{\alpha }_4{\alpha }_6-\frac{25}{2}{\alpha }_2{\alpha }_4^2-\frac{25}{2}{\alpha }_2^2{\alpha }_6+\frac{125}{6}{\alpha }_2^3{\alpha }_4-\frac{125}{24}{\alpha }_2^5\end{array}\right.$$E23

The coefficients in Eq. (23) are also easily expressed in a power series of the eccentricity by means of Mathematica. Omitting the main procedure, one arrives at

$$\left\{\begin{array}{l}{a}_2=\frac{3}{8}{e}^2+\frac{3}{16}{e}^4+\frac{213}{2048}{e}^6+\frac{255}{4096}{e}^8+\frac{20861}{524288}{e}^{10}\begin{array}{c}\end{array}\\ a_4=\frac{21}{256}{e}^4+\frac{21}{256}{e}^6+\frac{533}{8192}{e}^8+\frac{197}{4096}{e}^{10}\begin{array}{c}\end{array}\begin{array}{ccc}& & \begin{array}{c}\end{array}\end{array}\\ a_6=\frac{151}{6144}{e}^6+\frac{151}{4096}{e}^8+\frac{5019}{131072}{e}^{10}\begin{array}{cccc}\begin{array}{cc}& \end{array}& & & \end{array}\\ \begin{array}{l}{a}_8=\frac{1097}{131072}{e}^8+\frac{1097}{65536}{e}^{10}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}\begin{array}{cc}\begin{array}{ccc}\begin{array}{cc}& \end{array}& & \end{array}& \end{array}\\ a_{10}=\frac{8011}{2621440}{e}^{10}\end{array}\end{array}\right.$$E24

The results in Eq. (24) are consistent with the coefficients in Eq. (15), which substantiates the correctness of the derived formula.

2.4. The accuracy of the inverse expansion of the meridian arc

In order to validate the exactness of the inverse expansions of meridian arc derived by the author, one has examined its accuracy choosing the CGCS2000 reference ellipsoid with parameters$a=6378137$, $e=0.08181919104281579$. The accuracy of the inverse expansions derived by Yang (see [6, 7]) is also examined for comparison.

One makes use of Eq. (1) and Eq. (5) to obtain the theoretical value ${\psi }_0$at given geodetic latitude $B_0$. Then one makes use of the inverse expansions derived by Yang (see [6, 7]) to obtain the computation value $B_1$. Substituting ${\psi }_0$into Eq. (22), one arrives at the computation value $B_1^{\prime }$. The differences $\mathrm{\Delta }{B}_i=B_i-B_0$, $\mathrm{\Delta }{B}_i^{\prime }=B_i^{\prime }-B_0\left(i=1,2\right)$ indicate the accuracies of the inverse expansions derived by Yang (see [6, 7]) and the author respectively. These errors are listed in Table 1.

$B_0/\left({}^{\circ }\right)$	$20$	$40$	$60$	$80$
$\mathrm{\Delta }{B}_1/\left({}^{″}\right)$	$-3.2\times {10}^{-7}$	$-1.6\times {10}^{-6}$	$-1.8\times {10}^{-6}$	$-1.4\times {10}^{-6}$
$\mathrm{\Delta }{B}_1^{\prime }/\left({}^{″}\right)$	$2.7\times {10}^{-9}$	$8.6\times {10}^{-9}$	$1.3\times {10}^{-8}$	$1.7\times {10}^{-8}$

Table 1.

Errors of the inverse expansions of meridian arc.

From Table 1, one could find that the accuracy of the inverse expansion of meridian arc derived by Yang (see [6, 7]) is higher than 10⁻⁵″, and the accuracy of the inverse expansion Eq. (22) derived by the author is higher than 10⁻⁷″. The accuracy is improved by 2 orders of magnitude by means of computer algebra.

3. The singular integration in physical geodesy

Singular integrals associated with the reciprocal distance usually exist in the computations of physical geodesy and geophysics. For example, the integral expressions of height anomaly, deflections of the vertical and vertical gradient of gravity anomaly can be written in planar approximation as

$$\zeta =\frac{1}{2\mathit{\pi \gamma }}\iint \frac{\mathrm{\Delta }g}{r}\mathit{\text{dxdy}}$$E25

$$\xi =-\frac{1}{2\mathit{\pi \gamma }}\iint \frac{x\mathrm{\Delta }g}{r^3}\mathit{\text{dxdy}}$$E26

$$\eta =-\frac{1}{2\mathit{\pi \gamma }}\iint \frac{y\mathrm{\Delta }g}{r^3}\mathit{\text{dxdy}}$$E27

$$L=\frac{1}{2\pi }\iint \frac{\mathrm{\Delta }g-\mathrm{\Delta }{g}_0}{r^3}\mathit{\text{dxdy}}$$E28

where $r=\sqrt{x^2+y^2}$, $\mathrm{\Delta }{g}_0$ is the gravity anomaly at the computation point. When $r\to 0$, the above integrals become singular and need special treatment at the computation point. The past treatments are with respect to template computations, which regards the innermost area as a circle, see [3, 16]. But the gravity anomalies are given in the rectangular grids(such as $2^{\prime }\times 2^{\prime }$), if the approximate disposal is used, some significant error may be introduced. Sunkel and Wang expressed the gravity anomalies block by block as an interpolation polynomial and derived the analytic values of the integrals, see [17, 18]. However, the integrals of the rational functions are very complicated, especially when related interpolation polynomials contain many terms. One only can give the analytic values of the corresponding linear approximation. In this chapter, we use the nonsingular integration transformations proposed by Bian (see [19]) to compute the above integrals precisely.

As shown in Figure 1, let the innermost area be the rectangular $\sigma \in \left[-a<x<a,-b<y<b\right]$ ($a>0$, $b>0$) due to the convergence of meridian, and the gravity anomaly is expressed as double quadratic polynomial:

$$\mathrm{\Delta }g\left(x,y\right)=\sum _{i=0}^2{\left(\frac{x}{a}\right)}^i\sum _{j=0}^2{\alpha }_{\mathit{ij}}{\left(\frac{y}{b}\right)}^j$$E29

Figure 1.

Integrals in the rectangular area.

Inserting the innermost area into Eqs. (25)–(28), and let the contributions be $\mathrm{\Delta }\zeta$, $\mathrm{\Delta }\xi$, $\mathrm{\Delta }\eta$and $\mathrm{\Delta }L$, one arrives at

$$\mathrm{\Delta }\zeta =\frac{1}{2\mathit{\pi \gamma }}\underset{\sigma }{\iint }\frac{\mathrm{\Delta }g\left(x,y\right)}{r}\mathit{\text{dxdy}}$$E30

$$\mathrm{\Delta }\xi =-\frac{1}{2\mathit{\pi \gamma }}\underset{\sigma }{\iint }\frac{x\mathrm{\Delta }g\left(x,y\right)}{r^3}\mathit{\text{dxdy}}$$E31

$$\mathrm{\Delta }\eta =-\frac{1}{2\mathit{\pi \gamma }}\underset{\sigma }{\iint }\frac{y\mathrm{\Delta }g\left(x,y\right)}{r^3}\mathit{\text{dxdy}}$$E32

$$\mathrm{\Delta }L=\frac{1}{2\pi }\underset{\sigma }{\iint }\frac{\mathrm{\Delta }g\left(x,y\right)-\mathrm{\Delta }g\left(0,0\right)}{r^3}\mathit{\text{dxdy}}$$E33

The following transformation is introduced for $\sigma$

$$\left\{\begin{array}{c}u=\frac{x}{a}\\ v=\frac{y}{b}\end{array}\right.$$E34

Using the properties of the integration for even/odd functions and exploiting the symmetry of the integration area, one arrives at

$$\mathrm{\Delta }\zeta =\frac{4\mathit{ab}}{2\mathit{\pi \gamma }}{\int }_0^1{\int }_0^1\frac{{\alpha }_{00}+{\alpha }_{20}{u}^2+{\alpha }_{02}{v}^2+{\alpha }_{22}{u}^2{v}^2}{\sqrt{a^2{u}^2+b^2{v}^2}}\mathit{\text{dudv}}$$E35

$$\mathrm{\Delta }\xi =\frac{-4a^{2}b}{2\mathit{\pi \gamma }}{\int }_0^1{\int }_0^1\frac{u^2\left({\alpha }_{10}+{\alpha }_{12}{v}^2\right)}{{\left(a^2{u}^2+b^2{v}^2\right)}^{3/2}}\mathit{\text{dudv}}$$E36

$$\mathrm{\Delta }\eta =\frac{-4{\mathit{ab}}^2}{2\mathit{\pi \gamma }}{\int }_0^1{\int }_0^1\frac{v^2\left({\alpha }_{01}+{\alpha }_{21}{v}^2\right)}{{\left(a^2{u}^2+b^2{v}^2\right)}^{3/2}}\mathit{\text{dudv}}$$E37

$$\mathrm{\Delta }L=\frac{4\mathit{ab}}{2\pi }{\int }_0^1{\int }_0^1\frac{{\alpha }_{20}{u}^2+{\alpha }_{02}{v}^2+{\alpha }_{22}{u}^2{v}^2}{{\left(a^2{u}^2+b^2{v}^2\right)}^{3/2}}\mathit{\text{dudv}}$$E38

Drawing a line from the origin to the upper right corner, it divides the upper right quadrant.

into ${\sigma }_1\in \left[0<u<1,0<v<u\right]$, ${\sigma }_2\in \left[0<u<v,0<u<1\right]$.

The following nonsingular integration transformation is introduced for ${\sigma }_1$

$$\left\{\begin{array}{c}u=v\\ k=\frac{v}{u}\end{array}\right.$$E39

The following nonsingular integration transformation is introduced for ${\sigma }_2$

$$\left\{\begin{array}{c}v=v\\ \lambda =\frac{u}{v}\end{array}\right.$$E40

Inserting$v=\mathit{ku}$(or$u=\mathit{\lambda v}$)into Eqs. (35)–(37),one arrives at

$$\begin{array}{l}\mathrm{\Delta }\zeta =\frac{4\mathit{ab}}{2\mathit{\pi \gamma }}{\int }_0^1\left({\alpha }_{00}+\frac{1}{3}{\alpha }_{20}+\frac{1}{3}{\alpha }_{02}{k}^2+\frac{1}{5}{\alpha }_{22}{k}^2\right)\frac{\mathit{dk}}{\sqrt{a^2+b^2{k}^2}}\\ +\frac{4\mathit{ab}}{2\mathit{\pi \gamma }}{\int }_0^1\left({\alpha }_{00}+\frac{1}{3}{\alpha }_{02}+\frac{1}{3}{\alpha }_{20}{\lambda }^2+\frac{1}{5}{\alpha }_{22}{\lambda }^2\right)\frac{\mathit{d\lambda }}{\sqrt{b^2+a^2{\lambda }^2}}\end{array}$$E41

$$\begin{array}{l}\mathrm{\Delta }\xi =\frac{-4a^{2}b}{2\mathit{\pi \gamma }}{\int }_0^1\left({\alpha }_{10}+\frac{1}{3}{\alpha }_{12}{k}^2\right)\frac{\mathit{dk}}{{\left(a^2+b^2{k}^2\right)}^{3/2}}\\ -\frac{4a^{2}b}{2\mathit{\pi \gamma }}{\int }_0^1\left({\alpha }_{10}+\frac{1}{3}{\alpha }_{12}\right)\frac{{\lambda }^2\mathit{d\lambda }}{{\left(b^2+a^2{\lambda }^2\right)}^{3/2}}\end{array}$$E42

$$\begin{array}{l}\mathrm{\Delta }\eta =\frac{-4{\mathit{ab}}^2}{2\mathit{\pi \gamma }}{\int }_0^1\left({\alpha }_{01}+\frac{1}{3}{\alpha }_{21}\right)\frac{k^2\mathit{dk}}{{\left(a^2+b^2{k}^2\right)}^{3/2}}\\ -\frac{4{\mathit{ab}}^2}{2\mathit{\pi \gamma }}{\int }_0^1\left({\alpha }_{01}+\frac{1}{3}{\alpha }_{21}{\lambda }^2\right)\frac{\mathit{d\lambda }}{{\left(b^2+a^2{\lambda }^2\right)}^{3/2}}\end{array}$$E43

$$\begin{array}{l}\mathrm{\Delta }L=\frac{4\mathit{ab}}{2\pi }{\int }_0^1\left({\alpha }_{20}+{\alpha }_{02}{k}^2+\frac{1}{3}{\alpha }_{22}{k}^2\right)\frac{\mathit{dk}}{{\left(a^2+b^2{k}^2\right)}^{3/2}}\\ +\frac{4\mathit{ab}}{2\pi }{\int }_0^1\left({\alpha }_{00}+{\alpha }_{20}{\lambda }^2+\frac{1}{3}{\alpha }_{22}{\lambda }^2\right)\frac{\mathit{d\lambda }}{{\left(b^2+a^2{\lambda }^2\right)}^{3/2}}\end{array}$$E44

Now we can see that the denominators are greater than zero after transformation Eq. (39) and Eq. (40), and the singularities have been eliminated. The integrals in x and y directions are converted to the integrals of the powers of k and$\lambda$. This basically changes the double integrals to single variable integrals, which could easily be calculated in Mathematica as follows:

E45

E46

E47

E48

In case of a square grid with a unit length, Eqs. (45)–(48) can be simplified in Mathematica as.

E49

E50

E51

E52

One could find that it is greatly fussy to complete these integrals by hand, which could be easily realized using some commands of computer algebra system.

4. The series expansions of direct transformations between three anomalies in satellite geodesy

The determination of satellite orbit is one of the fundamental problems in satellite geodesy. A graphical representation of the Keplerian orbit is given in Figure 2, see [20].

Figure 2.

Keplerian orbit.

Eccentric, mean and true anomalies are used to describe the movement of satellites. Their transformations are often to be dealt with in satellite ephemeris computation and orbit determination of the spacecraft. In Figure 2, E is the Eccentric anomaly, $\upsilon$is the true anomaly. In order to realize the direct transformations between these anomalies, the series expansions of their transformations are derived using the power series method with the help of computer algebra system Mathematica. Their coefficients are expressed in a power series of the orbital eccentricity $e$ and extended up to eighth-order terms of the orbital eccentricity.

4.1. The series expansions of the direct transformation between eccentric and mean anomalies

Let the mean anomaly be $M$. $M$can be expressed by $E$ as follows:

$$M=E-esinE$$E53

Differentiating the both sides of Eq. (53) yields

$$\frac{\mathit{dE}}{\mathit{dM}}=\frac{1}{1-ecosE}$$E54

To expand Eq. (54) into a power series of $cosM$, we introduce the following new variable

$$t=cosM$$E55

therefore

$$\frac{\mathit{dM}}{\mathit{dt}}=-\frac{1}{sinM}$$E56

and then denote

$$f\left(t\right)=\frac{\mathit{dE}}{\mathit{dM}}=\frac{1}{1-ecosE}$$E57

Substituting $E_0=\frac{\pi }{2}$ into Eq. (53) yields

$$M_0=\frac{\pi }{2}-e$$E58

Substituting Eq. (59) into Eq. (55), one arrives at

$$t_0=sine$$E59

Making use of the chain rule of implicit differentiation

$$\begin{array}{l}{f}^{\prime }\left(t\right)=\frac{\mathit{df}}{\mathit{dE}}\frac{\mathit{dE}}{\mathit{dM}}\frac{\mathit{dM}}{\mathit{dt}}\\ f^{″}\left(t\right)=\frac{d{f}^{\prime }}{\mathit{dE}}\frac{\mathit{dE}}{\mathit{dM}}\frac{\mathit{dM}}{\mathit{dt}}+\frac{d{f}^{\prime }}{\mathit{dM}}\frac{\mathit{dM}}{\mathit{dt}}\\ \cdots \end{array}$$

It is easy to expand Eq. (58) into a power series of $t_0$

$$f\left(t\right)=\frac{\mathit{dE}}{\mathit{dM}}=f\left(t_0\right)+f^{\prime }\left(t_0\right)\left(t-t_0\right)+\frac{1}{2!}{f}^{″}\left(t_0\right){\left(t-t_0\right)}^2+\frac{1}{3!}{f}^{‴}\left(t_0\right){\left(t-t_0\right)}^3+\cdots$$E60

One can imagine that these procedures are too complicated to be realized by hand, but will become much easier and be significantly simplified by means of Mathematica. Omitting the detailed procedure in Mathematica, one arrives at

$$\begin{array}{l}\frac{\mathit{dE}}{\mathit{dM}}=1+b_1\left(cosM-sine\right)+\frac{b_2}{2!}{\left(cosM-sine\right)}^2+\frac{b_3}{3!}{\left(cosM-sine\right)}^3+\frac{b_4}{4!}{\left(cosM-sine\right)}^4\\ \begin{array}{ccc}& & +\frac{b_5}{5!}\end{array}{\left(cosM-sine\right)}^5+\frac{b_6}{6!}{\left(cosM-sine\right)}^6+\frac{b_7}{7!}{\left(cosM-sine\right)}^7+\frac{b_8}{8!}{\left(cosM-sine\right)}^8\end{array}$$E61

where

$$\left\{\begin{array}{l}{b}_1=e+\frac{1}{2}{e}^3+\frac{1}{24}e+5\frac{61}{720}{e}^7\\ b_2=4e^2+\frac{13}{3}{e}^4+\frac{47}{15}{e}^6+\frac{121}{63}{e}^8\\ b_3=27e^3+\frac{91}{2}{e}^5+\frac{1127}{24}{e}^7\\ b_4=256e^4+\frac{2937}{5}{e}^6+\frac{82771}{105}{e}^8\\ b_5=3125e^5+\frac{18173}{2}{e}^7\\ b_6=46656e^6+\frac{1150593}{7}{e}^8\\ b_7=823543e^7\\ b_8=16777216e^8\end{array}\right.$$E62

Integrating at the both sides of Eq. (62) gives the series expansion

$$\begin{array}{l}E=M+{\alpha }_{1}sinM+{\alpha }_{2}sin2M+{\alpha }_{3}sin3M+{\alpha }_{4}sin4M+{\alpha }_{5}sin5M\\ \begin{array}{cc}& \begin{array}{c}\end{array}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}+{\alpha }_{6}sin6M+{\alpha }_{7}sin7M+{\alpha }_{8}sin8M\end{array}$$E63

where

$$\left\{\begin{array}{l}{\alpha }_1=e-\frac{1}{8}{e}^3+\frac{1}{192}{e}^5-\frac{1}{9216}{e}^7\\ {\alpha }_2=\frac{1}{2}{e}^2-\frac{1}{6}{e}^4+\frac{1}{48}{e}^6-\frac{1}{720}{e}^8\\ {\alpha }_3=\frac{3}{8}{e}^3-\frac{27}{128}{e}^5+\frac{243}{5120}{e}^7\\ {\alpha }_4=\frac{1}{3}{e}^4-\frac{4}{15}{e}^6+\frac{4}{45}{e}^8\\ {\alpha }_5=\frac{125}{384}{e}^5-\frac{3125}{9216}{e}^7\\ {\alpha }_6=\frac{27}{80}{e}^6-\frac{243}{560}{e}^8\\ {\alpha }_7=\frac{16807}{46080}{e}^7\\ {\alpha }_8=\frac{128}{315}{e}^8\end{array}\right.$$E64

4.2. The series expansions of the direct transformation between eccentric and true anomalies

The true anomaly $\upsilon$ can be expressed by E as follows:

$$tan\frac{\upsilon }{2}=\sqrt{\frac{1+e}{1-e}}tan\frac{E}{2}$$E65

Therefore, it holds

$$\upsilon =2arctan\left(\sqrt{\frac{1+e}{1-e}}tan\frac{E}{2}\right)$$E66

One could expand $\upsilon$ as a power series of the eccentricity at $e=0$ in order to obtain the direct series expansion of the transformation from $E$ to $\upsilon$. Omitting the detailed procedure in Mathematica, one arrives at

$$\begin{array}{l}\upsilon =E+{\beta }_{1}sinE+{\beta }_{2}sin2E+{\beta }_{3}sin3E+{\beta }_{4}sin4E+{\beta }_{5}sin5E\\ \begin{array}{cc}& \begin{array}{c}\end{array}\begin{array}{c}\end{array}\end{array}+{\beta }_{6}sin6E+{\beta }_{7}sin7E+{\beta }_{8}sin8E\end{array}$$E67

where

$$\left\{\begin{array}{l}{\beta }_1=e+\frac{1}{4}{e}^3+\frac{1}{8}{e}^5+\frac{5}{64}{e}^7\\ {\beta }_2=\frac{1}{4}{e}^2+\frac{1}{8}{e}^4+\frac{5}{64}{e}^6+\frac{7}{128}{e}^8\\ {\beta }_3=\frac{1}{12}{e}^3+\frac{1}{16}{e}^5+\frac{3}{64}{e}^7\\ {\beta }_4=\frac{1}{32}{e}^4+\frac{1}{32}{e}^6+\frac{7}{256}{e}^8\\ {\beta }_5=\frac{1}{80}{e}^5+\frac{1}{64}{e}^7\\ {\beta }_6=\frac{1}{192}{e}^6+\frac{1}{128}{e}^8\\ {\beta }_7=\frac{1}{448}{e}^7\\ {\beta }_8=\frac{1}{1024}{e}^8\end{array}\right.$$E68

From Eq. (67), one knows

$$E=2arctan\left(\sqrt{\frac{1-e}{1+e}}tan\frac{\upsilon }{2}\right)$$E69

Expanding $E$ as a power series of the eccentricity at $e=0$ by means of Mathematica yields the direct series expansion of the transformation from $\upsilon$ to $E$

$$\begin{array}{l}E=\upsilon +{\gamma }_{1}sin\upsilon +{\gamma }_{2}sin2\upsilon +{\gamma }_{3}sin3\upsilon +{\gamma }_{4}sin4\upsilon +{\gamma }_{5}sin5\upsilon \\ \begin{array}{cc}\begin{array}{cc}& \end{array}& +\end{array}{\gamma }_{6}sin6\upsilon +{\gamma }_{7}sin7\upsilon +{\gamma }_{8}sin8\upsilon \end{array}$$E70

where

$$\left\{\begin{array}{l}{\gamma }_1=-e-\frac{1}{4}{e}^3-\frac{1}{8}{e}^5-\frac{5}{64}{e}^7\\ {\gamma }_2=\frac{1}{4}{e}^2+\frac{1}{8}{e}^4+\frac{5}{64}{e}^6+\frac{7}{128}{e}^8\\ {\gamma }_3=-\frac{1}{12}{e}^3-\frac{1}{16}{e}^5-\frac{3}{64}{e}^7\\ {\gamma }_4=\frac{1}{32}{e}^4+\frac{1}{32}{e}^6+\frac{7}{256}{e}^8\\ {\gamma }_5=-\frac{1}{80}{e}^5-\frac{1}{64}{e}^7\\ {\gamma }_6=\frac{1}{192}{e}^6+\frac{1}{128}{e}^8\\ {\gamma }_7=-\frac{1}{448}{e}^7\\ {\gamma }_8=\frac{1}{1024}{e}^8\end{array}\right.$$E71

4.3. The series expansions of the direct transformation between mean and true anomalies

The whole formulae for the transformation from $M$ to $\upsilon$ are as follows

$$\left\{\begin{array}{l}E=M+{\alpha }_{1}sinM+{\alpha }_{2}sin2M+{\alpha }_{3}sin3M+{\alpha }_{4}sin4M+{\alpha }_{5}sin5M\\ \begin{array}{cc}& \begin{array}{c}\end{array}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}+{\alpha }_{6}sin6M+{\alpha }_{7}sin7M+{\alpha }_{8}sin8M\\ \upsilon =E+{\beta }_{1}sinE+{\beta }_{2}sin2E+{\beta }_{3}sin3E+{\beta }_{4}sin4E+{\beta }_{5}sin5E\\ \begin{array}{cc}& \begin{array}{c}\end{array}\begin{array}{c}\end{array}\end{array}+{\beta }_{6}sin6E+{\beta }_{7}sin7E+{\beta }_{8}sin8E\end{array}\right.$$E72

Since the coefficients ${\alpha }_i$,${\beta }_i$ ($i=1,2,\cdots 8$) are expressed in a power series of the eccentricity, one could expand $\upsilon$ as a power series of the eccentricity at $e=0$ in order to obtain the direct expansion of the transformation from $M$ to $\upsilon$. Omitting the main operations by means of Mathematica, one arrives at the direct expansion of the transformation from $M$ to $\upsilon$

$$\begin{array}{l}\upsilon =M+{\delta }_{1}sinM+{\delta }_{2}sin2M+{\delta }_{3}sin3M+{\delta }_{4}sin4M+{\delta }_{5}sin5M\\ \begin{array}{cc}\begin{array}{cc}& \end{array}& +{\delta }_6\end{array}sin6M+{\delta }_{7}sin7M+{\delta }_{8}sin8M\end{array}$$E73

where

$$\left\{\begin{array}{l}{\delta }_1=2e-\frac{1}{4}{e}^3+\frac{5}{96}{e}^5+\frac{107}{4608}{e}^7\\ {\delta }_2=\frac{5}{4}{e}^2-\frac{11}{24}{e}^4+\frac{17}{192}{e}^6+\frac{43}{5760}{e}^8\\ {\delta }_3=\frac{13}{12}{e}^3-\frac{43}{64}{e}^5+\frac{95}{512}{e}^7\\ {\delta }_4=\frac{103}{96}{e}^4-\frac{451}{480}{e}^6+\frac{4123}{11520}{e}^8\\ {\delta }_5=\frac{1097}{960}{e}^5-\frac{5957}{4608}{e}^7\\ {\delta }_6=\frac{1223}{960}{e}^6-\frac{7913}{4480}{e}^8\\ {\delta }_7=\frac{47273}{32256}{e}^7\\ {\delta }_8=\frac{556403}{322560}{e}^8\end{array}\right.$$E74

The whole formulae for the transformation from $\upsilon$ to $M$are as follows

$$\left\{\begin{array}{l}E=\upsilon +{\gamma }_{1}sin\upsilon +{\gamma }_{2}sin2\upsilon +{\gamma }_{3}sin3\upsilon +{\gamma }_{4}sin4\upsilon +{\gamma }_{5}sin5\upsilon \\ \begin{array}{cc}\begin{array}{cc}& \end{array}& +\end{array}{\gamma }_{6}sin6\upsilon +{\gamma }_{7}sin7\upsilon +{\gamma }_{8}sin8\upsilon \\ M=E-esinE\end{array}\right.$$E75

Expanding $M$as a power series of the eccentricity at $e=0$ by means of Mathematica yields the direct series expansion of the transformation from $\upsilon$ to $M$

$$\begin{array}{l}M=\upsilon +{\epsilon }_{1}sin\upsilon +{\epsilon }_{2}sin2\upsilon +{\epsilon }_{3}sin3\upsilon +{\epsilon }_{4}sin4\upsilon \\ \begin{array}{cc}\begin{array}{cc}& \end{array}& +{\epsilon }_{5}sin5\upsilon +\end{array}{\epsilon }_{6}sin6\upsilon +{\epsilon }_{7}sin7\upsilon +{\epsilon }_{8}sin8\upsilon \end{array}$$E76

where

$$\left\{\begin{array}{l}{\epsilon }_1=-2e\\ {\epsilon }_2=\frac{3}{4}{e}^2+\frac{1}{8}{e}^4+\frac{3}{64}{e}^6+\frac{3}{128}{e}^8\\ {\epsilon }_3=-\frac{1}{3}{e}^3-\frac{1}{8}{e}^5-\frac{7}{16}{e}^7\\ {\epsilon }_4=\frac{5}{32}{e}^4+\frac{3}{32}{e}^6+\frac{15}{256}{e}^6\\ {\epsilon }_5=-\frac{3}{40}{e}^5-\frac{1}{16}{e}^7\\ {\epsilon }_6=\frac{7}{192}{e}^6+\frac{5}{128}{e}^6\\ {\epsilon }_7=-\frac{1}{56}{e}^7\\ {\epsilon }_8=\frac{9}{1024}{e}^8\end{array}\right.$$E77

4.4. The accuracy of the derived series expansions

In order to validate the exactness of the derived series expansions, one has examined their accuracies when the orbital eccentricity e is respectively equal to 0.01, 0.05, 0.1 and 0.2.

One makes use of Eq. (53) and Eq. (67) to obtain the theoretical value M₀ and ${\upsilon }_0$ at given geodetic latitude E₀. Substituting M₀ into Eq. (64) and Eq. (74), one arrives at the computation value E₁ and ${\upsilon }_1$. Substituting ${\upsilon }_0$ into Eq. (71) and Eq. (77), one arrives at the computation value E₂ and M₁. Substituting E₀ into Eq. (68), one arrives at the computation value ${\upsilon }_2$. The differences between the computation and theoretical values indicate the accuracies of the derived series expansions, which are denoted as $\mathrm{\Delta }{E}_1$_, $\mathrm{\Delta }{\upsilon }_1/\left({}^{″}\right)$_, $\mathrm{\Delta }{M}_1/\left({}^{″}\right)$_,$\mathrm{\Delta }{E}_2/\left({}^{″}\right)$_,$\mathrm{\Delta }{\upsilon }_2/\left({}^{″}\right)$. Due to limited space, these errors when e is equal to 0.05 are only listed in Table 2.

$E_0/\left({}^{\circ }\right)$	20	40	60	80
$\mathrm{\Delta }{E}_1/\left({}^{″}\right)$	$8.2\times {10}^{-8}$	$-1.7\times {10}^{-7}$	1.6 × 10^−7	$2.2\times {10}^{-8}$
$\mathrm{\Delta }{\upsilon }_1/\left({}^{″}\right)$	$3.5\times {10}^{-7}$	$-7.4\times {10}^{-7}$	$6.8\times {10}^{-7}$	$9.9\times {10}^{-8}$
$\mathrm{\Delta }{M}_1/\left({}^{″}\right)$	3.5 × 10^−8	$-5.2\times {10}^{-9}$	$-7.1\times {10}^{-9}$	$-3.7\times {10}^{-9}$
$\mathrm{\Delta }{E}_2/\left({}^{″}\right)$	$2.8\times {10}^{-8}$	$2.1\times {10}^{-8}$	$1.4\times {10}^{-8}$	$1.2\times {10}^{-8}$
$\mathrm{\Delta }{\upsilon }_2/\left({}^{″}\right)$	$-2.9\times {10}^{-8}$	$-2.4\times {10}^{-8}$	$-1.5\times {10}^{-8}$	$-1.3\times {10}^{-8}$

Table 2.

Errors of the derived series expansions.

From Table 2, one could find that the accuracy of derived series expansions is higher than ${10}^{-5\prime \prime }$, which could satisfy practical application. Other numerical examples indicate that when the orbital eccentricity $e$ is respectively equal to 0.01, 0.01 and 0.2, the accuracy of derived series expansions is correspondingly higher than ${10}^{-10\prime \prime }$, ${10}^{-3\prime \prime }$ and ${0.1}^{\prime \prime }$.

5. Conclusions

Some typical mathematical problems in geodesy are solved by means of computer algebra analysis method and computer algebra system Mathematica. The main contents and research results presented in this chapter are as follows:

1. The forward and inverse expansions of the meridian arc often used in geometric geodesy are derived. Their coefficients are expressed in a power series of the first eccentricity of the reference ellipsoid and extended up to its tenth-order terms.

2. The singularity existing in the integral expressions of height anomaly, deflections of the vertical and gravity gradient is eliminated using the nonsingular integration transformations, and the nonsingular expressions are systematically derived.

3. The series expansions of direct transformations between three anomalies in satellite geodesy are derived using the power series method. Their coefficients are expressed in a power series of the orbital eccentricity e and extended up to eighth-order terms of the orbital eccentricity.

Acknowledgments

This work was financially supported by National Natural Science Foundation of China (Nos. 41631072, 41771487, 41571441).

Conflict of interest

There are no conflicts of interest.

Thanks

We would like to express our great appreciation to the editor and reviewers. Thanks very much for the kind work and consideration of Ms. Kristina Jurdana on the publication of this chapter.