Parallel Transport Along a Curve on a Surface

We parallel transport along a nice, yet non-trivial, curve.

Here's we implement this fundamental operation in differential geometry on a surface. First we gather all imports & plot a surface.

# import necessary libraries
import plotly.graph_objects as go  # main interactive plotting library
import numpy as np # numerical operations
import os # operating system interactions
import plotly.io as pio # plotly input/output (saving, rendering)
import plotly.offline as pyo # plotly offline mode (exporting to html)
import sympy as sym # symbolic mathematics for derivatives
# import math # mathematical functions
# to reference variables within markdown
from myst_nb import glue
import matplotlib.pyplot as plt # some simple plots
# function to give components of curve
# use a scipy ODE solver
from scipy.integrate import solve_ivp

# Initialize Plotly offline mode, so we can display plots of the notebook in the .html file
pio.renderers.default = 'notebook'
pyo.init_notebook_mode()

# surface parameters ---
R = 1.0  # Radius of the surface/torus
# express the angles in terms of t
t_sym = sym.symbols('t')
phi_t_expr = 1-t_sym/np.pi
phi_t = sym.lambdify(t_sym, phi_t_expr, 'numpy')
theta_t_expr = t_sym  
theta_t = sym.lambdify(t_sym, theta_t_expr, 'numpy')

# --- 1.5 Curve Parameter ---
theta_0 = np.pi / 4 # Constant polar angle for the curve

# --- 2. Generate surface Surface Data 
phi = np.linspace(0, 2 * np.pi, 100)
theta = np.linspace(0, 2*np.pi, 100)
phi_sym, theta_sym = sym.symbols('phi, theta')
phi_grid, theta_grid = np.meshgrid(phi, theta)

# numeric parametrizations for surfaces - - - - - - - - - - - - -
# -- sphere --
x_surface_num = R * np.sin(phi_grid) * np.cos(theta_grid)
y_surface_num = R * np.sin(phi_grid) * np.sin(theta_grid)
z_surface_num = R * np.cos(phi_grid)

# -- torus --
x_surface_num = (R + 0.3 * np.cos(phi_grid)) * np.cos(theta_grid)
y_surface_num = (R + 0.3 * np.cos(phi_grid)) * np.sin(theta_grid)
z_surface_num = 0.3 * np.sin(phi_grid)

# symbolic parametrizations for surfaces - - - - - - - - - - - - -
# -- sphere --
x_surface = R * sym.sin(phi_sym) * sym.cos(theta_sym)
y_surface = R * sym.sin(phi_sym) * sym.sin(theta_sym)
z_surface = R * sym.cos(phi_sym)

# -- torus --
x_surface = (R + 0.3 * sym.cos(phi_sym)) * sym.cos(theta_sym)
y_surface = (R + 0.3 * sym.cos(phi_sym)) * sym.sin(theta_sym)
z_surface = 0.3 * sym.sin(phi_sym)

x_surface_func = sym.lambdify((phi_sym, theta_sym), x_surface, 'numpy')
y_surface_func = sym.lambdify((phi_sym, theta_sym), y_surface, 'numpy')
z_surface_func = sym.lambdify((phi_sym, theta_sym), z_surface, 'numpy')


# Create the surface Surface Trace
surface_trace = go.Surface(
    x=x_surface_num, y=y_surface_num, z=z_surface_num,
    colorscale='Earth',  # Uniform blue color
    opacity=0.6,
    showscale=False,  # Hide the color bar
    name='surface'
)

fig0 = go.Figure(data=[surface_trace])
fig0.update_layout(
    title='Parallel Transport of a Tangent Vector along a Curve in a surface',
    scene=dict(
        xaxis_title='X',
        yaxis_title='Y',
        zaxis_title='Z',
        # Enforce cube-like aspect ratio (important for spherical geometry)
        aspectmode='data',
        # Set the axis ranges
        xaxis=dict(range=[scale*xmin-buffer, scale*xmax+buffer]),
        yaxis=dict(range=[scale*ymin-buffer, scale*ymax+buffer]),
        zaxis=dict(range=[scale*zmin-buffer, scale*zmax+buffer]),
    ),
    margin=dict(l=0, r=0, b=0, t=30)
)
fig0.show()

Tangent & Normal Vectors

# symbolic partial derivatives
dx_dphi_surface   = sym.diff(x_surface, phi_sym)
dy_dphi_surface   = sym.diff(y_surface, phi_sym)
dz_dphi_surface   = sym.diff(z_surface, phi_sym)
dx_dtheta_surface = sym.diff(x_surface, theta_sym)
dy_dtheta_surface = sym.diff(y_surface, theta_sym)
dz_dtheta_surface = sym.diff(z_surface, theta_sym)

