A Physics-Informed Neural Network for Geodesic Computation and Visualization of Geodesic Spheres on Riemannian Manifolds

Introduction

Statement of the problem

Imagination comes when one visualize something and then extends it in mind. Visualization plays a very important role in mathematics and specially in geometry. Objects that appear in Euclidean geometry are easy to visualize. Any function $f:\mathbb{R}\to\mathbb{R}$ can be plotted as a subset $G(f)=\{(x,f(x)):x\in\mathbb{R}\}\subset\mathbb{R}^{2}$. In similar fashion a map $f:\mathbb{R}^{2}\to\mathbb{R}$ can be visualized using the graph $G(f)=\{(x,y,f(x,y)):(x,y)\in\mathbb{R}^{2}\}\subset\mathbb{R}^{3}$. Things gets interesting when dimension reaches to 4. For example, a map $f:\mathbb{R}^{3}\to\mathbb{R}$ or $f:\mathbb{R}^{2}\to\mathbb{R}^{2}$ can produce a graph $G(f)$ as a subset of $\mathbb{R}^{4}$. We have no direct access to the 4 dimensional space and hence they can not be visualized in the traditional manner. Another obstacle is the metric defined on the space. For example, we are able to plot those graphs on $\mathbb{R}^{2}$ or $\mathbb{R}^{3}$ under the Euclidean metric. Now what if the metric changes? Can we visualize them on the space again? One possible answer is that if $(\mathbb{R},d_{u})$ is the usual metric space then the space $(\mathbb{R},\lambda d_{u})$ is the scaled metric space whose open unit spheres are going to be $\lambda$ times bigger/smaller than the original space. But this explanation is going to work simple looking metrics. We are talking about visualizing open unit spheres of Riemannian manifolds. One may say that yes we can see it; the sphere $S^{2}$, the cube $[0,1]^{3}$, we can draw them and it has its own Riemannian metric. But my question is what if we have a Riemannian metric, for example $ds^{2}=x^{2}dx^{2}+\sin y\,dy^{2}$, how to “see” the open unit sphere under this metric?

In this article we are going to visualize the open unit spheres of a Riemannian manifold equipped with an abstract¹ Riemannian metric $g$ by solving the geodesic equations. We are also going to explore a volume comparison between the open sphere of the manifold and the open Euclidean sphere. The geodesic equations that are arising from the metric can be solved via Python-SciPy’s integrate IVP module, but it will produce a discrete solution, however if we configure SciPy then it will interpolate between the points and provide a continuous solution. But we create a Physics-Informed Neural Network (PINN) model which will solve the equations in a continuous manner. We are applying Neural Network to understand the nature of the solution in a continuous manner rather than having an interpolation of noisy data. This implementation helps us to study real world data and it also shows that PINN can be applied to solve geodesic equations over Riemannian manifolds. We also compare both the solution (PINN vs SciPy) to verify the accuracy and stability of the model.

Related works

This kind of work is not new in computational Riemannian geometry. In the past these works were done manually i.e., by hand. One can see the work of Hui and his coauthors [5, 7, 6, 8] for such investigations. But after the involvement of mathematical softwares like Mathematica, Matlab and Python especially after the introduction of symbolic computations several researchers have been investigating Riemannian manifolds and its various curvature tensors. One can see the work of Shaikh et al. [13, 12] for analysis on space-time metrics. For computations over vector field related equations we refer to the works of Blaga [4, 2, 3]. On the other hand several nonlinear differential equations were being solved using PINN, in this regard the works of Sur [15, 14] is worth mentioning. Other research leads to theoretical understanding of the space where as this article provides a direct visualization of the unseen open spheres of the Riemannian space. To the best of the authors knowledge and belief this kind of study is new to computational Riemannian geometry that brings a direct visualization of the open spheres of abstract Riemannian metrics.

Preliminaries

We start with recalling some important notions that have been used in this article. For detailed study we recommend the book of J. M. Lee [9] for Riemannian geometry and Geodesics in Riemannian manifolds, for PINN we recommend Bai [1] and Liu [10].

Riemannian manifold and Riemannian metric

An $n$-dimensional Riemannian manifold is an $n$-dimensional smooth manifold $M$ with a globally defined 2-form $g:\chi(M) \times\chi(M) \to C^\infty(M)$ which is

