Starting with, $ur^2 = u_0r_0^2$, where $u$ is my dependent variable, $u_0$ is the initial/constant value of $u$ at $r_0$, I can take two derivatives with respect to the radial distance, r, and obtain the ODE,
$$\frac{\partial^2 u}{\partial r^2} + \frac{4}{r}\frac{\partial u}{\partial r}+\frac{2}{r^2}u=0\tag{1}\label{1}$$
I want to write this in terms of Cartesian coordinates to be solved using the MATLAB PDE Toolbox, and thus far I've considered three ways to do it (my questions are at the end).
OPTION 1: Hand-wavy
Given $$x=r\cos\theta,y=r\sin\theta$$ it can be shown via the chain rule that, $$\frac{\partial u}{\partial r}=\cos\theta\frac{\partial u}{\partial x}+\sin\theta\frac{\partial u}{\partial y}$$ and the 2D Laplacian in Cartesian coordinates is, $$\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = \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}$$
And from the starting equation I know $\frac{\partial^2 u}{\partial \theta^2}=0$, thus, I get, $$\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{3}{r}\left(\cos\theta\frac{\partial u}{\partial x}+sin\theta\frac{\partial u}{\partial y}\right)+\frac{2}{r^2}u=0\tag{2}\label{2}$$ $$\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{3}{x^2+y^2}\left(x\frac{\partial u}{\partial x}+y\frac{\partial u}{\partial y}\right)+\frac{2}{x^2+y^2}u=0\tag{2i}\label{2i}$$ (Note: the equations labelled with an "i" represent the forms I will code/implement in MATLAB's PDE solver)
OPTION 2: Direct-Cartesian
Say I start with $u(x^2+y^2)=C$ and simply take the Laplacian using Cartesian coordinates. I get, $$\left(x^2+y^2\right)\nabla^2u+2\nabla\left(x^2+y^2\right)\cdot\nabla u+u\nabla^2\left(x^2+y^2\right)=0$$ which, provided my algebra is correct, becomes,
$$\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{4}{x^2+y^2}\left(x\frac{\partial u}{\partial x}+y\frac{\partial u}{\partial y}\right)+\frac{4}{x^2+y^2}u=0\tag{3i}\label{3i}$$
Please note the difference in coefficients between \eqref{2i} and \eqref{3i}.
OPTION 3: Canonical?
Following another example, an honest conversion of \eqref{1} from polar to cartesian, with no short-cuts, should result in,
$$\cos^2\theta\frac{\partial^2 u}{\partial x^2}+ 2\cos\theta\sin\theta\frac{\partial^2 u}{\partial x \partial y} + \sin^2\theta\frac{\partial^2 u}{\partial y^2} + \frac{4}{r}\left(\cos\theta\frac{\partial u}{\partial x}+\sin\theta\frac{\partial u}{\partial y}\right)+\frac{2}{r^2}u=0\tag{4}\label{4}$$
$$\frac{x^2}{x^2+y^2}\frac{\partial^2 u}{\partial x^2}+ 2\frac{xy}{x^2+y^2}\frac{\partial^2 u}{\partial x \partial y} + \frac{y^2}{x^2+y^2}\frac{\partial^2 u}{\partial y^2} + \frac{4}{x^2+y^2}\left(x\frac{\partial u}{\partial x}+y\frac{\partial u}{\partial y}\right)+\frac{2}{x^2+y^2}u=0\tag{4i}\label{4i}$$
Once again, different coefficients. Here are where my questions come in...
QUESTIONS
(1) (a) Does everyone generally agree that \eqref{3i} is more rigorous/accurate than \eqref{2i}? If not, why not? (b) What are your thoughts on the correctness of \eqref{4i} vs \eqref{3i}? (c) As a side note, in terms of accuracy, MATLAB seems to work better with \eqref{2i} using linear triangular elements, but I'm not sure why. Any insight on that would also be appreciated.
(2) I would like to implement \eqref{4} but MATLAB does not implement cross derivatives directly. Instead, they must take the form,
$$\nabla\cdot\left(c\nabla u\right)$$
If $c$ is a constant, we get back the Laplacian with a constant coefficient. Although that's easy to implement, it does not strictly apply to \eqref{4}. Here's why: If $c$ is a non-constant scalar, we know,
$$\nabla\cdot\left(c\nabla u\right)=\nabla c\cdot\nabla u + c\nabla^2u$$
If $c$ is a matrix, I can do something like,
$$\nabla\cdot\left(c\nabla u\right)=\left(\nabla\cdot c\right)\cdot \nabla u + c\cdot\cdot\left(\nabla\nabla u\right)$$
plus rearranging terms to produce an equation equivalent to \eqref{4} that the solver can handle ($\nabla\nabla u$ is the Hessian). HOWEVER, \eqref{4} was obtained using the fact that $\partial()/\partial r = \cos\theta\partial()/\partial x+\sin\theta\partial()/\partial y$. That is, $\theta$ is held constant. Thus, I cannot express $c$ in terms of sines and cosines because (a) treating them as the constants that they are prohibits me from implementing them in MATLAB's format, (b) if I attempt to force MATLAB to use sines and cosines, the MATLAB solver will attempt to differentiate them anyway, which is equivalent to replacing $\cos\theta$ with $x/r$. The problem I see with using $\cos\theta=x/r$ is that I am essentially creating a circular argument, which undermines the logic behind \eqref{4}. In other words, \eqref{4i} is okay as-is but if I try to rearrange terms, causing $\nabla$ to operate on the second order term coefficients, such as $x^2/r^2$, that will introduce a fallacy into my equation and the results will be garbage. Question -> Does everyone agree that it would be wrong to implement,
$$\nabla\cdot\left(\begin{bmatrix}x^2/r^2 & xy/r^2\\xy/r^2 & y^2/r^2\\\end{bmatrix}\nabla u\right)=\left(\nabla\cdot \begin{bmatrix}x^2/r^2 & xy/r^2\\xy/r^2 & y^2/r^2\\\end{bmatrix}\right)\cdot \nabla u + \begin{bmatrix}x^2/r^2 & xy/r^2\\xy/r^2 & y^2/r^2\\\end{bmatrix}\cdot\cdot\left(\nabla\nabla u\right)$$
in an attempt to create something equivalent to \eqref{4}?
Here are some related stack questions I looked into:
Laplacian from cartesian to polar
Laplacian of a Function depending on r in Polar Coordinates
How to transform a Laplacian operator from (x,y) coordinate system to polar system?
Transforming the Laplace operator from Polar to Cartesian coordinates
See also:
http://wwwf.imperial.ac.uk/~jdg/AECHAIN.PDF
http://planetmath.org/derivationofthelaplacianfromrectangulartosphericalcoordinates
