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

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 expan- sions 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.


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.

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.

The forward expansion of the meridian arc
The meridian arc from the equator where the latitude is from where X is the meridian arc; Bis 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 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].

The inverse expansion of the meridian arc using the Hermite interpolation method
Differentiation to the both sides of Eq. (1) yields: Define ψ as Substituting Eq. (5) into Eq. (4) yields: Suppose that the inverse solution of Eq. (6) has the following form: B ¼ ψ þ a 2 sin 2ψ þ a 4 sin 4ψ þ a 6 sin 6ψ þ a 8 sin 8ψ þ a 10 sin 10ψ 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 The further derivation of Eq. (6) with respect to ψ yields B 00 ψ ð Þ. Unfortunately, B 00 ψ ð Þ is equal to Hence, one differentiates Eq. (6) twice and it yields B ‴ ψ ð Þ. Omitting the derivative procedure by means of Mathematica, one arrives at The solution to Eq. (13) is À64 À216 À512 1000 32 1024 7776 32768 100000 Omitting the main operations by means of Mathematica, one arrives at

The inverse expansion of the meridian arc using Lagrange's theorem method
Suppose that with f x ð Þ j j<< x j j and y ≈ x. Lagrange's theorem states that in a suitable domain the solution of Eq. (16) is Supposef x ð Þis defined by the following finite trigonometric series: 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 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 where Substituting x for Band y for ψ in Eq. (16), the coefficients α, β, γ, δ, ε are consistent with α 2 , α 4 , α 6 , α 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 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 The results in Eq. (24) are consistent with the coefficients in Eq. (15), which substantiates the correctness of the derived formula.

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 parametersa ¼ 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 ψ 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 ψ 0 into Eq. (22), one arrives at the computation value B 0 1 .
Þ indicate the accuracies of the inverse expansions derived by Yang (see [6,7]) and the author respectively. These errors are listed in Table 1.
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 À500 , and the accuracy of the inverse expansion Eq. (22) derived by the author is higher than 10 À700 . The accuracy is improved by 2 orders of magnitude by means of computer algebra.

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 where r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , Δg 0 is the gravity anomaly at the computation point. When r ! 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 0 Â 2 0 ), 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 σ ∈ Àa < x < a; Àb < y < b ½ (a > 0, b > 0) due to the convergence of meridian, and the gravity anomaly is expressed as double quadratic polynomial: Inserting the innermost area into Eqs. (25)-(28), and let the contributions be Δζ, Δξ, Δηand ΔL, one arrives at The following transformation is introduced for σ Using the properties of the integration for even/odd functions and exploiting the symmetry of the integration area, one arrives at Δη ¼ À4ab 2 2πγ ΔL ¼ 4ab 2π Drawing a line from the origin to the upper right corner, it divides the upper right quadrant.
The following nonsingular integration transformation is introduced for σ 1 The following nonsingular integration transformation is introduced for Insertingv ¼ ku(oru ¼ λv)into Eqs. (35)-(37),one arrives at 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λ. This basically changes the double integrals to single variable integrals, which could easily be calculated in Mathematica as follows: 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.

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].
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, υ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.

The series expansions of the direct transformation between eccentric and mean anomalies
Let the mean anomaly be M. Mcan be expressed by E as follows: Differentiating the both sides of Eq. (53) yields Integrating at the both sides of Eq. (62) gives the series expansion

The series expansions of the direct transformation between eccentric and true anomalies
The true anomaly υ can be expressed by E as follows: One could expand υ as a power series of the eccentricity at e ¼ 0 in order to obtain the direct series expansion of the transformation from E to υ. Omitting the detailed procedure in Mathematica, one arrives at From Eq. (67), one knows 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 υ to E E ¼ υ þ γ 1 sin υ þ γ 2 sin 2υ þ γ 3 sin 3υ þ γ 4 sin 4υ þ γ 5 sin 5υ þ γ 6 sin 6υ þ γ 7 sin 7υ þ γ 8

The series expansions of the direct transformation between mean and true anomalies
The whole formulae for the transformation from M to υ are as follows Since the coefficients α i ,β i (i ¼ 1, 2, ⋯8) are expressed in a power series of the eccentricity, one could expand υ as a power series of the eccentricity at e ¼ 0 in order to obtain the direct expansion of the transformation from M to υ. Omitting the main operations by means of Mathematica, one arrives at the direct expansion of the transformation from M to υ υ ¼ M þ δ 1 sin M þ δ 2 sin 2M þ δ 3 sin 3M þ δ 4 sin 4M þ δ 5 sin 5M The whole formulae for the transformation from υ to Mare as follows E ¼ υ þ γ 1 sin υ þ γ 2 sin 2υ þ γ 3 sin 3υ þ γ 4 sin 4υ þ γ 5 sin 5υ þ γ 6 sin 6υ þ γ 7 sin 7υ þ γ 8 sin 8υ Expanding Mas a power series of the eccentricity at e ¼ 0 by means of Mathematica yields the direct series expansion of the transformation from υ to M M ¼ υ þ ε 1 sin υ þ ε 2 sin 2υ þ ε 3 sin 3υ þ ε 4 sin 4υ þε 5 sin 5υþ ε 6 sin 6υ þ ε 7 sin 7υ þ ε 8 sin 8υ (76) where ε 1 ¼ À2e  68), one arrives at the computation value υ 2 . The differences between the computation and theoretical values indicate the accuracies of the derived series expansions, which are denoted as ΔE 1, Δυ 1 = 00 ð Þ , ΔM 1 = 00 ð Þ , ΔE 2 = 00 ð Þ , Δυ 2 = 00 ð Þ. Due to limited space, these errors when e is equal to 0.05 are only listed in Table 2.
From Table 2, one could find that the accuracy of derived series expansions is higher than 10 À5 00 , 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 00 , 10 À3 00 and 0:1 00 .

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.