1. positive definite i.e., $g(X,X) \ge0$, $\forall X \in\chi(M)$ and $g(X,X) = 0 \Longleftrightarrow X = 0$.

2. symmetric i.e., $g(X,Y) = g(Y,X)$, $\forall X,Y \in\chi(M)$.

3. $C^\infty(M)$-bilinear i.e., $g(fX,Y) = g(X,fY) = fg(X,Y)$, $\forall f \in C^\infty(M)$, $\forall X,Y \in\chi(M)$.

This 2-form is actually an inner product in the tangent spaces $T_p(M)$ for each $p \in M$. Let $p,q \in M$ and let $\gamma:[0,1] \to M$ be a smooth curve, then the length of the curve is given by

This length functional is independent of parametrization. So we can measure the distance between any two points $p,q$ on the manifold $M$ by the metric induced by the length structure given by

As the distance metric $d_\ell$ depends solely on the coefficients of $g$, so $g$ is called the Riemannian metric. If one fixes an orthonormal basis $\{e_i:i=1,2,\ldots,n\}$ for $\chi(M)$ then we have a local expression (also known as local form) of $g$ given by

These $g_{ij}$ coefficients are real valued functions of $n$ variables $x_1,x_2,\cdots,x_n$ and it forms a symmetric matrix which we denote as $g$ (the same notation) and we write

For our numerical computation we will be using this form to consider our matrix of the Riemannian metric. For any two points $p,q\in M$, a curve $\gamma:[0,1]\to M$ joining them is said to be the shortest path if $d_{\ell}(p,q)=\ell(\gamma)$. In the next section we will be discussing geodesics and we will show how the shortest paths becomes the geodesics.

1.3.2. Geodesics in Riemannian manifolds. A Riemannian metric contains important information about the geometric structure of the underlying manifold, particularly its curvature. The curvature of a Riemannian manifold is determined by its Levi-Civita connection, whose local representation is given by the Christoffel symbols. These symbols are computed from the metric tensor and its first-order partial derivatives. The Christoffel symbols associated with a Riemannian metric $g$ are defined in local coordinates by

where $g^{il}$ denotes the components of the inverse metric tensor. The Christoffel symbols describe the local behavior of the Levi-Civita connection and are used to express covariant derivatives in local coordinates. In particular, they determine the Riemann curvature tensor, which measures the failure of covariant derivatives to commute.

The Riemann curvature tensor is defined by

where $\chi(M)$ denotes the space of smooth vector fields on $M$, $\nabla$ is the Levi-Civita connection, and $[X,Y]$ is the Lie bracket of vector fields. The Riemann curvature tensor is a tensor field of type $(1,3)$, and its local components, under the convention $R(\partial_k,\partial_l)\partial_j=R^{i}_{jkl}\partial_i$, are given by

Contracting the first and third tensor indices yields the Ricci curvature tensor. Its local expression is

where the repeated index $i$ is summed over as per Einstein’s summation convention. A further contraction with respect to the indices $j$ and $k$ gives the scalar curvature:

These local expressions are useful for computing the components of the Christoffel symbols, the Riemann curvature tensor, the Ricci tensor, and the scalar curvature symbolically.

The Christoffel symbols also play a fundamental role in determining the geodesics of a Riemannian manifold. Geodesics are curves whose tangent vectors are parallel transported along the curves themselves. Let $\gamma(t) = (x^1(t), \ldots, x^n(t))$ be a smooth curve on $M$, with tangent vector

A curve is a geodesic if its covariant acceleration vanishes, that is,

Expressing this condition in local coordinates and using the Christoffel symbols gives the geodesic equations

Thus, the Christoffel symbols determine how the coordinate components of the velocity of a geodesic change along the curve. Since the Christoffel symbols are determined by the metric and its first-order derivatives, the geodesic equations describe how the metric governs the local behavior of geodesics. Consequently, both curvature computations and geodesic calculations can be carried out symbolically once the metric tensor is specified.

Physics-Informed Neural Network

Physics-Informed Neural Network is a particular type of neural network from Deep Learning which helps us to solve mathematical problems by incorporating known physical laws, such as differential equations, into the training process. In ordinary machine learning, a neural network generally learns from a given set of data. In a PINN, the network learns from data and from the governing equations of the problem. Some times we often let the system learn from the governing equations only to produce a data independent model and then compare the model’s solution to known experimental data for accuracy. This makes PINNs particularly useful for solving ordinary differential equations (ODEs), partial differential equations (PDEs), and certain inverse problems. The framework was popularized by Raissi, Perdikaris, and Karniadakis in their 2019 paper [11].

