Differential geometry

Directional derivative

use math.la: vector

function diffvh(h)
   return |v,f,x| (f(x+h*v)-f(x-h*v))/(2*h)
end

diffv = diffvh(0.001)

# The standard basis
e0 = vector(1,0)
e1 = vector(0,1)


# A scalar field
f = |[x,y]| x*y

# Value and partial derivatives at (x=2,y=4)
p = vector(2,4)
print(f(p))
print(diffv(e0,f,p))
print(diffv(e1,f,p))


# A vector field
F = |[x,y]| vector(x*y,x+y)

# Value and partial derivatives at (x=2,y=4)
p = vector(2,4)
print(F(p))
print(diffv(e0,F,p))
print(diffv(e1,F,p))


# Tangential plane of f at p
T = |f,p| |x| f(p)+diffv(x-p,f,p)

p = vector(2,4)
g = T(f,p)
G = T(F,p)


# Note that diffv may be used as Gâteaux differential.
# One just needs to define the arithmetic operations
# for functions:
Function.add = |f;g| |x| f(x)+g(x)
Function.sub = |f;g| |x| f(x)-g(x)
Function.rmul= |r;f| |x| r*f(x)
Function.div = |f;r| |x| f(x)/r

Gradient

use math.la: unit, vector

# Nabla operator
function nabla(f,p)
   n = p.shape[0]
   E = list(0..n-1).map(|k| unit(n,k))
   return vector(*E.map(|e| diffv(e,f,p)))
end

# Tangential plane of f at p
function T(f,p)
   a = nabla(f,p)
   y = f(p)
   return |x| y+a*(x-p)
end


f = |[x,y]| x*y
g = T(f,vector(2,4))


# Taking the gradient of a vector valued function
# is possible, because we use the polymorphic array
# data type and thus can have vectors of vectors
# (i.e. jagged arrays).
F = |[x,y]| vector(x*y,x+y)
G = T(F,vector(2,4))

Jacobian matrix

use math.la: unit, vector, matrix


# The Jacobian matrix is the transposed gradient
function Jacobian(F,p)
   return matrix(*nabla(F,p).list.map(|v| v.list)).T
end


# After simplification we obtain
function Jacobian(F,p)
   n = p.shape[0]
   E = list(0..n-1).map(|k| unit(n,k))
   return matrix(*E.map(|e| diffv(e,F,p).list)).T
end


# Tangential plane of F at p
function T(F,p)
   J = Jacobian(F,p)
   y = F(p)
   return |x| y+J*(x-p)
end

Volume

use math.na: integral
use math.la.inversion: det
use math: pi, sin, cos

function volume(F,[x0,x1],[y0,y1],[z0,z1])
   return integral(x0,x1,|x| integral(y0,y1,|y| integral(z0,z1,|z|
      abs(det(Jacobian(F,vector(x,y,z))))
   )))
end

torus = |{R}| |[r,p,t]| vector(
   (R+r*cos(p))*cos(t),
   (R+r*cos(p))*sin(t),
   r*sin(p)
)

F = torus(R=2)
V = volume(F,[0,1],[0,2*pi],[0,2*pi])
print(V)

# R=2; r=1
# V = 2*pi^2*r^2*R

Area

use math.na: integral
use math.la.inversion: det
use math: sqrt, pi, sin, cos

function metric_tensor(F)
   return fn|p|
      J = Jacobian(F,p)
      return J.T*J
   end
end

function area(F,[x0,x1],[y0,y1])
   g = metric_tensor(F)
   return integral(x0,x1,|x| integral(y0,y1,|y|
      sqrt(det(g(vector(x,y))))
   ))
end

torus = |{R,r}| |[p,t]| vector(
   (R+r*cos(p))*cos(t),
   (R+r*cos(p))*sin(t),
   r*sin(p)
)

F = torus(R=2,r=1)
A = area(F,[0,2*pi],[0,2*pi])
print(A)

# R=2; r=1
# A = 4*pi^2*r*R