# unit tangent symbolic expressions (normalize using symbolic sqrt)
v_phi_surface   = 1/sym.sqrt(  dx_dphi_surface**2 + dy_dphi_surface**2   + dz_dphi_surface**2)
v_theta_surface = 1/sym.sqrt(dx_dtheta_surface**2 + dy_dtheta_surface**2 + dz_dtheta_surface**2)

# symbolic unit tangent components dx/dphi, dy/dphi, dz/dphi
T_x_phi_surface = v_phi_surface * dx_dphi_surface
T_y_phi_surface = v_phi_surface * dy_dphi_surface
T_z_phi_surface = v_phi_surface * dz_dphi_surface

# symbolic unit tangent components dx/dtheta, dy/dtheta, dz/dtheta
T_x_theta_surface = v_theta_surface * dx_dtheta_surface
T_y_theta_surface = v_theta_surface * dy_dtheta_surface
T_z_theta_surface = v_theta_surface * dz_dtheta_surface

# numeric functions for unit tangent components dx/dphi, dy/dphi, dz/dphi
T_x_phi_func = sym.lambdify((phi_sym, theta_sym), T_x_phi_surface, 'numpy')
T_y_phi_func = sym.lambdify((phi_sym, theta_sym), T_y_phi_surface, 'numpy')
T_z_phi_func = sym.lambdify((phi_sym, theta_sym), T_z_phi_surface, 'numpy')

# numeric functions for unit tangent components dx/dtheta, dy/dtheta, dz/dtheta
T_x_theta_func = sym.lambdify((phi_sym, theta_sym), T_x_theta_surface, 'numpy')
T_y_theta_func = sym.lambdify((phi_sym, theta_sym), T_y_theta_surface, 'numpy')
T_z_theta_func = sym.lambdify((phi_sym, theta_sym), T_z_theta_surface, 'numpy')

# tangent basis vectors
T_theta = sym.Matrix([T_x_theta_surface, T_y_theta_surface, T_z_theta_surface])
T_phi = sym.Matrix([T_x_phi_surface, T_y_phi_surface, T_z_phi_surface])

# numeric functions for tangent basis vectors
T_theta_func = sym.lambdify((phi_sym, theta_sym), T_theta, 'numpy')
T_phi_func = sym.lambdify((phi_sym, theta_sym), T_phi, 'numpy')

# Evaluate tangent vectors on grid
Tx_theta= T_x_theta_func(phi_grid, theta_grid)
Ty_theta= T_y_theta_func(phi_grid, theta_grid)
Tz_theta= T_z_theta_func(phi_grid, theta_grid)*np.ones(phi_grid.shape)

# Evaluate tangent vectors on grid
Tx_phi= T_x_phi_func(phi_grid, theta_grid)
Ty_phi= T_y_phi_func(phi_grid, theta_grid)
Tz_phi= T_z_phi_func(phi_grid, theta_grid)*np.ones(phi_grid.shape)

# Symbolic normal vector components
Nx_sym = T_y_phi_surface * T_z_theta_surface - T_z_phi_surface * T_y_theta_surface
Ny_sym = T_z_phi_surface * T_x_theta_surface - T_x_phi_surface * T_z_theta_surface
Nz_sym = T_x_phi_surface * T_y_theta_surface - T_y_phi_surface * T_x_theta_surface

# Normal vector is cross product of tangents --  --  --  --  --  --  --  --  --  --
Nx_num = Ty_phi * Tz_theta - Tz_phi * Ty_theta
Ny_num = Tz_phi * Tx_theta - Tx_phi * Tz_theta
Nz_num = Tx_phi * Ty_theta - Ty_phi * Tx_theta

# normalize
norm_N = np.sqrt(Nx_num**2 + Ny_num**2 + Nz_num**2)
Nx_num = Nx_num / norm_N
Ny_num = Ny_num / norm_N
Nz_num = Nz_num / norm_N
N = sym.Matrix([Nx_sym, Ny_sym, Nz_sym])
# make a numeric function for N
N_func = sym.lambdify((phi_sym, theta_sym), N, 'numpy')

# verify orthogonality
dot_product_1 = Tx_phi * Nx_num + Ty_phi * Ny_num + Tz_phi * Nz_num
dot_product_2 = Tx_theta * Nx_num + Ty_theta * Ny_num + Tz_theta * Nz_num
print("Max dot product (T_phi . N): ", np.max(np.abs(dot_product_1)))
print("Max dot product (T_theta . N): ", np.max(np.abs(dot_product_2)))

# verify orthogonality at a point
print('should be zero = ',np.dot(N_func(1,1).flatten(),T_theta_func(1,1).flatten()))
print('should be zero = ',np.dot(N_func(3,1).flatten(),T_phi_func(3,1).flatten()))

