Monday, February 27, 2017

Derivation of Divergence in Cylindrical Coordinates

Derivation of Divergence in Cylindrical Coordinates $$ \newcommand{\pard}[2]{\frac{\partial #1}{\partial #2}} $$ $$ \newcommand{\vec}[1]{\boldsymbol{#1}} $$ $$ \newcommand{\unitv}[1]{\boldsymbol{\hat{#1}}} $$ The divergence of a vector field in Cartesian coordinates is fairly straight forward: $$ \nabla\cdot \vec{F} = \pard{F_x}{x} + \pard{F_y}{y} + \pard{F_z}{z} $$ But before we undertake this derivation, it will help to try and understand the meaning of divergence. It is tempting to just take the partial derivatives in whatever coordinate system you're using and add them together. But alas, this doesn't quite work.

First, let us consider flux.

Flux is the amount of something passing through a surface multiplied by the area of that surface. It can be a fluid, gas, electric field, or anything that can be described by a vector field. In the picture below, there are F liters/minute of water flowing in the circular pipe.

If we were to consider an imaginary surface slicing perpendicularly through the pipe (surface A), the flux through that circular surface would be F liters/minute multiplied by the area of the circle. It is important to note that that if we were to consider a surface slicing through at a bias (like surface B), the flux would be the same. Even though the elliptical area is larger than the circular one the amount of flux per unit area is less because the flux is no longer perpendicular to the surface. Flux through an area is then defined as $$ \Psi = \vec{F}\cdot\vec{A} = \|\vec{F}\|\|\vec{A}\|\cos\theta $$ where the vector area \( \vec{A} \) is defined to be normal to the surface. For a surface that is open, the choice of which of the two normal vectors we mean is arbitrary (we just have to be consistent). For a closed surface, the normal is defined to point to the outside of the surface. The angle \(\theta\) is the angle between the area-normal and the vector \(\vec{F}\). You can see as the angle between the flow and the surface gets closer to \(\frac{\pi}{2}\) that the flux gets smaller and smaller until it vanishes completely when the two are perpendicular.

Now, consider the volume defined by surfaces A, B and the walls of the pipe in the figure above. We can define the total flux through the volume as $$ \Psi_{total} = \Psi_{pipe} + \Psi_A + \Psi_{B} $$ Unless we have a leaky pipe, \(\Psi_{pipe} = 0\). Because we define the surface normal of a closed surface to point outward, the flux through surface A is negative (the flow and the normal are in opposite directions, so the dot-product is negative). The flow through B is equal in magnitude to that passing through surface B but positive. Hence, \(\Psi_{total} = 0 \). This makes intuitive sense when dealing with conservation of mass as otherwise there would be an accumulation of mass inside the pipe (\(\Psi_{total} < 0\)) or loss of mass (\(\Psi_{total} > 0\)).

Divergence is the amount of flux leaving a closed surface as the volume of the closed surface approaces zero. So, we now look at the general sitution of the pipe above in cylindrical coordinates. We consider a small enough closed volume that we can consider it a cube. In order to compute the total flux, we have to do 6 integrations as there are 6 sides to a cube. But we will group them in pairs, one pair for each direction \(\rho, \phi, z\).

The vector field in cylindrical coordinates is given by $$ \vec{F} = F_{\rho} \unitv{\rho} + F_{\phi} \unitv{\phi} + F_{z} \unitv{z} $$ We place our tiny cube-like volume in the field \(\vec{F}\) at some point \(\vec{P}=(\rho,\phi,z)\).

First we consider the flux passing through the surfaces normal to the \(z\) direction (the top and bottom faces of the cube). Since the distances involved are small, we can approximate the field passing through these as $$ F_{top} = F_z(\rho,\phi,z + \Delta z) \approx F_z(\rho,\phi,z) + \pard{F_z}{z}(\rho,\phi,z)\frac{\Delta z}{2} $$ It should be noted that the other components of the field don't play a part here because they are parallel to these surfaces so nothing going in the, say, \(\rho\) direction leaves via the \(z\) direction. Similarly for the bottom surface, $$ F_{bottom} = F_z(\rho,\phi,z - \Delta z) \approx F_z(\rho,\phi,z) - \pard{F_z}{z}(\rho,\phi,z)\frac{\Delta z}{2} $$ The areas of these surfaces are equal and independent of \(z\). $$ A_z = \rho\Delta\phi\Delta\rho $$ Therefore, the flux through these z-surfaces is (remembering that the flux through the bottom surface is into the volume and so is negative) $$ \Psi_z = F_{top}A_z - F_{bottom}A_z = (F_{top} - F_{bottom})A_z $$ $$ \Psi_z = (F_z + \pard{z}{F_z}\frac{\Delta z}{2} - (F_z - \pard{F_z}{z}\frac{\Delta z}{2}))\rho\Delta\phi\Delta\rho $$ $$ \Psi_z = \pard{F_z}{z}\rho\Delta \phi \Delta\rho \Delta z $$ Recognize that the term \( \rho\Delta \phi \Delta\rho \Delta z \) is the differential volume in cylindrical coordinates, so we can write the previous equation as $$ \begin{equation} \label{eq:fluxz} \Psi_z = \pard{F_z}{z}\Delta V \end{equation} $$ Where \(\Delta V\) is the differential volume element. I should note that if you were to do this derivation in Cartesian coordinates, the x and y component derivations would be identical (except for the subscripts). But in cylindrical coordinates, things get a little more interesting from here on out. Looking now at the \(\phi\) direction we consider the front and back faces of the volume. (I will omit the function arguments to the partial derivatives in the service of brevity with the understanding that we are talking about the point \(P=(\rho,\phi,z)\)). $$ F_{back} = F_\phi(\phi + \Delta \phi) \approx F_\phi + \pard{F_\phi}{\phi}\frac{\Delta \phi}{2} $$ $$ F_{front} = F_\phi(\phi - \Delta \phi) \approx F_\phi - \pard{F_\phi}{\phi}\frac{\Delta \phi}{2} $$ $$ A_\phi = \Delta\rho\Delta z $$ Because the positive \(\phi\) direction is toward the back of the volume, $$ \Psi_{\phi} = (F_{back} - F_{front})A_{\phi} $$ $$ \Psi_{\phi} = (F_\phi + \pard{F_\phi}{\phi}\frac{\Delta \phi}{2} - (F_\phi - \pard{F_\phi}{\phi}\frac{\Delta \phi}{2}))\Delta\rho\Delta z $$ $$ \Psi_{\phi} = \pard{F_\phi}{\phi}\Delta\phi\Delta\rho\Delta z $$ The term \(\Delta\phi\Delta\rho\Delta z\) is almost the differential volume. And maybe because we know where we are headed, we bend the equation to our will $$ \Psi_{\phi} = \pard{F_\phi}{\phi}\Delta\phi\Delta\rho\Delta z = \pard{F_\phi}{\phi}\frac{\rho}{\rho}\Delta\phi\Delta\rho\Delta z $$ $$ \begin{equation} \label{eq:fluxphi} \Psi_\phi = \frac{1}{\rho}\pard{F_\phi}{\phi}\Delta V \end{equation} $$ Now we try to handle the \(\rho\) component. We need to be careful here. But first things first. The vector field is defined similarly to previous components. $$ F_{right} = F_\rho(\rho + \Delta \rho) \approx F_\rho + \pard{F_\rho}{\rho}\frac{\Delta \rho}{2} $$ $$ F_{left} = F_\phi(\rho - \Delta \rho) \approx F_\rho - \pard{F_\rho}{\rho}\frac{\Delta \rho}{2} $$ Here is where it gets tricky. The area definitions are different because the two surfaces are not equal in size and so must be treated differently. $$ A_{right} = (\rho+\frac{\Delta\rho}{2})\Delta\phi\Delta z $$ $$ A_{left} = (\rho-\frac{\Delta\rho}{2})\Delta\phi\Delta z $$ First computing the flux on the right side. $$ \Psi_{right} = F_{right}A_{right}=(F_\rho + \pard{F_\rho}{\rho}\frac{\Delta \rho}{2})((\rho+\frac{\Delta\rho}{2})\Delta\phi\Delta z) $$ $$ \Psi_{right} = (F_\rho + \pard{F_\rho}{\rho}\frac{\Delta \rho}{2})(\rho\Delta\phi\Delta z + \frac{1}{2}\Delta\rho\Delta\phi\Delta z) $$ $$ \Psi_{right} = F_\rho\rho\Delta\phi\Delta z+\frac{1}{2}F_\rho\Delta\rho\Delta\phi\Delta z + \frac{1}{2}\pard{F_\rho}{\rho}\rho\Delta\rho\Delta\phi\Delta z + \frac{1}{4}\pard{F_\rho}{\rho}\Delta\rho^2\Delta\phi\Delta z $$ Similarly on the left hand side $$ \Psi_{left} = F_{left}A_{left}=(F_\rho - \pard{F_\rho}{\rho}\frac{\Delta \rho}{2})((\rho-\frac{\Delta\rho}{2})\Delta\phi\Delta z) $$ $$ \Psi_{left} = (F_\rho - \pard{F_\rho}{\rho}\frac{\Delta \rho}{2})(\rho\Delta\phi\Delta z - \frac{1}{2}\Delta\rho\Delta\phi\Delta z) $$ $$ \Psi_{left} = F_\rho\rho\Delta\phi\Delta z-\frac{1}{2}F_\rho\Delta\rho\Delta\phi\Delta z - \frac{1}{2}\pard{F_\rho}{\rho}\rho\Delta\rho\Delta\phi\Delta z + \frac{1}{4}\pard{F_\rho}{\rho}\Delta\rho^2\Delta\phi\Delta z $$ Computing \(\Psi_\rho\) (finally, after all that algebra) by subtracting and cancelling terms $$ \Psi_\rho = \Psi_{right} - \Psi_{left} $$ $$ \Psi_\rho = F_\rho\Delta\rho\Delta\phi\Delta z + \pard{F_\rho}{\rho}\rho\Delta\rho\Delta\phi\Delta z $$ We make the same modification to the first term that we did with \(\Psi_\phi\) and substituting \(\Delta V\) $$ \begin{equation} \label{eq:fluxrho} \Psi_\rho = \frac{1}{\rho}F_\rho\Delta V + \pard{F_\rho}{\rho}\Delta V \end{equation} $$ This last equation is correct but we normally don't see it written this way, so we factor out a \(\frac{1}{\rho}\) and invoke the reverse chain rule. $$ \Psi_\rho = \frac{1}{\rho}(F_\rho + \rho\pard{F_\rho}{\rho})\Delta V = \frac{1}{\rho}\pard{}{\rho}(\rho F_\rho)\Delta V $$ And so the total flux through our tiny volume is $$ \Psi_{total} = \Psi_\rho + \Psi_\phi + \Psi_z = \left(\frac{1}{\rho}\pard{}{\rho}(\rho F_\rho) + \frac{1}{\rho}\pard{F_\phi}{\phi} + \pard{F_z}{z}\right)\Delta V $$ The divergence is the ratio of the flux through a closed surface as the volume of that closed surface approaces 0. $$ \nabla\cdot\vec{F} = \lim\limits_{\Delta V \to 0} \frac{\Psi_{total}}{\Delta V} = \frac{1}{\rho}\pard{}{\rho}(\rho F_\rho) + \frac{1}{\rho}\pard{F_\phi}{\phi} + \pard{F_z}{z} $$

The derivation of divergence in spherical coordinates procedes along similar lines as the \(\phi\) and \(\rho\) terms. It's a little more complicated because of the added \(\sin\theta\) terms.

Tuesday, February 21, 2017

Geodesic on a Plane in Polar Coordinates

Geodesic on a Plane in Polar Coordinates This is a simple problem in the calculus of variations where we prove that the shortest distance between two points is a straight line. But, instead of the normal Cartesian coordinate system, we will use polar coordinates \( \left( r, \phi \right) \). Some people use \( \theta \) for the angle but to be consistent with the cylindrical and spherical coordinate nomenclature, I choose to use \( \phi \). I was somewhat stumped at the end because there is a difficult integral to evaluate but we'll cross that bridge when we get to it.

In typical fashion, we wish to minimize the functional $$ I = \int_a^b ds $$ And the first thing we need to do is figure out what \( ds \) is. Looking at the figure below,

we can see that the differential length of the curve is the hypotenuse of the right triangle (well, not quite a triangle but it approaches one as \(d\phi \rightarrow 0\)) and is given by $$ ds = \sqrt{dr^2 + (rd\phi)^2} $$ In order to get something we can integrate, we have a decision to make: we need to decide which differential to factor out of the radical and integrate by either \(dr\) or \(d\phi\). Do we go with $$ ds = \sqrt{\left(\frac{dr}{d\phi}\right)^2 + r^2} d\phi $$ or $$ ds = \sqrt{1 + (r\frac{d\phi}{dr})^2} dr $$ In order to decide, we need to look at the Euler-Lagrange equation knowing that the \(F\) in the equation is the integrand in our functional. $$ F = F(x,y,y') : \frac{d}{dx}\left(\frac{\partial F}{\partial y'}\right) - \frac{\partial F}{\partial y} = 0 $$ The equation becomes significantly simpler to solve if we can choose an integrand that is not a function of the middle argument, i.e. a function only of the independent variable and the derivative. If this is the case, then \(\frac{\partial F}{\partial y} = 0 \) leaving the first term to be 0 as well. If the first term is 0, then inner partial derivative must be a constant because the derivative of a constant is 0. $$ \frac{\partial F}{\partial y'} = K $$ In our case,that leads us to use the second equation $$ ds = \sqrt{1 + (r\frac{d\phi}{dr})^2} dr = \sqrt{1 + (r\phi')^2} dr $$ because the first form contains a dependency on \(r\) as well as \(r'\).

So we now have to solve $$ F = F(r,\phi,\phi') : \frac{\partial F}{\partial \phi'} = \frac{\partial}{\partial \phi'}\left(\sqrt{1 + (r\phi')^2}\right) = K $$ Taking the partial derivative with respect to \(\phi'\), we have $$ {{r^2\phi'} \over \sqrt{1 + r^2\phi'^2}} = K $$ Squaring both sides gives $$ {{r^4\phi'^2} \over {1 + r^2\phi'^2}} = K^2 $$ Then $$ {r^4\phi'^2} = K^2 \left( 1 + r^2\phi'^2 \right) $$ and solving for \(\phi'2\) gives us $$ \phi'^2 = \frac{K^2}{r^2\left( r^2 - K^2 \right)} $$ Taking the square root and moving \( dr \) to the right side, gives $$ d\phi = \frac{K}{r\sqrt{ r^2 - K^2 }}dr $$ So, here is where it got tricky for me. If you type the above integral into Wolfram-Alpha you'll get something ridiculously complicated. If you're smarter than me, you might be able to simplify it but I had to look at this integral for a while before I came up with an answer. The trick with this stuff is knowing what spell to cast.

