Laplace equation in polar coordinates, using matrices I saw that the method shown below could be used to derive the Laplace equation for polar coordinates using less calculations.
\begin{aligned}
&\nabla^{2} u=\frac{\partial^{2} u}{\partial r^{2}}+\frac{1}{r} \frac{\partial u}{\partial r}+\frac{1}{r^{2}} \frac{\partial^{2} u}{\partial \theta^{2}}\\
&\left(\begin{array}{c}
\frac{\partial}{\partial r} \\
\frac{\partial}{\partial \theta}
\end{array}\right)=\left(\begin{array}{ll}
\frac{\partial x}{\partial r} & \frac{\partial y}{\partial r} \\
\frac{\partial x}{\partial \theta} & \frac{\partial x}{\partial \theta}
\end{array}\right)\left(\begin{array}{c}
\frac{\partial}{\partial x} \\
\frac{\partial}{\partial y}
\end{array}\right) \quad \longrightarrow \quad\left(\begin{array}{c}
\frac{\partial}{\partial x} \\
\frac{\partial}{\partial y}
\end{array}\right)=\left(\begin{array}{ll}
\frac{\partial x}{\partial r} & \frac{\partial y}{\partial r} \\
\frac{\partial x}{\partial \theta} & \frac{\partial x}{\partial \theta}
\end{array}\right)^{-1}\left(\begin{array}{c}
\frac{\partial}{\partial r} \\
\frac{\partial}{\partial \theta}
\end{array}\right)\\
&\nabla^{2}=\left(\begin{array}{c}
\frac{\partial}{\partial x} \\
\frac{\partial}{\partial y}
\end{array}\right) \cdot\left(\begin{array}{c}
\frac{\partial}{\partial x} \\
\frac{\partial}{\partial y}
\end{array}\right)
\end{aligned}
Even after computing the following matrix:
\begin{pmatrix}\frac{\partial x}{\partial r} & \frac{\partial y}{\partial r} \\
\frac{\partial x}{\partial \theta} & \frac{\partial x}{\partial \theta}
\end{pmatrix}
(using $x=r\cos\theta, y=r\sin\theta$), I still do not know how the Laplacian $\nabla^{2} u=\frac{\partial^{2} u}{\partial r^{2}}+\frac{1}{r} \frac{\partial u}{\partial r}+\frac{1}{r^{2}} \frac{\partial^{2} u}{\partial \theta^{2}}$ is gotten. What especially confuses me is the dot product and how the $2\times1$ operator gets applied to a $2 \times2$ matrix.
 A: First of all, the idea behind the operator being multiplied on the right by a matrix is that it will still output an operator. For example, operating on a function $f$ we would find
$$\begin{aligned}
\left(\begin{array}{c}
\frac{\partial}{\partial r} \\
\frac{\partial}{\partial \theta}
\end{array}\right)[f]=\left(\begin{array}{ll}
\frac{\partial x}{\partial r} & \frac{\partial y}{\partial r} \\
\frac{\partial x}{\partial \theta} & \frac{\partial x}{\partial \theta}
\end{array}\right)\left(\begin{array}{c}
\frac{\partial}{\partial x} \\
\frac{\partial}{\partial y}
\end{array}\right)[f] = \left(\begin{array}{ll}
\frac{\partial x}{\partial r} & \frac{\partial y}{\partial r} \\
\frac{\partial x}{\partial \theta} & \frac{\partial x}{\partial \theta}
\end{array}\right)\left(\begin{array}{c}
\frac{\partial f}{\partial x} \\
\frac{\partial f}{\partial y}
\end{array}\right) 
\end{aligned} \; .$$
Just note that this vector is really a functional that takes an input.
As for your main question, it looks like you're starting to derive it using the metric tensor $g_{ij}.$ It can be thought of as a matrix. You have already written down the Jacobian matrix (or it's transpose depending on who you ask). We will say
$$J^T=\begin{bmatrix}
\frac{\partial x}{\partial r} & \frac{\partial y}{\partial r} \\
\frac{\partial x}{\partial \theta} & \frac{\partial x}{\partial \theta}
\end{bmatrix}
$$
and define $$g_{ij} = J^TJ$$ to be a symmetric, indexable object that gives the components of the resulting matrix from the multiplication $J^TJ\;.$ Now let $g^{ij}$ denote the matrix inverse of $g_{ij}$ and let $g$ denote the matrix determinant of $g_{ij}.$
From this framework, we can use this "metric" $g_{ij}$ to define $$\Delta f = \frac{1}{\sqrt{g}} \sum_i \sum_j \frac{\partial }{\partial x_i}\bigg[  \sqrt{g} \; g^{ij} \frac{\partial f}{\partial x_j}  \bigg] \; $$ as the Laplacian of $f.$ Do note that you should then take $(x_1,x_2) = (r,\theta) .$ Also, $g_{ij}$ works out quite nicely for polar coordinates and when all is said and done is
$$g_{ij} = \begin{bmatrix}  1 & 0\\ 0 & r^2  \end{bmatrix} \; .$$
Let me know if this needs any further clarification.
$\textbf{EDIT:}\;$ clarification to the OP
The Laplacian is defined as I have written it above for a general coordinate system on any pseudo-Riemannian manifold. This can be worked out for general coordinates as defined above or for a specific coordinate system. As an example (which I will not fully work out) we can use polar coordinates as is relevant to your question.
Consider a scalar function of polar coordinates $f:(r,\theta)\rightarrow\mathbb{R} \;.$ We know that in Cartesian coordinates that the Laplacian is defined as
$$\Delta f = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} \;.$$
We can use the chain and product rules to expand this as
$$\Delta f = \frac{\partial^2 f}{\partial r^2}\bigg[ \frac{\partial r}{\partial x} \bigg]^2 + \frac{\partial f}{\partial r}\frac{\partial^2 r}{\partial x^2} + \frac{\partial^2 f}{\partial \theta^2}\bigg[ \frac{\partial \theta}{\partial x} \bigg]^2 + \frac{\partial f}{\partial \theta}\frac{\partial^2 \theta}{\partial x^2} \\+ 
\frac{\partial^2 f}{\partial r^2}\bigg[ \frac{\partial r}{\partial y} \bigg]^2 + \frac{\partial f}{\partial r}\frac{\partial^2 r}{\partial y^2} + \frac{\partial^2 f}{\partial \theta^2}\bigg[ \frac{\partial \theta}{\partial y} \bigg]^2 + \frac{\partial f}{\partial \theta}\frac{\partial^2 \theta}{\partial y^2} \; .
$$
Using the standard polar coordinates $$\begin{align}
x&=r\cos(\theta)\\
y&=r\sin(\theta)
\end{align}$$
we can evaluate the expression above and it had better come out to be what we expect (and it does). We can do something similar if you prefer the Laplacian definition $$\Delta f = \nabla \cdot \nabla f = \text{div}(\nabla f)$$
but then we would have to work out the polar gradient and polar divergence. If you feel unconvinced, I recommend that you work one of these out and see that it yields the same results as the summation definition I gave above.