Mathematical formulation of PINN:

Consider a general differential equation of the form

where $\Omega\subset\mathbb{R}^n$ denotes the computational domain, $u$ is the unknown function, $\mathcal{D}$ is a differential operator, and $f$ is a given source term.

A Physics-Informed Neural Network (PINN) approximates the unknown solution $u$ by a neural network $u_\theta$, where $\theta$ represents the trainable parameters of the network. Thus, we write

To incorporate the governing differential equation into the training process, we define the physics residual as

If the neural network satisfies the differential equation exactly, then $p_\theta(x)=0$ throughout the domain $\Omega$. In practice, the residual is evaluated at a finite collection of collocation points $\{x_i\}_{i=1}^{N_p}\subset\Omega$, where $N_p$ denotes the number of collocation points. The corresponding physics loss is defined by

Suppose that the differential equation is subject to the boundary condition

where $\partial\Omega$ denotes the boundary of the domain and $\phi$ is a prescribed boundary function. The boundary condition is incorporated into the training process through the boundary loss

where $\{x_i^b\}_{i=1}^{N_b}$ are points on the boundary $\partial\Omega$, and $N_b$ is the number of boundary points.

The total loss function is then constructed as a weighted combination of the physics loss and the boundary loss:

where $\lambda_p$ and $\lambda_b$ are nonnegative weighting parameters that control the relative contributions of the two loss terms. The optimal network parameters are obtained by solving the minimization problem

The minimization is typically performed using gradient-based optimization algorithms. The derivatives of the neural network output with respect to the input coordinates, which are required to evaluate the differential operator $\mathcal{D}$, can be computed efficiently using automatic differentiation.

If observational data are available, an additional data-fitting term can be incorporated into the loss function. For example, given observations $\{(x_i^d,u_i^d)\}_{i=1}^{N_d}$, the data loss can be defined as

In this case, the total loss may be written as

where $\lambda_d$ is the weight assigned to the data loss.

Thus, a PINN combines the governing differential equation, boundary conditions, and, when available, observational data within a unified optimization framework. This approach allows the neural network to approximate solutions of differential equations even when exact solutions or extensive training data are unavailable.

In our case, we are going to implement PINN into solving geodesic equations (11) formed by the metric $g$ in (2) with suitably chosen initial conditions. It should be mentioned here that we are not going to use experimental data for our PINN model, so that the network will learn the system by itself and we will use the SciPy solution only to compare the accuracy of the model’s solution. We specify all the working methodologies in the next section.

Working Method and Configurations

Working method

Step-I (Dimension): First we fix the dimension $n$. For our case it is $n=3$. We are taking a 3 dimensional manifold so that we can have a relatable visualization. Jumping into higher dimension will be easy once we can see what is happening in the lower dimension. One may also consider the same experiment on dimension 2. We also intent to investigate dimension 4 with space-time metrics with possibly singularities in near future.

Step-II (Domain and Coordinate): Next we fix a domain $\Omega\subset\mathbb{R}^{3}$ with coordinate $x,y,z$ in which the Riemannian metric will be defined. In general case the unit cube $[0,1]^{3}$ serves the purpose. But in our case we relax this case by choosing globally defined metric which has no singularity in the whole of $\mathbb{R}^{3}$, thus in our case we specify $\Omega=\{(x,y,z)\in\mathbb{R}^{3}:x>0\}\subset\mathbb{R}^{3}$. However, the parameter $\tau$ that controls the geodesic trajectory is taken in the interval $[0,2]$ so that half length of the geodesic can be traced out in $[0,1]$. The reason behind choosing $x>0$ is specified in the next section.

Step-III (The Riemannian metric): After domain specification we specify the Riemannian metric $g$ in matrix form. In our case we are dealing with a 3 dimensional matrix which is given by

We prefer to start with a simple metric and just changing one element may seem to be not so interesting but the results are very impressive as this will directly show us the effect of changing one parameter in the metric or in the line element. In our case the line element is

We are going to visualize the open spheres under this metric $g$. We keep the metric to be positive definite for now but in future we may extend this to negative definite or indefinite metrics as well.