# time parameter
t_sym = sym.symbols('t')

# express n(t) in terms of t
n_t_expr = sym.Matrix([Nx_sym.subs({phi_sym:phi_t_expr, theta_sym:theta_t_expr}),
                       Ny_sym.subs({phi_sym:phi_t_expr, theta_sym:theta_t_expr}),
                       Nz_sym.subs({phi_sym:phi_t_expr, theta_sym:theta_t_expr})])
n_t_fn = sym.lambdify(t_sym, n_t_expr, 'numpy')

# verify orthogonality at a random t
rng = np.random.default_rng()
random_t = rng.uniform(low=0, high=2*np.pi)
print('should be zero = ',np.dot(n_t_fn(random_t).flatten(),
                                T_theta_func(phi_t(random_t),theta_t(random_t)).flatten()))

# create fast numeric function for n(t)
def n_fast(t):
    return np.asarray(n_t_fn(t)).astype(float).flatten()

# matrix of n'(t): differentiating n_t_expr wrt t_sym
n_t_prime_expr = sym.Matrix([sym.diff(n_t_expr[0], t_sym),
                             sym.diff(n_t_expr[1], t_sym),
                             sym.diff(n_t_expr[2], t_sym)])

# create fast numeric function for n'(t)
n_t_prime_fn = sym.lambdify(t_sym, n_t_prime_expr, 'numpy')
def n_prime_fast(t):
    return np.asarray(n_t_prime_fn(t)).astype(float).flatten()

# extract initial point coordinates

x0 = x_surface_func(phi_t(0), theta_t(0))
y0 = y_surface_func(phi_t(0), theta_t(0))
z0 = z_surface_func(phi_t(0), theta_t(0))

# extract initial tangent vectors
tx0 = T_phi_func(phi_t(0), theta_t(0))[0][0]
ty0 = T_phi_func(phi_t(0), theta_t(0))[1][0]  
tz0 = T_phi_func(phi_t(0), theta_t(0))[2][0]

tangents_trace_1 = go.Cone(
    x=[x0], y=[y0], z=[z0], # Start point
    u=[tx0],v=[ty0], w=[tz0], # vector components
    sizemode='absolute',
    sizeref=.1, # Adjust arrow size
    colorscale=[[0, 'blue'], [1, 'yellow']],
    name='Initial Vector',
    showscale=False
)

# extract initial tangent vectors
tx1 = T_theta_func(phi_t(0), theta_t(0))[0][0]
ty1 = T_theta_func(phi_t(0), theta_t(0))[1][0]  
tz1 = T_theta_func(phi_t(0), theta_t(0))[2][0]

tangents_trace_2 = go.Cone(
    x=[x0], y=[y0], z=[z0], # Start point
    u=[tx1],v=[ty1], w=[tz1], # vector components
    sizemode='absolute',
    sizeref=.1, # Adjust arrow size
    colorscale=[[0, 'yellow'], [1, 'yellow']],
    name='Initial Vector',
    showscale=False
)

# extract initial normal vectors
nx0 = n_fast(0)[0]
ny0 = n_fast(0)[1]  
nz0 = n_fast(0)[2]

# normal vector
normal_trace = go.Cone(
    x=[x0], y=[y0], z=[z0], # Start point
    u=[-nx0],v=[-ny0], w=[-nz0], # vector components
    sizemode='absolute',
    sizeref=.05, # Adjust arrow size
    colorscale=[[0, 'red'], [1, 'red']],
    name='Normal Vector',
    showscale=False
)

# prepare to plot
xmin, xmax = min(x_surface_num.flatten()), max(x_surface_num.flatten())
ymin, ymax = min(y_surface_num.flatten()), max(y_surface_num.flatten())
zmin, zmax = min(z_surface_num.flatten()), max(z_surface_num.flatten())
buffer = 0.1
scale = 1.5

fig0 = go.Figure(data=[surface_trace, tangents_trace_1,tangents_trace_2,normal_trace])

# Set up the layout
fig0.update_layout(
    title='Parallel Transport of a Tangent Vector along a Curve in a surface',
    scene=dict(
        xaxis_title='X',
        yaxis_title='Y',
        zaxis_title='Z',
        # Enforce cube-like aspect ratio (important for spherical geometry)
        aspectmode='data',
        # Set the axis ranges
        xaxis=dict(range=[scale*xmin-buffer, scale*xmax+buffer]),
        yaxis=dict(range=[scale*ymin-buffer, scale*ymax+buffer]),
        zaxis=dict(range=[scale*zmin-buffer, scale*zmax+buffer]),
    ),
    margin=dict(l=0, r=0, b=0, t=30)
)
fig0.show()
pio.write_html(fig0, 'torus_tangents.html', auto_open=True)