Whenever you see square-roots with squares and quotients, think right-triangles.

You can see that the \( \frac{K}{\sqrt{r^2 - K^2}} \) term is \(\tan \alpha\). So, we make that substitution giving us $$ d\phi = \frac{\tan \alpha}{r} dr $$ Now, what do we do about the \( r \)? Well, looking at the triangle, $$ \sin \alpha = \frac{K}{r} $$ This also means that $$ r = K\csc \alpha $$ so $$ dr = -K\cot(\alpha) \csc(\alpha) d\alpha $$ Making those substitutions leaves us with $$ d\phi = \frac{1}{K}\frac{K}{r}\tan(\alpha)dr = -\frac{1}{K}\sin(\alpha)\tan(\alpha)K\cot(\alpha)\csc(\alpha)d\alpha $$ And lo and behold, the entire thing becomes $$ d\phi = -d\alpha $$ which is something even I can solve. $$ \phi = -\alpha + C $$ Solving for \(\alpha\) in terms of \(r\), $$ \sin(\alpha) = \frac{K}{r} \Rightarrow \alpha = \sin^{-1}(\frac{K}{r}) $$ Plugging this back into our result gives us $$ \phi = -\sin^{-1}\left(\frac{K}{r}\right) + C $$ $$ \sin(C-\phi) = \frac{K}{r} $$ $$ r\sin(C-\phi) = K $$ If you don't recognize the above as a line in polar coordinates, we can go one step further using the sin-of-the-sum-of-angles identity $$ r\left(\sin(C)\cos(\phi)-\cos(C)\sin(\phi)\right) = K $$ and gathering up constants, we are left with $$ Ar\cos(\phi) + Br\sin(\phi) = K $$ and since \(x = r\cos(\phi), y = r\sin(\phi)\) are the transformations for polar to Cartesian coordinates $$ Ax + By = K $$ which is the general form of a line in Cartesian coordinates.

It's always satisfying to see something work out.