Step-IV (Study the geometry of $g$): Now we use symbolic calculations to find the expressions for Christoffel symbols using (4), Riemann curvature tensor using (6), Ricci curvature using (7), Scalar curvature using (8).

Step-V (Einstein manifold and Flatness checking) [Optional]: It is not a necessary step but we like to keep it in the algorithm part so that we can report this additional information. A manifold is said to be Einstein if $\operatorname{Ric} = \lambda g$ for some scalar $\lambda$. A manifold is said to be flat if its Riemann curvature vanishes identically, i.e., all of the Riemann curvature components are zero. Similarly, a Ricci flat manifold is the one with vanishing Ricci curvature. So Einstein manifolds can be considered as an extension of Ricci flat manifolds. It has various applications in Astrophysics and Cosmology. We can check for all nonzero terms of Riemann and Ricci tensors after the computation. Depending on their vanishing we will present our verdict. For Einstein manifold checking, we simply compute $\operatorname{Ric} - \lambda g$ and compare it with the $O_{3 \times3}$ (the null matrix), here we solve for $\lambda$ and depending on the number of solutions we will say whether it is Einstein manifold or not.

Step-VI (Geodesic equations): Next using the Christoffel symbols and the coordinates, we form the geodesic equations. We solve them using SciPy’s integrate module. In our case we will have the solution of the form $G = \{(x(\tau), y(\tau), z(\tau)) : \tau\in[0, 2]\}$. This will form the geodesics in a parametric form. Here $x(\tau), y(\tau), z(\tau)$ satisfies the geodesic equations (11) for our specified metric $g$.

Step-VII (PINN): Next we implement PINN to solve the equations and compare the results with SciPy’s outcome. After running the optimizer for 10 different seeds we calculate the overall root mean square error (RMSE) using this formula

where

Please refer to Table 2 for the numerical values.

Table 2. Error metrics for ten independent PINN runs with different random seeds.

Seed	Coordinate	MAE	RMSE	MaxAE	Overall RMSE
1234	$x$	$1.475333 \times10^{-4}$	$1.801110 \times10^{-4}$	$3.399890 \times10^{-4}$	$1.144660 \times10^{-4}$
	$y$	$3.439841 \times10^{-5}$	$4.341200 \times10^{-5}$	$7.044621 \times10^{-5}$
	$z$	$6.393794 \times10^{-5}$	$7.058888 \times10^{-5}$	$1.122085 \times10^{-4}$
1235	$x$	$5.113784 \times10^{-3}$	$5.295075 \times10^{-3}$	$7.063421 \times10^{-3}$	$3.638967 \times10^{-3}$
	$y$	$1.526544 \times10^{-3}$	$1.569376 \times10^{-3}$	$1.985233 \times10^{-3}$
	$z$	$2.936206 \times10^{-3}$	$3.037348 \times10^{-3}$	$3.953934 \times10^{-3}$
1236	$x$	$4.181847 \times10^{-4}$	$4.715432 \times10^{-4}$	$7.168459 \times10^{-4}$	$3.653336 \times10^{-4}$
	$y$	$7.321597 \times10^{-5}$	$9.054659 \times10^{-5}$	$2.007484 \times10^{-4}$
	$z$	$4.063160 \times10^{-4}$	$4.121339 \times10^{-4}$	$4.953508 \times10^{-4}$
1237	$x$	$5.677972 \times10^{-4}$	$6.236564 \times10^{-4}$	$9.560724 \times10^{-4}$	$3.917300 \times10^{-4}$
	$y$	$1.602889 \times10^{-4}$	$1.783923 \times10^{-4}$	$2.461938 \times10^{-4}$
	$z$	$1.728306 \times10^{-4}$	$1.989624 \times10^{-4}$	$3.289656 \times10^{-4}$
1238	$x$	$5.162852 \times10^{-3}$	$5.186230 \times10^{-3}$	$5.921158 \times10^{-3}$	$3.713242 \times10^{-3}$
	$y$	$1.331224 \times10^{-3}$	$1.351524 \times10^{-3}$	$1.674414 \times10^{-3}$
	$z$	$3.552884 \times10^{-3}$	$3.555404 \times10^{-3}$	$3.729424 \times10^{-3}$