Max dot product (T_phi . N):  1.6653345369377348e-16
Max dot product (T_theta . N):  1.1102230246251565e-16
should be zero =  2.7755575615628914e-17
should be zero =  0.0
should be zero =  5.551115123125783e-17

The derivatives of the surface tangent and normal vectors.

In case we wish to plot a vector field later, these will be useful.

# derivative of the unit tangent vector along phi direction
dTx_dphi_surface = sym.diff(T_x_phi_surface, phi_sym)
dTy_dphi_surface = sym.diff(T_y_phi_surface, phi_sym)
dTz_dphi_surface = sym.diff(T_z_phi_surface, phi_sym)
dTx_dtheta_surface = sym.diff(T_x_theta_surface, theta_sym)
dTy_dtheta_surface = sym.diff(T_y_theta_surface, theta_sym)
dTz_dtheta_surface = sym.diff(T_z_theta_surface, theta_sym)

# numeric functions for derivative of unit tangent vector along phi direction
dTx_dphi_func = sym.lambdify((phi_sym, theta_sym), dTx_dphi_surface, 'numpy')
dTy_dphi_func = sym.lambdify((phi_sym, theta_sym), dTy_dphi_surface, 'numpy')
dTz_dphi_func = sym.lambdify((phi_sym, theta_sym), dTz_dphi_surface, 'numpy')

# numeric functions for derivative of unit tangent vector along theta direction
dTx_dtheta_func = sym.lambdify((phi_sym, theta_sym), dTx_dtheta_surface, 'numpy')
dTy_dtheta_func = sym.lambdify((phi_sym, theta_sym), dTy_dtheta_surface, 'numpy')
dTz_dtheta_func = sym.lambdify((phi_sym, theta_sym), dTz_dtheta_surface, 'numpy')

# symbolic derivative of normal vector components
dNx_dphi_surface = sym.diff(Nx_sym, phi_sym)
dNy_dphi_surface = sym.diff(Ny_sym, phi_sym)
dNz_dphi_surface = sym.diff(Nz_sym, phi_sym)
dNx_dtheta_surface = sym.diff(Nx_sym, theta_sym)
dNy_dtheta_surface = sym.diff(Ny_sym, theta_sym)
dNz_dtheta_surface = sym.diff(Nz_sym, theta_sym)

# numeric functions for derivative of normal vector components
dNx_dphi_func = sym.lambdify((phi_sym, theta_sym), dNx_dphi_surface, 'numpy')
dNy_dphi_func = sym.lambdify((phi_sym, theta_sym), dNy_dphi_surface, 'numpy')
dNz_dphi_func = sym.lambdify((phi_sym, theta_sym), dNz_dphi_surface, 'numpy')
dNx_dtheta_func = sym.lambdify((phi_sym, theta_sym), dNx_dtheta_surface, 'numpy')
dNy_dtheta_func = sym.lambdify((phi_sym, theta_sym), dNy_dtheta_surface, 'numpy')
dNz_dtheta_func = sym.lambdify((phi_sym, theta_sym), dNz_dtheta_surface, 'numpy')

Plot Curve & Initial Vector

Next we plot the curve we wish to parallel transport over and the initial vector we wish to parallel transport. We display the initial tangent vector as a cone.

# components of the curve
def alpha_x(t):
    # -- sphere curves --
    # return R * sym.sin(theta_0) * sym.cos(t) # og
    # loxodrome return R * sym.cos(1-t/np.pi) * sym.cos(t)
    # nonsimple 5-petal return R * sym.sin(beta) * sym.cos(t)
    # return R * sym.cos(1-t/np.pi) * sym.cos(5*t) # nice loxodrome
    # -- torus curves --
    phi = 1-t/np.pi
    theta = t
    dphi_dt = -1/np.pi
    dtheta_dt = 1
    return (R + 0.3 * sym.cos(phi)) * sym.cos(theta), dphi_dt, dtheta_dt
def alpha_y(t):
    # -- sphere curves --
    # return R * sym.sin(theta_0) * sym.sin(t) # og
    # loxodrome return R * sym.cos(1-t/np.pi) * sym.sin(t)
    # non-simple 5-petal return R * sym.sin(beta) * sym.sin(t)
    # return R * sym.cos(1-t/np.pi) * sym.sin(5*t) # nice loxodrome
    # -- torus curves --
    phi = 1 - t / np.pi
    theta = t
    dphi_dt = -1 / np.pi
    dtheta_dt = 1
    return (R + 0.3 * sym.cos(phi)) * sym.sin(theta), dphi_dt, dtheta_dt

