Calculating curvature of a contour

Question

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$.

Answer 1

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

Answer 2

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<n$ 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 ν.

Answer 3

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.