1239	$x$	$1.905061 \times10^{-4}$	$2.423644 \times10^{-4}$	$4.548552 \times10^{-4}$	$1.499765 \times10^{-4}$
	$y$	$4.927924 \times10^{-5}$	$5.595751 \times10^{-5}$	$8.532256 \times10^{-5}$
	$z$	$5.685171 \times10^{-5}$	$7.488056 \times10^{-5}$	$1.519464 \times10^{-4}$
1240	$x$	$2.048747 \times10^{-4}$	$2.624104 \times10^{-4}$	$4.853309 \times10^{-4}$	$1.611155 \times10^{-4}$
	$y$	$4.898126 \times10^{-5}$	$5.419302 \times10^{-5}$	$8.166822 \times10^{-5}$
	$z$	$5.842658 \times10^{-5}$	$7.796503 \times10^{-5}$	$1.593409 \times10^{-4}$
1241	$x$	$2.516110 \times10^{-3}$	$2.576858 \times10^{-3}$	$3.237878 \times10^{-3}$	$1.726738 \times10^{-3}$
	$y$	$5.398933 \times10^{-4}$	$5.412438 \times10^{-4}$	$5.860503 \times10^{-4}$
	$z$	$1.397981 \times10^{-3}$	$1.418353 \times10^{-3}$	$1.695860 \times10^{-3}$
1242	$x$	$1.551430 \times10^{-4}$	$1.901201 \times10^{-4}$	$3.664752 \times10^{-4}$	$1.207510 \times10^{-4}$
	$y$	$4.965467 \times10^{-5}$	$5.578220 \times10^{-5}$	$8.464544 \times10^{-5}$
	$z$	$5.297984 \times10^{-5}$	$6.697122 \times10^{-5}$	$1.326309 \times10^{-4}$
1243	$x$	$4.394631 \times10^{-2}$	$4.430123 \times10^{-2}$	$4.949229 \times10^{-2}$	$3.176514 \times10^{-2}$
	$y$	$1.305206 \times10^{-2}$	$1.314090 \times10^{-2}$	$1.490195 \times10^{-2}$
	$z$	$2.951794 \times10^{-2}$	$2.986286 \times10^{-2}$	$3.465152 \times10^{-2}$

Step-VIII (Spike Diagram): After that we implement our custom made python code to plot the geodesic spike diagram. The spike diagram is basically tracing all the points on the surface of a sphere (finite number of points) from the center of a sphere derived from the geodesic equations. In the later section we will discuss in detail about this diagram.

Step-IX (Mesh creation for the geodesic open sphere): Using python’s inbuilt codes we are able to construct the surface formed by the end points of the geodesics derived in the earlier step. We call this closed surface a geodesic sphere or geodesic open sphere.

Step-X (Analysis on the newly formed geodesic sphere): Now we do certain analysis by answering questions like: How much it is different from the Euclidean sphere? Which side of the geodesic sphere is compressed and which side is expanded compared to the Euclidean sphere? What is the relative volume / surface area difference between them?

Configurations

Device Configurations

We are using free version of Google Colab (Python 3 Google Compute Engine backend (GPU)) T4 GPU with system RAM 12.7 GB, GPU RAM 15 GB, Disk 112.6 GB. However, the code can be executed through Python 3 CPU, but it would take a little long time for the whole evaluation.

Architecture of the Physics-Informed Neural Network

Our Physics-Informed Neural Network consists of 6 layers in total in which 1 input layer (with 1 neuron $\tau$), 4 hidden layers each equipped with 64 neurons each and 1 output layer (with 3 neurons $x(\tau)$, $y(\tau)$, $z(\tau)$). Only two types of derivatives are involved $\frac{d}{d\tau}$ and $\frac{d^2}{d\tau^2}$ and they are used under automatic differentiation. Three ODEs are involved with 6 initial conditions (3 conditions for initial points and 3 conditions for initial velocities foe each of the output variables). We used Adam optimizer with learning rate $10^{-3}$ and number of iterations are 10000 for each execution. Our loss function is given by