def alpha_z(t):
    # -- sphere curves --
    # return R * sym.cos(theta_0) # og
    # loxodrom return R * sym.sin(1-t/np.pi)
    # non-simple 5-petal return R * sym.cos(beta) # nice loxodrome
    # return R * sym.sin(1-t/np.pi)
    # -- torus curves --
    phi = 1 - t / np.pi
    dphi_dt = -1 / np.pi
    dtheta_dt = 0
    return (0.3 * sym.sin(phi)), dphi_dt, dtheta_dt

# symbolic variable
t = sym.symbols('t')

axt, dphi_dt, dtheta_dt = alpha_x(t)
ayt, _, _ = alpha_y(t)
azt, _, _ = alpha_z(t)

# the curve in numpy
# the numberic (numpy) components of the curve
cx = sym.lambdify(t,axt,'numpy')
cy = sym.lambdify(t,ayt,'numpy')
cz = sym.lambdify(t,azt,'numpy')

# the symbolic derivatives of the components
v_x = sym.diff(axt,t)
v_y = sym.diff(ayt,t)
v_z = sym.diff(azt,t)
v = 1/sym.sqrt(v_x**2 + v_y**2 + v_z**2)

# symbolic unit tangent vector
T_x = v*v_x
T_y = v*v_y
T_z = v*v_z

# numeric unit tangent vector components
Tx = sym.lambdify(t,T_x,'numpy')
Ty = sym.lambdify(t,T_y,'numpy')
Tz = sym.lambdify(t,T_z,'numpy')

# numeric velocity components
vx = sym.lambdify(t,v_x,'numpy')
vy = sym.lambdify(t,v_y,'numpy')
vz = sym.lambdify(t,v_z,'numpy')

# need to compute n(t) = d/dt of v(t)
# should be surface normal, not accelleration vector
acc_x = sym.diff(v_x,t)
acc_y = sym.diff(v_y,t)
acc_z = sym.diff(v_z,t)

## the numeric components of accelleration
nx = sym.lambdify(t,acc_x,'numpy')
ny = sym.lambdify(t,acc_y,'numpy')
nz = sym.lambdify(t,acc_z,'numpy')

# derivative of accelleration
np_x = sym.diff(acc_x,t)
np_y = sym.diff(acc_y,t)
np_z = sym.diff(acc_z,t)

# the numberic (numpy) derivatives of the components of the full
# #accelleration vector
npx = sym.lambdify(t,np_x,'numpy')
npy = sym.lambdify(t,np_y,'numpy')
npz = sym.lambdify(t,np_z,'numpy')

# parameterize curve
s = np.linspace(0, 2 * np.pi, 100)

# for example, the following contains all x-coords of the velocity vector
# print(v_x(s))

# convert the array s to a list before feeding the lambdified function
x_curve = cx(s) 
y_curve = cy(s)
z_curve = cz(s) * np.ones(s.shape[0]) # sometimes cz returns a scalar

# Create the Curve Trace
curve_trace = go.Scatter3d(
    x=x_curve, y=y_curve, z=z_curve,
    mode='lines',
    line=dict(color='red', width=5),
    name=r"Curve"
)

# --- 4. Generate Initial Tangent Vector Data (Yellow) ---
t_point = 0  # Point at phi = 45 degrees
# choose the starting point as a base point
xp = x_curve[0]
yp = y_curve[0]
zp = z_curve[0]


# Scale the vector for to get the right length visualization
scale_factor = 1
norm_Vx = Tx(0) * scale_factor
norm_Vy = Ty(0) * scale_factor
norm_Vz = Tz(0) * scale_factor

# Create the Vector Trace
vector_trace = go.Cone(
    x=[xp], y=[yp], z=[zp], # Start point
    u=[norm_Vx], v=[norm_Vy], w=[norm_Vz], # Vector components
    sizemode='absolute',
    sizeref=0.2, # Adjust arrow size
    colorscale=[[0, 'green'], [1, 'green']],
    name='Initial Vector',
    showscale=False
)



fig0 = go.Figure(data=[surface_trace, tangents_trace_1,tangents_trace_2,normal_trace,
                      curve_trace, vector_trace])
# Set up the layout
fig0.update_layout(
    title='Parallel Transport of a Tangent Vector along a Curve in a surface',
    scene=dict(
        xaxis_title='X',
        yaxis_title='Y',
        zaxis_title='Z',
        # Enforce cube-like aspect ratio (important for spherical geometry)
        aspectmode='data',
        # Set the axis ranges
        xaxis=dict(range=[scale*xmin-buffer, scale*xmax+buffer]),
        yaxis=dict(range=[scale*ymin-buffer, scale*ymax+buffer]),
        zaxis=dict(range=[scale*zmin-buffer, scale*zmax+buffer]),
    ),
    margin=dict(l=0, r=0, b=0, t=30)
)
fig0.show()

Compute parallel transported vector

The condition for a vector to be parallel transported is $D_u w = 0$. This in turn gives

