Hallo there,
I also saw the equation in Griffiths. After some searching I found a satisfactory explanation.
Here is a short derivation.
There is a more elaborate explanation in the attachment!
This equation is the
Minimal surface equation derived by
Lagrange in 1762 using the
Euler-Lagrange equation.
One can also derive the equation directly by using
calculus of variations to minimize surface area. The simplest derivation with this method I found
here. (I have also done this in the attachment.)
Physical interpretation
Griffiths talks about an equation describing a stretched rubber sheet. Though, traditionally this equation is used to describe soap films. Namely, surface tension ##\sigma \doteq FL^-1 \doteq JL^-2## is in units of energy per area. As energy is minimized, the surface area of soap films is minimized.
Rubber approximates the equation if it is not to thick and shearing stresses are neglected. Namely, thick rubber will have additional resistance to bending and surface tension is only a "pulling" and not a "turning" force. If we assume the rubber is flat enough bending forces and twisting forces will be negligible.
Surface area equation
Firstly, we want to find an equation for the total surface area given some height function / surface function ##h(x,y)## on a domain ##x,y\in\Omega\subset\mathbb{R}## with a boundary ##\partial \Omega##. (You can think of the domain as a shadow directly below the surface and the boundary as the edge of that shadow.) A nice visualization and derivation you can find
here.
Essentially, you chop up the the surface into area elements ##\Delta A## with a base of sides ##\Delta x## by ##\Delta y##. Then, the area is spanned by two vectors ##\Delta \textbf{u} = (\Delta x, 0, \frac{\partial f}{\partial x})^\top## and##\Delta \textbf{v} = (0, \Delta y, \frac{\partial f}{\partial x})^\top##. Thus,
$$\Delta A = \Delta |u\times v|= \sqrt{1+(\frac{\partial f}{\partial x})^2+(\frac{\partial f}{\partial y})^2} $$ and therefore the area of a function ##f(x,y)## over a domain ##x,y\in \Omega \subset \mathbb{R}## is given by $$A_\Omega (f) = \int_\Omega \sqrt{1+(\frac{\partial f}{\partial x})^2+(\frac{\partial f}{\partial y})^2} dxdy$$.
Minimizing functional by solving the Euler-Lagrange equation
The
Euler-Lagrange equation for a function ##f(x,y)## with two variables that minimizes a
functional of the form $$I[f] = \int_\Omega \mathcal{L}(x,y,f,\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y})dxdy$$ is $$\frac{\partial \mathcal{L}}{\partial f} - \frac{d}{dx}(\frac{\partial \mathcal{L}}{\partial f_x})-\frac{d}{dy}(\frac{\partial \mathcal{L}}{\partial f_y})=0 \; \; | \; \; f_x = \frac{\partial f}{\partial x}, \; f_y = \frac{\partial f}{\partial y}$$
Identifying that the solution to the area functional ##A_\Omega## should be the optimal surface \ height function ##h(x,y)##.
Deriving minimum surface equation
Suppose we have some function ##h: \Omega \rightarrow \mathbb{R}^2## defined on a domain ##\Omega \subset \mathbb{R^2}## with boundary ##\partial \Omega## and this function minimizes the area functional ##A_\Omega (h)##. Then, for the set of all possible functions on this domain ##f:\Omega \rightarrow \mathbb{R}^2## with the same boundary conditions ##h(x_b,y_b)=f(x_b,y_b) \; \forall x_b,y_b\in \partial \Omega## we have $$\textrm{min}\; A_\Omega(f) = A_\Omega (h)$$.
Now, we identify ##\mathcal{L} = \sqrt{1+(\frac{\partial h}{\partial x})^2+(\frac{\partial h}{\partial y})^2}## must satisfy the Euler-Lagrange equation. As $$\frac{\partial \mathcal{L}}{\partial h} - \frac{d}{dx}(\frac{\partial \mathcal{L}}{\partial h_x})-\frac{d}{dy}(\frac{\partial \mathcal{L}}{\partial h_y})=0 \; \; | \; \; h_x = \frac{\partial h}{\partial x}, \; h_y = \frac{\partial h}{\partial y}$$ $$\frac{\partial \mathcal{L}}{\partial h}=0, \; \frac{\partial \mathcal{L}}{\partial h_x}=\frac{h_x}{\mathcal{L}}, \; \frac{\partial \mathcal{L}}{\partial h_y}=\frac{h_y}{\mathcal{L}}$$ we conclude $$ \frac{\partial}{\partial x} \left( \frac{\frac{\partial h}{\partial x}}{\sqrt{1+(\frac{\partial h}{\partial x})^2+ (\frac{\partial h}{\partial y})^2}} \right) + \frac{\partial}{\partial y} \left( \frac{\frac{\partial h}{\partial y}}{\sqrt{1+(\frac{\partial h}{\partial x})^2+ (\frac{\partial h}{\partial y})^2}} \right) = 0$$.
or more compactly $$\nabla \cdot \left( \frac{\nabla h}{\sqrt{1+|\nabla h|^2}} \right) = 0$$ which is even generalizable for a gradient of any order, as explained in the attachment.
Two dimensional equivalent
We can similarly solve for the two dimensional case. Assume we have a function ##f(x,y)## in domain ##x,y\in \Omega \subset \mathbb(R)## with boundary ##(x_0,x_1) = \partial \Omega##. Now we will minimize the line functional $$L(f)=\int_{x_0}^{x_1} \sqrt{1+(\frac{\partial f}{\partial x})^2} dx$$.
Defining ##\textrm{min}\; L_\Omega(f) = L_\Omega (y)## we find $$\frac{\partial}{\partial x} \left( \frac{y_x}{\sqrt{1+(y_x)^2}} \right)=0$$ Now, we can solve this for the given boundary conditions $$\frac{y_x}{\sqrt{1+(y_x)^2})} = c$$ $$y_x=c\sqrt{1+(y_x)^2} \rightarrow (y_x)^2 =\frac{c^2}{1-c^2}$$ concluding $$y(x)=c_0 x+c_1 \rightarrow y(x) = (y(x_1)-y(x_0))x+x_0$$ The line functional / arc length is optimized for a straight line between the two boundary points.
Similarity to Laplace's equation
We can obtain Laplace's equation from the minimum surface equation by assuming the surface is almost planar. Namely, then ##(\frac{\partial h}{\partial x})^2 +(\frac{\partial h}{\partial y})^2\approx 0## and the equation will reduce to ##\nabla^2 h=0##