where $\mathcal{L}_p$ denotes the loss in physics (i.e., loss in learning the governing equations) and $\mathcal{L}_b$ is the loss in learning the boundary/initial conditions. The $\lambda_p,\lambda_b$ are the weights for the physics loss and boundary loss functions. We are letting the network learn the ODEs by the physics and boundary/initial conditions only, and we do not invoke any experimental data points. We only use the SciPy data to compare the results with the PINN. We are executing the neural network for 10 random seeds starting from 1234 and up to 1243. This makes the network add different weights to the loss function and from here we find the best seed by choosing the minimum overall root mean square error (RMSE). We used Tanh as activation function and Glorot Uniform as initializer. The network architecture diagram is given here for readers convenience.

Results and Discussions

Computation time

First let us mention the computation time for the execution of the Python script on Google Colab T4 GPU using the same configuration as mentioned in § 2.2.

• Calculating symbolic Christoffel symbol, Riemann, Ricci and Scalar curvature, Einstein manifold checking and Flatness checking: 6.220 sec

• Finding geodesic equations: 0.240 sec

• Installation of DeepXDE: 10.100 sec

• PINN 10 seeds each with 10k iterations total time: 315.777 sec

• Volume comparison : 60.05 sec

Estimated Total time: 6 minutes 32.387 seconds.

Remark 1. It should be mentioned here that the parameters we have chosen for this experiment are not random. They are chosen after multiple execution of the program and several failures. Sometimes the method took long time to respond and sometimes it was faster but the errors were drastic. After a careful observation and keeping the time of execution in mind we have finalized the parameter values.

Results on the geometry of the metric

We are considering $(M,g)$ as our manifold where $M = \{(x,y,z) \in\mathbb{R}^{3}, x > 0\}$, $x,y,z$ as the coordinates in the manifold and the metric $g$ is defined by

So the line element becomes

Remark 2. *We have chosen this region so that $g$ remains positive definite, which is necessary for our experiment. Also the reason behind choosing such simple metric is that we had no idea about the outcome. To make sure we have very good control over the metric we are changing just one diagonal element and keeping the remaining the same as Euclidean metric. In future, we have plan to investigate over complex metrics, degenerate and indefinite metrics along with metrics with singularity.

The nonzero components of the Christoffel symbols for $g$ are

All of the Riemann, Ricci and Scalar curvature for the metric $g$ are zero. So the manifold is a flat manifold.

The geodesic equations for $g$ are given by

To solve the system (22) we assume that at $\tau\in[0,2]$ and at $\tau=0$ we have

Results of PINN

Now we apply both SciPy Integrate Solve IVP with DOP853 method (which is 8th order Runge-Kutta Numerical Solver with 5th order error estimator) and our PINN to get the solution. We have done random seed experiments starting seen number from 1234 up to 1243 and the following graphs shows the test-train loss for each experiment.

The Table 2 shows the error metrics of each experiment. As you can see the best performance on the basis of minimum overall RMSE is achieved by the first seed 1234 (See Table 2). So interested reader can use this seed to reconstruct the same result or may find another seed for better outcome. The Table 3 gives a comparison of final position at $\tau=2$ predicted by Numerical solution and PINN solution. Tables are placed at the end of the section for better readability.

Table 3. Comparison at $\tau= 2$ (geodesic end points PINN vs Numerical).

Method	$x(2)$	$y(2)$	$z(2)$
Numerical	3.60555128	2.58800260	3.00000000
PINN	3.60543630	2.58794780	3.00008150

Now let us see the solution comparison between PINN and Numerical method. See Figure 3. The PINN solution lines almost coincides with te numerical solution.

To understand how close is the PINN prediction we show the box plots of the errors generated at each seed.

Honestly, we did not expected the model to produce such a small box plot for each run. It is very good to get such small errors for each of the run except for the last one. This shows the model’s stability over multiple runs. We suspect that further tuning the model’s configuration may lead to a more stabilized prediction. For this run the overall RMSE is $1.144660 \times10^{-4}$. If one runs the code again, they may get a different RMSE and different seed.

Results on Geodesic Spike Diagram

Now we are ready to visualize the geodesic spike diagram. Let $p$ be a point fixed in $\mathbb{R}^{3}$ and let $S(p,1)$ be the unit sphere around the point $p$ i.e., $S(p,1)=\{x\in\mathbb{R}^{3}: \lVert x-p\rVert_{2}=1\}$. First we discretize the sphere into finite number of points. Then draw the geodesics (produced by $g$) starting from the center $p$ to all the other points on the surface. This will create a spike like diagram. For the Euclidean metric, the shortest paths are straight lines. If we draw the same for Euclidean metric then all we are going to see a bunch of finite sticks glues at one point. Topologically speaking it is equivalent to $[0,1]^{100}/\sim$ where $\sim$ identifies the 0 of each interval. This is what we mean by spike diagram. For arbitrary Riemannian metric we are interested is the diagram, this will give us an idea what does the unit sphere will look like under the metric $g$.