$\nabla_v \bar{w} = (\nabla_v \bar{w} \cdot \bar{n}) \bar{n}$ where $v$ is the unit tangenct vector.

Using the shape operator this is equivalent to:

$\nabla_v \bar{w} = (-\bar{w} \cdot \nabla_v \bar{n} ) \bar{n}$

The derivatives are just regular time derivatives and this gives on ODE.

$\bar{w}'(t) = (-\bar{w}(t)) \cdot \bar{n}'(t)) \bar{n}$

To set-up an ODE that matches our given constraints we need to compute $\bar{n}'(t)$ and $\bar{n}(t)$ along the curve.

# plot w0 and n_fast
# Create the Vector Trace
curve_normal_trace = go.Cone(
    x=[xp], y=[yp], z=[zp], # Start point
    u=[n_fast(0)[0]], v=[n_fast(0)[1]], w=[n_fast(0)[2]], # Vector components
    sizemode='absolute',
    sizeref=0.2, # Adjust arrow size
    colorscale=[[0, 'purple'], [1, 'purple']],
    name='Normal Vector',
    showscale=False
)
fig0 = go.Figure(data=[surface_trace, tangents_trace_1,tangents_trace_2,normal_trace,
                      curve_trace, vector_trace,curve_normal_trace])

# Set up the layout
fig0.update_layout(
    title='Parallel Transport of a Tangent Vector along a Curve in a surface',
    scene=dict(
        xaxis_title='X',
        yaxis_title='Y',
        zaxis_title='Z',
        # Enforce cube-like aspect ratio (important for spherical geometry)
        aspectmode='data',
        # Set the axis ranges
        xaxis=dict(range=[scale*xmin-buffer, scale*xmax+buffer]),
        yaxis=dict(range=[scale*ymin-buffer, scale*ymax+buffer]),
        zaxis=dict(range=[scale*zmin-buffer, scale*zmax+buffer]),
    ),
    margin=dict(l=0, r=0, b=0, t=30)
)

fig0.show()

Define the ODE

# Initial Conditions
#w0 = np.array([Tx(0),Ty(0),Tz(0)])                # Initial state vector = unit tangent vector but could be anything
w0 = np.array([T_x_phi_func(phi_t(0),theta_t(0)),
               T_y_phi_func(phi_t(0),theta_t(0)),
               T_z_phi_func(phi_t(0),theta_t(0))])              

# visualize w0
x0 = float(x_surface.subs({phi_sym:phi_t(0), theta_sym:theta_t(0)}))
y0 = float(y_surface.subs({phi_sym:phi_t(0), theta_sym:theta_t(0)}))
z0 = float(z_surface.subs({phi_sym:phi_t(0), theta_sym:theta_t(0)}))


w0_trace = go.Cone(
    x=[x0], y=[y0], z=[z0], # Start point
    u=[w0[0]], v=[w0[1]], w=[w0[2]], # Vector components
    sizemode='absolute',
    sizeref=0.2, # Adjust arrow size
    colorscale=[[0, 'orange'], [1, 'orange']],
    name='Initial w Vector',
    showscale=False
)

fig0 = go.Figure(data=[surface_trace, tangents_trace_1,tangents_trace_2,normal_trace,
                      curve_trace, vector_trace,curve_normal_trace, w0_trace])
# Set up the layout
fig0.update_layout(
    title='Parallel Transport of a Tangent Vector along a Curve in a surface',
    scene=dict(
        xaxis_title='X',
        yaxis_title='Y',
        zaxis_title='Z',
        # Enforce cube-like aspect ratio (important for spherical geometry)
        aspectmode='data',
        # Set the axis ranges
        xaxis=dict(range=[scale*xmin-buffer, scale*xmax+buffer]),
        yaxis=dict(range=[scale*ymin-buffer, scale*ymax+buffer]),
        zaxis=dict(range=[scale*zmin-buffer, scale*zmax+buffer]),
    ),
    margin=dict(l=0, r=0, b=0, t=30)
)
fig0.show()

# Define the ODE function
def vector_ode_derivative(t, w, n, n_prime):
    """
    Defines the ODE system w'(t) = -dot(w(t), n'(t)) * n.

    Args:
        t (float): Current time.
        w (np.ndarray): Current state vector w(t).
        n_prime_func (callable): Function returning n'(t).
        n_func (callable): Function returning n(t).

    Returns:
        np.ndarray: The derivative w'(t).
    """
    n_prime_t = n_prime(t)
    n_t = n(t)
    dot_product = -np.dot(w, n_prime_t)
    dw_dt = dot_product * n_t
    return dw_dt

Solve the ODE

