How does a Pringle lose its curvature?

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]
• I'm getting $-\frac{4}{\left(4 x^2+4 y^2+1\right)^2}$ if I use the Gaussian curvature routine from here. May 16 '16 at 9:26
• @J.M. Ah yes, that was a typo, I've corrected it from - to + May 16 '16 at 11:35