For this visualization we are creating 15 points along the $\theta$ coordinate and 30 points in the $\phi$ coordinate on the unit sphere. As a total we have 450 points in the surface to which geodesics are going to be connected. However, we will see 422 geodesics in the plot. Why? To visualize the geodesics emanating from the center (fixed point), we discretize the unit sphere of initial directions using $N_{\theta}=15$ polar subdivisions and $N_{\phi}=30$ azimuthal subdivisions. Since all azimuthal angles represent the same direction at each pole, only one direction is retained at each pole. Hence, the total number of distinct initial directions is

Each initial direction determines a geodesic initial-value problem, which is numerically solved using SciPy’s DOP853 method. Thus, at most 422 geodesic trajectories are obtained, depending on the success of numerical integration. This number depends on the chosen angular resolution and does not represent the total number of geodesics on the manifold.

Construction of Geodesic Sphere

Now we make a triangular mesh to construct the surface formed by the displaced points of the deformed sphere. We call this the geodesic sphere. Since the metric we have considered in (2) is very close to Euclidean metric, so plotting the open sphere is of no use as the human eye can not tell the difference between this image and the Euclidean sphere. To explicitly visualize the change we compare the newly generated geodesic sphere with Euclidean unit sphere $S^2$. First comparison comes in form of compression and decompression. Let us define the displacement at the $i$-th point as

where $x_i^{\mathrm{geodesic}}$ is the endpoint of the $i$-th geodesic, $x_i^{\mathrm{Euclidean}}$ is the corresponding point on the reference Euclidean sphere, $\Delta_i$ is the Euclidean coordinate distance between the two corresponding points. Plotting the heat map on the geodesic sphere gives us a nice visualization.

Spatial distance comparison

This graph (Figure 6) gives an idea about the parts of the open sphere of the Riemannian metric where it has a difference with the Euclidean unit sphere. Next we consider the signed radial deviation. The signed radial deviation measures how far a geodesic endpoint lies from the reference sphere in the radial direction. It is defined as

where $p=(2,2,2)$ is the common center of the two constructions and $R$ is the radius of the reference Euclidean sphere, here $R=1$. The sign of $\delta_i$ indicates the endpoint’s radial position:

1. $\delta_i>0$: the geodesic endpoint lies outside the reference sphere.

2. $\delta_i<0$: the geodesic endpoint lies inside the reference sphere.

3. $\delta_i=0$: the geodesic endpoint lies on the reference sphere.

Before going into the signed distance deviation let us see the effect of $\Delta_i$ in the Euclidean sphere.

This image is on the Euclidean sphere, and this represents which part has a change when the metric is invoked. Now for the signed distance deviation heatmap on the Euclidean sphere we present the following diagram.

It can be seen from here that contraction is happening when the metric is involved in $\mathbb{R}^3$. This is what we are looking for, this visualization answers one of the fundamental questions which is “what does the Riemannian metric does to the unit sphere?” Answer: Figure 8 is what is does to the sphere.

The next question comes to mind is, Can we compare the change in volume and surface for this newly generated geodesic sphere with the Euclidean sphere? Yes we can and now we are going to do the same.

Since the surface of the geodesic sphere is made of triangular meshes so we can numerically compute all the necessary information.

Intrinsic Riemannian Area and Volume Comparison

Till now we have $(M,g)$ as a three-dimensional Riemannian manifold with metric tensor $g=(g_{ij})$. For a geodesic sphere of radius $R$ centered at $p \in M$, we compare its intrinsic area and volume with those of a Euclidean sphere and sphere of the same radius. Unlike Euclidean mesh calculations, which use the vertex coordinates and triangular faces, intrinsic Riemannian quantities are determined by the metric tensor $g$.

The intrinsic area of the geodesic sphere is calculated using the induced metric:

where

and $X(\theta,\phi)$ parameterizes the geodesic sphere.

The intrinsic volume of the geodesic sphere is calculated by integrating the Riemannian volume density:

