What is the best way to make notation in multi-variable calculus rigorous? (Warning: long post worrying about formalism, may seem totally pointless to many.)
Related to "Can the idea of a 'function of a variable' be made rigorous?".
There is a disconnect in how mathematics is usually formalized and how mathematicians tend to write about derivatives. For example, it's not uncommon to see notation like:
\begin{alignat}{1}    v &= a^2b + 3ab
\\[5pt] \frac{\partial v}{\partial a} &= 2ab + 3b
\\[5pt] \frac{\partial v}{\partial b} &= a^2 + 3a
\\[5pt] a &= x^2 + y^2
\\[5pt] b &= xy
\\[5pt] \frac{\partial v}{\partial x} &= \frac{\partial v}{\partial a}\frac{\partial a}{\partial x} + \frac{\partial v}{\partial b}\frac{\partial b}{\partial x}
  \end{alignat}
This notation assumes that

*

*A function 'knows' its variables, so that $\frac{\partial v}{\partial a}$ and $\frac{\partial v}{\partial b}$ are both well-defined and different

*Different Variables refer to different objects, so that $v(a,b)$ and $v(x,y)$ are both well-defined and different

*The same object can be both a function and a parameter, i.e., $a$ appears in the formulas $v=a^2b + 3ab$ and $a=x^2+y^2$
However, in the classical formalism, a function $f : \mathbb{R} \rightarrow \mathbb{R}$, say the one that doubles all inputs, is just a set of tuples $\{(1,2), (2,4), (3,6), ...\}$. (Alternatively, it may be a triple (domain, codomain, set-of-tuples), but this difference doesn't matter for this post.) This formalism doesn't support any of the above things:

*

*If $v = \{((1,1),4), ...\}$, then $a$ and $b$ aren't objects of any kind. The function $v$ may be defineed using notation where 'a' an 'b' appear as bound variables, such as $v = \{((a,b),a^2b + 3ab) \;|\; (a,b) \in \mathbb R^2\}$, but they're only bound within the expression that defines the set $v$. This makes it "impossible" for $a = x^2 + y^2$ to refer to the $a$ that was used in $v = a^2b + 3ab$.

*For the same reason, $v(a,b)$ and $v(x,y)$ cannot refer to different objects.

*If a function's domain consists of real numbers, it doesn't consist of functions, so if $a$ and $b$ are functions, one can never evaluate $v(a,b)$.

I see three ways to deal with this disconnect

*

*Ignore it and continue to do calculations that lead to correct results. That is, until one gets confused about things like 'wait $\frac{\partial v}{\partial a}$ is a function from what to what?'. In such a case, always stare at the expression long enough to figure it out. Proceed to pretend that there is no problem. Alternatively, never get confused, somehow (do people really do this?)

*Try to interpret all notation as shorthands that talk about sets

*Think of the notation as being grounded in a different formalism

I've been doing #1.1 ever since learning about set theory, but I'm no longer satisfied with that.
The problem with #2 is that it doesn't seem to correspond to how people think about what they're doing. Consider the fact that serious mathematicians write things like $v = v(a,b)$ or "$f$ is a function of $x$, and $g$ is a function of $y$". These statements have no meaning in terms of the underlying formalism (except to state that the domain of $f$ is one-dimensional): if one writes $f(x) = x^2$, then $f$ is a well-defined set, but $x$ is nothing whatsover (except a bound variable in a formula that defines the set $f$). This would make the statement "$f$ is a function of $x$" merely a hint for how to interpret future notation, similar to a statement like "we write $p(x)$ for $p(\{x\})$". Again, this does not seem to correspond to how people are thinking about the material.
#3 is something I'm very interested in but don't know much about.
My question: is there a satisfying solution for #2 or #3? Could $\lambda$-calculus be a candidate?
 A: I find the best way is to do what differential geometers do. Here’s a short sketch in two dimensions:
I’ll write $\mathbb{A}$ (for ‘affine’) for the plane $\mathbb{R}^2$. We introduce functions $x, y : \mathbb{A} \to \mathbb{R}$ that assign to each point $P \in \mathbb A$ its coordinates $x(P)$ and $y(P)$. Because each point $P$ is determined by these coordinates, each other function $f : \mathbb A \to \mathbb R$ can be written in terms of the values $x(P)$ and $y(P)$. For example, we might write
$$
    f(P) = \sin(x(P)) + (x(P))^2 + 2 x(P) y(P) + (y(P))^2.
$$
However, we won’t; instead, we write it as
$$
    f = \sin x + x^2 + 2xy + y^2
$$
with the implicit understanding that $\sin x$ really means the composition ${\sin} \circ x$ and so on (but because $x$ is now a function, it really can’t be anything else). (Sometimes people will write
$$
    f(x, y) = \sin x + x^2 + 2xy + y^2,
$$
but I dislike that. The $f(x, y)$ serves as a reminder that $f$ is expressed in terms of $x$ and $y$ and mostly tells you that $f(a, b)$ (for actual numbers $a, b$) is the value of $f$ at the point $P$ with coordinates $x(P) = a$, $y(P) = b$.)
Now I can change the coordinate functions. For example, if I set $u = x$ and $v = x + y$, each point is also determined by the values $u(P)$ and $v(P)$. I can rewrite the function $f$ above as
$$
    f = \sin u + v^2.
$$
Note that even though $u$ is the same function as $x$, some occurrences of $x$ “vanished”. What the representation of a function $f$ looks like really depends on the pair of coordinate functions used. This is the reason why I renamed $x$: it is better to think of each coordinate function as belonging to one specific pair of coordinate functions. Here $x$ is always paired with $y$ and $u$ is always paired with $v$.
The notation $\frac{\partial f}{\partial x}$ now means: Express $f$ using the coordinate functions $x$ and $y$ (because $x$ is always paired with $y$) and then compute the derivate as usual. Technically, you compose $f$ with the inverse of the function $(x, y) : \mathbb{A} \to \mathbb{R}^2$ to get $f \circ (x, y)^{-1} : \mathbb{R}^2 \to \mathbb{R}$, compute the derivative with respect to the first coordinate and compose with $(x, y)$ to get a function on $\mathbb A$ again. Here, it is really important to know which other coordinate functions belong to $x$ because the inverse $(x, y)^{-1}$ depends on it. For example,
$$
  \frac{\partial f}{\partial x} = \cos x + 2x + 2y
$$
but
$$
  \frac{\partial f}{\partial u} = \cos u
$$
and these functions are different even though $x = u$.
In this context, the chain really holds in the form
$$
  \frac{\partial f}{\partial x} = \frac{\partial f}{\partial u} \frac{\partial u}{\partial x} + \frac{\partial f}{\partial v} \frac{\partial v}{\partial x}
$$
(if everything is differentiable). In the example,
$$
  \frac{\partial u}{\partial x} = 1,\qquad \frac{\partial f}{\partial v} = 2v,\qquad \frac{\partial v}{\partial x} = 1
$$
and combining these gives
$$
  \frac{\partial f}{\partial x} = \cos u + 2v = \cos x + 2x + 2y,
$$
as expected.
To summarize, if we work in this setting, the notation $\frac{\partial f}{\partial x}$ only makes sense if $x$ is part of a complete set of coordinates on the domain of $f$ (even though the other coordinates are not explicit in the notation). To calculate it, express $f$ in terms of these coordinates and compute the derivative as usual.
A: I'm going to take a shot at #2 myself, i.e., trying to make sense of everything in terms of set theory (where functions are triples).

*

*Function Definition: the equation $v = a^2b + 3ab$ is a shorthand for $v := [(a,b) \mapsto a^2b + 3ab]$, which is itself a shorthand for $v := \big(\mathbb R^2, \mathbb R, \{((a,b),a^2b + 3ab) \;|\; (a,b) \in \mathbb R^2\}\big)$, and analogously for any other function definition. If the equation is instead written as $v(a,b) = a^2b + 3ab$, it is still a shorthand for the same thing.


*Derivatives: We can define derivatives in terms of the $k$-th argument, i.e., we can define operators $D_k$ such that, if $f : \mathbb R^n \rightarrow \mathbb R$ (this is a well-defined restriction on $f$ since we defined it as a triple) with $k \le n$, then $D_k(f)$ is the function obtained by taking th derivative with respect to the $k$-th argument. In particular, if $f : \mathbb R \rightarrow \mathbb R$, then
$$ D_1(f) := [x \mapsto \lim_{t \rightarrow 0} \frac{f(x+t) - f(x)}{t}]$$
I believe that the $D_j$ are class functions, which would make them rigorously defined objects, but not sets.


*One-dimensional chain rule: if $f,g : \mathbb R \rightarrow \mathbb R$, then $D_1(f \circ g) = [x \mapsto D_1(f)(g(x)) \cdot D_1(g)(x)]$. Actually, in the one-dimensional case, the notation with $f'$ and $g'$ is pretty good and you wouldn't really need anything else. The problem is that it doesn't scale to several dimensions.


*Multi-dimensional chain rule: (this needs to cover the notation $v = a^2b + 3ab$ and $a = x^2 + y^2$ and $b = xy$). If $v : \mathbb R^2 \rightarrow \mathbb R$ and $a,b : \mathbb R^2 \rightarrow \mathbb R$, we define
$$f \circ (a,b) := [(x,y) \mapsto f(a(x,y), b(x,y))]$$
and have the rules
$$D_1(f \circ (a,b)) = [(x,y) \mapsto \big(D_1(f) \circ (a,b)\big)(x,y) \cdot D_1(a)(x,y) + \big(D_2(f) \circ (a,b)\big)(x,y) \cdot D_1(b)(x,y)] $$
and
$$D_2(f \circ (a,b)) = [(x,y) \mapsto \big(D_1(f) \circ (a,b)\big)(x,y) \cdot D_2(a)(x,y) + \big(D_2(f) \circ (a,b)\big)(x,y) \cdot D_2(b)(x,y)] $$
(This reminds me that I've seen similar notation in the book 'Computability and Logic' at the part where recrusive functions and substitutions have been defined, but I've never realized that it's just the formal version of what people do in regular multi-variable calculus.) The same could be extended to the case where the domain has dimension $n > 2$, using sums $\sum_{i = 1}^n$.
We can now translate the variable notation thus
\begin{alignat}{3}    v &= a^2b + 3ab & \hspace{30pt}V &:= [(a,b) \mapsto a^2b + 3ab]&
\\[5pt] \frac{\partial v}{\partial a} &= 2ab + 3b & D_1(V) &= [(a,b) \mapsto 2ab + 3b]
\\[5pt] \frac{\partial v}{\partial b} &= a^2 + 3a & D_2(V) &= [(a,b) \mapsto a^2 + 3a]
\\[5pt] a &= x^2 + y^2 & A &:= [(x,y) \mapsto x^2 + y^2]
\\[5pt] b &= xy & B &:= [(x,y) \mapsto xy]
  \end{alignat}
The important difference now is that we no longer ascribe any meaning to bound variables outside of their quantor. Any variable in the translated statements could be changed to a different previously unused symbol, and the meaning wouldn't change. For example, we could have defined $V$ as $V := [(\phi,q) \mapsto \phi^2q + 3\phi q]$ and nothing would change as $V$ is the same set either way. Conversely, writing $v = \phi^2q + 3\phi q$ in the original equations would change the interpretation. Similarly, I've defined the function-that-takes-the-role-of-$a$ with the uppercase $A$, but I could have used a different symbol as well.
Anyway, and the final equation
$$\frac{\partial v}{\partial x} = \frac{\partial v}{\partial a}\frac{\partial a}{\partial x} + \frac{\partial v}{\partial b}\frac{\partial b}{\partial x}$$
becomes
\begin{align}D_1(V \circ (A,B)) &= [(x,y) \mapsto \big(D_1(V) \circ (A,B)\big)(x,y) \cdot D_1(A)(x,y) + \big(D_2(f) \circ (A,B)\big)(x,y) \cdot D_1(B)(x,y)]\end{align}
I come out positively surprised at how close to the original notation I was able to keep this. I think this may relieve some of the existential dread I feel when I look at instances of the usual notation. But I'm still curious if there is a formalism where $v$ and $a$ and $x$ are all objects and one doesn't need to differentiate between $a$ and $A$.
A: This is more of an extended comment than an alternative answer. It all depends what you want to calculate. If you want to deal with implicit regions given by multiple equations without biased pairings/groupings of variables, then you'd need a bit more than are covered by the two answers at the time of writing.

I'll steal $a=x^{2}+y^{2}$ and $b=xy$ from the question and some notation from silver's answer. If you want to keep $a$ and $b$ paired and $x$ and $y$ paired, then Eike Schulte's answer seems like the tidiest approach, and I believe it covers the calculations that would appear in essentially all math texts.
With $x$ and $y$ paired, then there is no question that $\dfrac{\partial b}{\partial y}=x$ since $b=xy$ is suggesting a function $B:(x,y)\mapsto xy$ (or similarly $\widetilde{B}:\left(y,x\right)\mapsto xy$) and $\dfrac{\partial b}{\partial y}$ really means something like $\left(D_{2}B\right)\left(x,y\right)$ where $D_{2}$ means “partial with respect to the second argument” since $x$ is paired with $y$ (or similarly $\left(D_{1}\widetilde{B}\right)\left(y,x\right)$).
But we can do more than one calculation if we are agnostic about variable relationships. Suppose we just have $S=\left\{ \left(x,y,a,b\right):a=x^{2}+y^{2}\text{ and }b=xy\right\}$. This is a 2D surface in 4D space, and all the partial derivatives we might care about are slopes of relevant slices at various points.
One “partial derivative $\dfrac{\partial b}{\partial y}$” would be obtained as follows: Fix $x=x_{0}$ so that we have $\left\{ \left(x_{0},y,x_{0}^{2}+y^{2},x_{0}y\right):x\in\mathbb{R}\right\}$. Write the $b$ coordinate locally as a function of $y$ (i.e. $B(y)=x_{0}y$), and differentiate so that at any point $\left(x_{0},y_{0},x_{0}^{2}+y_{0}^{2},x_{0}y_{0}\right)$ we have $\dfrac{\partial b}{\partial y}=x_{0}$.
But another “partial derivative $\dfrac{\partial b}{\partial y}$” would be something like the following : Fix $a=a_{0}$ so that in a neighborhood of a point $\left(x,y,a_{0},b\right)\in S$ with $x\ne0$, we have (in generality by the implicit function theorem) $x=\pm\sqrt{a_{0}-y^{2}}$ (for an appropriate sign choice). Then we can write the $b$ coordinate locally as a function of $y$ (i.e. $B(y)=\pm y\sqrt{a_{0}-y^{2}}$) and differentiate. Thus, we would find that at $\left(x_{0},y_{0},a_{0},b_{0}\right)\in S$, $\dfrac{\partial b}{\partial y}=\pm\sqrt{a_{0}-y_{0}^{2}}\mp\dfrac{y_{0}^{2}}{\sqrt{a_{0}-y_{0}^{2}}}=\dfrac{a_{0}-2y_{0}^{2}}{\pm\sqrt{a_{0}-y_{0}^{2}}}=\dfrac{x_{0}^{2}-y_{0}^{2}}{x_{0}}=x_{0}-\dfrac{y_{0}^{2}}{x_{0}}$.
In thermodynamics (and some related fields?), it is common to distinguish these calculations by using subscripts to show the variable(s) being fixed initially. One might write $\left(\dfrac{\partial b}{\partial y}\right)_{x}=x$ and $\left(\dfrac{\partial b}{\partial y}\right)_{a}=x-\dfrac{y^{2}}{x}$.