# constants
k = 2*np.pi                    # Upper bound for time
L = 100                        # number of points to derive a solution
t_span = (0, k)                # Time span (start, end)
t_eval = np.linspace(t_span[0], t_span[1], L) # Points at which to store the solution

# solve the ODE
solution = solve_ivp(vector_ode_derivative, t_span, w0, 
    args=(n_fast,n_prime_fast), t_eval=t_eval, atol=1e-6, rtol=1e-6)

# Time points
time_points = solution.t

# Solution vectors at each time point
w_solution = solution.y.T # Transpose to get w(t) as rows

Incorporate Solution

Next, we use the solution in our plots. We use a shaft and cone to display the vector. The shaft part is created using Scatter3D in line mode. Because the way Scatter3D seems to work, we need to great a different trace for each shaft we draw. Because there are many vectors in the vector field solution we have many segments (portions of vectors) to plot. We do that in a loop.

wx = w_solution[:,0]    # the solution x,y,z to the ODE sytem
wy = w_solution[:,1]    
wz = w_solution[:,2]

xpn = cx(time_points) # the new basepoint along curve
ypn = cy(time_points)
zpn = cz(time_points) * np.ones(np.size(xpn))

norm_w = np.sqrt(wx**2 + wy**2 + wz**2) # uniformly scale the solution vex
norm_w= np.where(norm_w == 0, 1.0, norm_w) # avoid division by zero
scale_factor = 1
wx = wx / norm_w * scale_factor
wy = wy / norm_w * scale_factor
wz = wz / norm_w * scale_factor

e = 0.3                            # arrow shaft scale factor
dx = e*wx; dy = e*wy; dz = e*wz;   # arrow shaft vector components

# Build single Scatter3d shaft trace using NaN separators (fast)
n_pts = xpn.size
xs = np.empty(n_pts * 3)
ys = np.empty(n_pts * 3)
zs = np.empty(n_pts * 3)

xs[0::3] = xpn
xs[1::3] = xpn + e * wx
xs[2::3] = np.nan

ys[0::3] = ypn
ys[1::3] = ypn + e * wy
ys[2::3] = np.nan

zs[0::3] = zpn
zs[1::3] = zpn + e * wz
zs[2::3] = np.nan

shaft_trace = go.Scatter3d(
    x=xs, y=ys, z=zs,
    mode='lines',
    line=dict(color='blue', width=2),
    showlegend=False
)

# vectorized dot products (no loop)
axp_all = vx(time_points)
ayp_all = vy(time_points)
azp_all = vz(time_points)
dots = wx[0] * axp_all + wy[0] * ayp_all + wz[0] * azp_all



# list to hold all traces of segments
segment_traces = []   
dots = np.ones(xpn.shape[0])
for i in range(xpn.shape[0]):
    x_start = xpn[i] # the x-coordinates of the i-th basepoints 
    y_start = ypn[i]
    z_start = zpn[i]

    x_end = x_start + e*wx[i] # the x-coords of the i-th endpoints
    y_end = y_start + e*wy[i]
    z_end = z_start + e*wz[i]

    segment_traces.append(go.Scatter3d(
        
        x=[x_start, x_end],
        y=[y_start, y_end],
        z=[z_start, z_end],
        mode='lines', # to plot lines not points
        line=dict(color='blue', width=2),
        showlegend=False
    ))
    # dots[i] = np.dot([wx[i],wy[i],wz[i]],[v_x(i),v_y(i),v_z(i)])
    dots[i] = np.dot([wx[i],wy[i],wz[i]],[w0[0],w0[1],w0[2]])

Angle between transport and velocity

The angles between the transported vector and the tangent vector does not remain constant when the path is a loxodrome, but is oscillatory.

plt.title("Angle between PT-vector & Unit Tangent over Time")
plt.plot(time_points, 180/np.pi*np.arccos(dots), color='blue')
# --- 7. Compute the Holonomy of the loop
h = 180/np.pi * 2*np.pi*(1-np.sin(np.pi/2 - theta_0))
mydot = np.dot([wx[-1],wy[-1],wz[-1]],[w0[0],w0[1],w0[2]])
print('dot = ', mydot)
H = np.arccos(mydot)
#plt.text(2,20,
#    r'Holonomy = , %.2f degrees' % (180/np.pi*H), 
#    ha = 'center', va = 'bottom', fontsize=10)
#plt.text(3,30,
#    r'Holonomy Est: R(L) = 2pi(1 - sin(pi/2 - phi))) = , %.2f degrees' % h, 
#    ha = 'center', va = 'bottom', fontsize=10)

plt.show()

dot =  -0.4161459680827364

The Parallel Transported Vector Field Tip

This is done using a Cone from ?graphical-objects. The syntax is self-explanatory. The (1-e) scalar is chosen so that the total length of the cone+shaft equals the scaled length of w.

