Nom!
As part of a bigger project, I've was writing some code to calculate the scalar curvature of surfaces of the form $z = f(x,y)$. This uses a general calculation of the scalar curvature to produce a function ScalarCurvature that calculates curvature in for surfaces with this specific form. However, I'm finding that some (but not all) of my test cases fail...
I suspect the issue comes from how I built a function from the output of the calculations, but the more standard way was failing to evaluate derivatives completely. I'm pretty sure the formula I obtain from the calculations is correct, it is just a matter of how it is put into a function.
Behaviour
What I'm finding is that for a sphere
ScalarCurvature[Sqrt[1 - (x^2 + y^2)], x, y] // FullSimplify
gives
2
as it should. But, for a Pringle
Plot3D[x^2 - y^2, {x, -1, 1}, {y, -1, 1},
RegionFunction -> Function[{x, y, z}, x^2 + y^2 < 1]]
ScalarCurvature[x^2 - y^2, x, y]
yields
0
when it should, by my calculation, be
- 8 / (1 + 4x^2 + 4y^2)^2
Code
Here's the code, I'm pretty sure it is something to do with the bottom few lines
(* Formula for intrinsic curvature of a function *)
(* STEP 1 - metric tensor *)
(* Find g such that: g dY dY = δ dX dX *)
dX = {dx, dy, D[f[x, y], x] dx + D[f[x, y], y] dy};
dY = {dx, dy};
long = dX.dX // FullSimplify;
a = Coefficient[long, dx^2];
b = Coefficient[long, dx dy]/2;
c = Coefficient[long, dy^2];
g = {{a, b}, {b, c}}
gInv = Inverse[g]
short = dY.g.dY;
(* Check *)
long == short // FullSimplify
(* STEP 2 - Christoffel symbols *)
cristoffelLower = Module[{v = {x, y}},
Table[
D[g[[c]][[a]], v[[b]]] +
D[g[[c]][[b]], v[[a]]] -
D[g[[a]][[b]], v[[c]]],
{c, 1, 2}, {a, 1, 2}, {b, 1, 2}]]/2;
(* Seconds kind - these appear to be correct*)
cristoffelHigher =
Table[
Plus @@ Table[
cristoffelLower[[d]][[a]][[b]] gInv[[c]][[d]], {d, 1, 2}],
{c, 1, 2}, {a, 1, 2}, {b, 1, 2}] // FullSimplify
(* STEP 3 - Curvature Tensors *)
reimann = Module[{v = {x, y}},
Table[
D[cristoffelHigher[[ρ]][[ν]][[σ]], v[[μ]]] -
D[cristoffelHigher[[ρ]][[μ]][[σ]], v[[ν]]] +
Plus @@ Table[
cristoffelHigher[[ρ]][[μ]][[λ]]
cristoffelHigher[[λ]][[ν]][[σ]], {λ, 1, 2}]
- Plus @@ Table[
cristoffelHigher[[ρ]][[ν]][[λ]]
cristoffelHigher[[λ]][[μ]][[σ]], {λ, 1, 2}],
{ρ, 1, 2}, {σ, 1, 2}, {μ, 1, 2}, {ν, 1,
2}]] // FullSimplify;
(* Ricci curvature, 1st position is contravariant, 3rd is covariant *)
ricci = Table[
Plus @@ Table[
reimann[[k]][[i]][[k]][[j]], {k, 1, 2}],
{i, 1, 2}, {j, 1, 2}];
(* Scalar curvature, both positions are covariant *)
scalar =
Plus @@ Table[
Plus @@ Table[
ricci[[i]][[j]] gInv[[j]][[i]],
{j, 1, 2}], {i, 1, 2}] // FullSimplify
(* Define a function to do the calculation *)
ScalarCurvature[fun_, xx_, yy_] :=
scalar /. Derivative[i_, j_][f][x, y] -> D[fun, {xx, i}, {yy, j}]
GaussianCurvature[f_, x_, y_] := ScalarCurvature[f, x, y]
