# Calculating curvature of a contour

I have an equation of a scalar field in the form

$$f(x, y) = x^2 + y^2 + xy + c$$

I want to find the curvature of the contour of the curve at $$f_c = f(0.5, 0.5)$$.

So I need to calculate the derivative $$\frac{\mathrm dy}{\mathrm dx}$$ and $$\frac{\mathrm d^2y}{\mathrm dx^2}$$

I can solve the equation $$f(x, y) = f_c$$ and get the derivative of $$f(x, y)$$ with respect to $$x$$.

On paper we do,

\begin{align*}\frac{d}{dx} f(x,y)&=\frac{d}{dx}(f_c)\\ 2x+2y\frac{dy}{dx}+y+x\frac{dy}{dx}&=0\\\frac{dy}{dx}&=-\frac{2x+y}{x+2y}\end{align*}

and further, do $$\frac{d}{dx}\left(\frac{dy}{dx}\right)$$ for a curvature approximation.

How can I rearrange the equation such that I can get the value of $$\frac{dy}{dx}$$ on Mathematica?

EDIT: This is what I tried so far:

Clear["Global*"]

ftest[x_, y_] = x*x + y*y + x*y + c;

dftest = (D[ftest[x, y], x]/D[ftest[x, y], y]);
d2ftest = D[dftest, x];


Solution:

-((2 x + y)/(x + 2 y)^2) + 2/(x + 2 y)


However, this will be wrong as I will be taking the ratio of partials with respect to $$x$$ and $$y$$.

• You could solve your equation for y, then differentiate. – b.gates.you.know.what Nov 2 '18 at 10:51
• Have you tried anything? Rules require you to show a minimal familiarity with Mathematica. – Kuba Nov 2 '18 at 10:51
• @b.gatessucks How do I do that.. I used D[Solve[f(x,y)=f(c)] but it returns {{0 -> InverseFunction[fs,2,2]^(1,0)[t,0]}} – Some_guy Nov 2 '18 at 11:01
• Start by solving your ftest for y. – b.gates.you.know.what Nov 2 '18 at 11:09
• So your implicit function is something like  f[x,y]:= x^2 + y^2 + x y - 3/4==0 ? – Ulrich Neumann Nov 2 '18 at 11:23

Clear[f]
f[x_, y_] := x^2 + y[x]^2 + x y[x] + c == fc

Solve[D[f[x, y], {x, 1}], y'[x]][[1, 1]] // FullSimplify
Solve[D[f[x, y], {x, 2}], y''[x]][[1, 1]] /. % // FullSimplify

y'[x] -> -(2 x + y[x])/(x + 2 y[x])
y''[x] -> -6 (x^2 + y[x] (x + y[x])/(x + 2 y[x])^3

• Thanks for this... It works great... But can you please tell me what the [[1, 1]] is for – Some_guy Nov 3 '18 at 9:38
• @Some_guy Glad it helps. You can look up Part ([[...]]). Solve here returns something in the form of {{a -> b}}, but I think we only need a -> b, so I use [[1, 1]] to pick it out. – Αλέξανδρος Ζεγγ Nov 3 '18 at 9:40
• @Some_guy One can interpret [[1, 1]] as "the first part of the first part (of something)" due to the double-fold braces (List ({...})). – Αλέξανδρος Ζεγγ Nov 3 '18 at 9:49

IIRC, the second fundamental form $$I\!I$$ of a levelset $$M = \varPhi^{-1}(\{0\})$$ of a mapping $$\varPhi \colon \mathbb{R}^n \to \mathbb{R}^m$$, $$m at point $$x \in \mathbb{R}^n$$ is given (up to sign that I use to mix up) by

$$I\!I(x) = \pm D\varPhi(x)^\dagger \, D^2\varPhi(x),$$

where $$D\varPhi(x)^\dagger$$ denotes the Moore-Penrose pseudoinverse of the Jacobian $$D\varPhi(x)$$ of $$\varPhi$$. Here, I assume that $$D\varPhi(x)$$ is surjective (but a similar formula can be derived if $$D\varPhi$$ has constant rank in a neighborhood of $$x$$).

More precisely, we have

$$I\!I(x)(u,v) = \pm D\varPhi(x)^\dagger \, D^2\varPhi(x)(u,v) \quad \text{for all u,\,v \in T_xM = \operatorname{ker}(D\varPhi(x))}$$

as $$I\!I(x)(u,v)$$ is not really meaningful for non-tangent vectors $$u, \, v \not \in T_xM$$.

This can be obtained as follows:

Φ = {x, y} \[Function] {x^2 + y^2 + x y + c};
DΦ = D[Φ[x, y], {{x, y}, 1}];
DDΦ = D[Φ[x, y], {{x, y}, 2}];
II = Transpose[LinearSolve[DΦ. Transpose [DΦ], DΦ]].DDΦ


In this 2-dimensional example, the curvature $$\kappa$$ of the level set at point $$(x,y)$$ should be (up to sign) computable as

ν = DΦ[[1]]/Sqrt[Total[DΦ^2, 2]];
τ = RotationMatrix[Pi/2].ν
κ = ν.(II.τ).τ // Simplify


(6 (x^2 + x y + y^2))/(5 x^2 + 8 x y + 5 y^2)^(3/2)

Here, ν is the normal to the levelset (point upwards with respect to the function Φ) and τ is the tangent obtained by counterclockwise 90-degree-rotation of ν.

• Thanks for the answer. Although I did get the right answer, I am lost trying to understand the math. Could you please point me to some reference to read about this. – Some_guy Nov 4 '18 at 6:57

Altho Αλέξανδρος's answer, which uses the partial derivative function D[] works, it is more appropriate for this situation to use Dt[], the function for total derivatives.

f = x^2 + y^2 + x y + c;

{yp, ypp} = {Dt[y, x], Dt[y, {x, 2}]} /.
First[Solve[Thread[{Dt[f, x], Dt[f, {x, 2}]} == 0 /. Dt[c, _] -> 0],
{Dt[y, x], Dt[y, {x, 2}]}]]
{(-2 x - y)/(x + 2 y), -((6 (x^2 + x y + y^2))/(x + 2 y)^3)}


These can now be easily plugged into the usual expression for curvature:

FullSimplify[ypp/PowerExpand[(Together[1 + yp^2])^(3/2)]]
-((6 (x^2 + x y + y^2))/(5 x^2 + 8 x y + 5 y^2)^(3/2))
`

which, up to a sign difference, matches the expression in Henrik's answer.