# vector field
new_new_vector_trace = go.Cone(
    x=xpn + e*wx, y=ypn+e*wy, z=zpn+e*wz, # Start point
    u=(1-e)*wx, v=(1-e)*wy, w=(1-e)*wz, # Vector components
    sizemode='absolute',
    sizeref=.4, # Adjust arrow size
    colorscale='Viridis',
    name='Parallel Transported Vector',
    showscale=False,
    showlegend=True,
    anchor = 'tail',
    opacity = 0.6
)

The accelleration vector

In case we might want to see it is here. Just include the trace in the final fig() at the end.

# --- 4.5 Add accelleration vector
time_point = 0  # Point at phi = 45 degrees

# Vector components V (Vx, Vy, Vz)
Vx = -R * np.sin(theta_0) * np.sin(time_point) 
Vy = R * np.sin(theta_0) * np.cos(time_point)
Vz = 0

# Unit Tangent
vector_magnitude = np.sqrt(Vx**2 + Vy**2 + Vz**2)
norm_Vx = Vx / vector_magnitude
norm_Vy = Vy / vector_magnitude
norm_Vz = Vz / vector_magnitude

# Time-derivative of the vector w (in this case = V)
# Full accelleration vector
Vxg = -R * np.sin(theta_0) * np.cos(time_point) 
Vyg = -R * np.sin(theta_0) * np.sin(time_point)
Vzg = 0

# Create the Vector Trace using Cone
new_vector_trace = go.Cone(
    x=[xpn], y=[ypn], z=[zpn], # Start point
    u=[Vxg], v=[Vyg], w=[Vzg], # Vector components
    sizemode='absolute',
    sizeref=0.2, # Adjust arrow size
    colorscale=[[0, 'steelblue'], [1, 'steelblue']],
    name='Accelleration Vector',
    showlegend=True
)

Final Plot Commands

The segment_traces is a list of traces and not a trace itself, so needs to be combined as shown (concatenated?) instead of separated by commas.

xmin, xmax = min(x_surface_num.flatten()), max(x_surface_num.flatten())
ymin, ymax = min(y_surface_num.flatten()), max(y_surface_num.flatten())
zmin, zmax = min(z_surface_num.flatten()), max(z_surface_num.flatten())
buffer = 0.1
scale = 1.5

# --- 5. Create Figure and Layout ---
fig = go.Figure(data=[surface_trace, curve_trace, vector_trace, new_new_vector_trace]+segment_traces)
#fig.update_zaxes(range = [zmin, zmax])

# efficient fig = go.Figure(data=[surface_trace, curve_trace, vector_trace, new_new_vector_trace,shaft_trace])
# for trouleshooting
#fig = go.Figure(data=[surface_trace] + segment_traces)

# Set up the layout
fig.update_layout(
    title='Parallel Transport of a Tangent Vector along a Curve in a surface',
    scene=dict(
        xaxis_title='X',
        yaxis_title='Y',
        zaxis_title='Z',
        # Enforce cube-like aspect ratio (important for spherical geometry)
        aspectmode='data',
        # Set the axis ranges
        xaxis=dict(range=[scale*xmin-buffer, scale*xmax+buffer]),
        yaxis=dict(range=[scale*ymin-buffer, scale*ymax+buffer]),
        zaxis=dict(range=[scale*zmin-buffer, scale*zmax+buffer]),
    ),
    margin=dict(l=0, r=0, b=0, t=30)
)

fig.show()

# --- 7. Compute the Holonomy of the loop in case of a sphere
#h = 180/np.pi * 2*np.pi*(1-np.sin(np.pi/2 - theta_0))
#print('holonomy = R(L) = 2pi(1 - sin(pi/2 - phi))) = ', h)
#glue("holonomy", h)

Compute the Holonomy of the Loop

The holonomy is {glue:}holonomy

# Define a path in a known writable location (e.g., user's home directory)
output_directory = os.path.expanduser("~/g/teaching/ma541/python/") # Or any other writable path
output_file_path = os.path.join(output_directory, "surface_pt.html")

try:
    #pio.renderers.default = "jupyterlab" # or "notebook", "iframe", "vscode"
    pio.renderers.default = "notebook" # "iframe", "vscode"
    pio.write_html(fig, output_file_path, auto_open=True)
    print(f"Plotly figure saved successfully to: {output_file_path}")
except OSError as e:
    print(f"Error saving Plotly figure: {e}")
    print("Please ensure you have write permissions to the specified directory.")
# Save the figure as an HTML file

# - - - - run in terminal to convert notebook to html with dark theme
# jupyter nbconvert --execute --to html /path/to/examplej.ipynb --HTMLExporter.theme=dark

Plotly figure saved successfully to: /Users/josh/g/teaching/ma541/python/surface_pt.html