where $F(s,\theta,\phi)$ is the geodesic polar map and

is its coordinate Jacobian matrix. These formulas assume that the geodesic polar map is regular and injective over the sphere.

For comparison, the Euclidean reference values are ($R = 1$),

The percentage differences are

Thus, Euclidean mesh area and volume describe geometry in ambient Euclidean coordinates, whereas intrinsic Riemannian area and volume incorporate the metric tensor $g$.

Parameter	Value
Center	(2, 2, 2)
Geodesic radius	1.0
Base grid ($N_\theta$, $N_\phi$, $N_r$)	(30, 60, 40)
Surface Area
Euclidean reference area	12.5663706144
Riemannian sphere area	12.5208785515
Difference	−0.0454920628
Percentage difference	−0.362014%
Volume
Euclidean reference volume	4.1887902048
Riemannian sphere volume	4.1735465537
Difference	−0.0152436511
Percentage difference	−0.363915%
Grid-Refinement Diagnostics
Refined grid ($N_\theta$, $N_\phi$, $N_r$)	(45, 90, 60)
Refined sphere area	12.5461407142
Refined sphere volume	4.1820115974
Area change relative to base	+0.201354%
Volume change relative to base	+0.202416%

This shows that at the base resolution, the computed Riemannian sphere area and sphere volume are approximately 0.362014% and 0.363915% smaller than their respective Euclidean reference values. Grid refinement changes the area and volume by approximately 0.201354% and 0.202416%, respectively. These changes indicate numerical sensitivity. Interested readers may investigate on further refinement to establish the convergence. This numerical computations justifies why there are so many blue regions in the Figure 8.

Mathematical Assumptions we have taken.

• The metric must remain smooth and positive-definite throughout the sampled region.

• The exponential map must remain injective over the geodesic sphere.

• The coordinate Jacobian is approximated using finite differences, so the results should be checked at successively finer resolutions.

Numerical limitations. The computed area and volume depend on the chosen grid resolution and the numerical approximation of the coordinate Jacobian. Although grid refinement produces relatively small changes in the computed quantities, these results alone do not establish convergence. Further refinement and error analysis are required to assess the numerical accuracy of the intrinsic area and volume estimates.

Future Works and Conclusion

Future Works

Several extensions of this work are possible. The PINN framework may be investigated for more general Riemannian metrics and higher-dimensional manifolds. Further improvements in neural network architecture, training strategies, and numerical integration may also help improve the accuracy and stability of the computed geodesics. Future work will additionally focus on systematic grid-refinement studies to establish the convergence of the computed intrinsic area and volume as well as on degenerate metrics. For our experiment we solely focused on positive definite metric, thus experimenting with indefinite metrics, or metrics with singularities will lead us to another dimension of this study.

Conclusion

In this work, we presented a computational framework for visualizing geodesics and geodesic spheres associated with an abstract Riemannian metric on $\mathbb{R}^{3}$. The geodesic equations were solved using both SciPy and a Physics-Informed Neural Network (PINN), and their solutions were compared to assess the accuracy of the neural network approximation. The resulting geodesic sphere was compared with the Euclidean reference sphere using displacement visualization and intrinsic area and volume calculations. It was found that, at the base grid resolution, the computed intrinsic area and volume were approximately 0.362% and 0.364% smaller than their Euclidean reference values, respectively. Grid refinement produced small but nonzero changes in these estimates, indicating the need for further numerical convergence studies.

Overall, this work illustrates how numerical differential equation solvers and PINNs can be combined to study and visualize geodesic structures associated with Riemannian metrics. The approach provides a computational framework for exploring geometric properties that may be difficult to visualize directly.

Data availability: The author(s) declare that all the data generated during the study are in the manuscript and no additional data was generated. The python script used for this study is available in the GitHub repository [16]

Funding: The author(s) declare that no funding was received for this work.

Declaration of Generative AI and AI-assisted technologies: The author(s) declares that no generative AI or AI assisted technology is used for the completion of the manuscript.

Conflict of interest: The author(s) declares no conflict of interest.

Acknowledgement: The author(s) acknowledges all the reviewers for giving their valuable suggestions towards the improvement of the manuscript. Sujit Bhattacharyya

Department of Mathematics, Sister Nivedita University, DG 1/2 New Town, Action Area 1, Kolkata - 700156, India

Email: [email protected]