In this post we will run through regular white noise and iid noise and how to identify them in your residuals after applying for example a ARMA or GARCH model.

**In time series sample data there is always noise. But there exist different types of noise. In this post we will run through regular white noise and iid noise and how to identify them in your residuals after for example applying an ARMA or GARCH model. First we need to define weak stationarity before we define white noise.**

## Definition of Weak Stationarity

A time series $X$ is weakly stationary if

- $\mathbb{E}[X_t] = \mu = \text{const}$ for all $t \in \mathbb{Z}$
- $\text{Var}[X_t] = \sigma^2 = \text{const} < \infty$ for all $t \in \mathbb{Z}$
- $\text{Cov}[X_r, X_s] = \text{Cov}[X_{r+h}, X_{s+h}]$ for all $s,r,h \in \mathbb{Z}$

## Defintion of White Noise

A time series $X$ is said to be white noise if

- $X$ is weakly stationary
- $\gamma(h) = \text{Cov}[X_t, X_{t+h}] = 0$ where $h \neq 0$

## Definition of IID Noise

What is iid noise? A time series $X$ is said to be iid noise if it meets the conditions for regular white noise above but also $X_r$ and $X_s$ ($r \neq s$) is independent of each other and identically distributed for all $r,s \in \mathbb{Z}$, that is why it is called iid noise where iid stands for independent and identically distributed. 

Before going through the tricks to identifying the two different kind of noises we need to go through the definitions of the autocorrelation function (ACF) and the partial autocorrelation function (PACF).

## The Autocorrelation Function (ACF)

The autocorrelation function for the time series $X$ is defined as:

$$
    \rho_X(h) = \frac{\gamma(h)}{\gamma(0)} = \frac{\text{Cov}[X_t, X_{t+h}]}{\text{Var}[X_t]}.
$$

## The Partial Autocorrelation Function (ACF)

The autocorrelation function $\alpha(h)$ for the time series $X$ is defined by $\alpha(0) = 1$ and $\alpha(h) = \phi_{hh}$ for $h \geq 1$ where $\phi_{hh}$ is the last component of

$$
    \phi_h = ([\gamma(i-j)]_{i,j=1}^{h})^{-1} [\gamma(1), ... , \gamma(h)]'.
$$

## How to Identify White Noise

The first simple thing you could do to see if your data is just white noise is if it looks like it has no structure. The second thing is looking at the sample ACF and PACF graphs for the data. Lets look at the definition of white noise and the ACF. We know that $\gamma(h) = 0$ for $h \geq 1$ for white noise which means $\rho_X(h) = 0$ for $h \geq 1$. So all lags greater than zero should be within the red striped confidence band ($\pm 1.96/\sqrt{n}$, where $n$ is the sample size) of the ACF in the plot below. So the following sample ACF plot indicates white noise:

![Sample ACF of the sample data](https://i.imgur.com/3exyzTL.png "Sample ACF of the sample data")

The same can be said for the PACF. If $X$ is white noise then the vector $[\gamma(1), ... , \gamma(h)]'$ will just be filled with zeros and thus $\alpha(h) = 0$ for $h \geq 1$.

## How to Identify IID Noise

To identify iid noise we will have to look at the sample ACF and PACF of the data squared, $X^2$. Lets start with the ACF. We have:

$$
    \gamma_{X^2}(h) = \text{Cov}[X^2_t, X^2_{t+h}] = \mathbb{E}[(X_t^2 - \mathbb{E}[X_t^2])(X_{t+h}^2 - \mathbb{E}[X_{t+h}^2])]
$$

and now we use the property of iid noise that $X_t$ and $X_{t+h}$ is independent of each other so we get the following:

$$
    (\mathbb{E}[X_t^2 - \mathbb{E}[X_t^2]]) (\mathbb{E}[X_{t+h}^2 - \mathbb{E}[X_{t+h}^2]]) = (\mathbb{E}[X_t^2] - \mathbb{E}[\mathbb{E}[X_t^2]]) (\mathbb{E}[X_{t+h}^2] - \mathbb{E}[\mathbb{E}[X_{t+h}^2]]).
$$

Lastly we use the general property that $\mathbb{E}[\mathbb{E}[Y]] = \mathbb{E}[Y]$ so we get:

$$
    (\mathbb{E}[X_t^2] - \mathbb{E}[X_t^2]) (\mathbb{E}[X_{t+h}^2] - \mathbb{E}[X_{t+h}^2]) = 0.
$$

So the sample ACF for $X^2$ should be in the confidence band for all lags greater than zero. For the PACF we can see that the vector $[\gamma_{X^2}(1), ... , \gamma_{X^2}(h)]'$ would be zero for all $h$. Therefore the PACF should also be within the confidence band for lags greater than zero if $X$ were to be iid noise which we can see in the following plot:

![Sample PACF of the sample data squared](https://i.imgur.com/4ywfylN.png "Sample ACF of the sample data squared")

So to conclude, for iid noise $X$ and $X^2$'s sample ACF and PACF should all have to be within the confidence band for lags greater than zero.

## Ljung-Box Test for IDD Noise

Lastly one can do a Ljung-Box test to determine whether the data is independently distributed. The test is defined by the following hypothesis:

$$
        \begin{aligned}
            & H_0 : X \sim \text{IID}(\mu, \sigma^2) \\
            & H_1 : X \not\sim \text{IID}(\mu, \sigma^2) 
        \end{aligned}
$$

The Ljung-Box test has a test statistic which is:

$$
    \lambda = n (n + 2) \sum_{i=1}^h \frac{\hat{\rho}(i)^2}{n-i},
$$

where $\hat{\rho}(i)$ is the sample ACF at lag $i$, $n$ is the sample size and $h$ is the number of lags being tested. The null-hypothesis, $H_0$, is rejected at the $\alpha$-level, $\alpha \in (0,1)$, if $\lambda > \chi_{1-\alpha, h}^2$, where $\chi_{1-\alpha, h}^2$ denotes the $1 - \alpha$-quantile of the $\chi^2$-distribution with $h$ degrees of freedom.

Time Series White Noise

The continuity equation is an important tool that one can use whether it is in a fluid dynamics problem or an electromagnetic problem. The continuity equation for charges relates the charge density to the current density in a volume.

**The continuity equation is an important tool that one can use whether it is in a fluid dynamics problem or an electromagnetic problem. The continuity equation for charges relates the charge density to the current density in a volume.**

## Derivation of Continuity Equation for Charges

Lets derive the continuity equation for electric charges. If you are just searching for it, this is what it looks like:

$$
    \frac{\partial \rho}{\partial t} = - \nabla \cdot \vec{j},
$$

where $\rho$ is the charge density and $\vec{j}$ is the current density (current divided by an area). Lets get a few things straight first. When we say the charge is conserved we are talking about the charge in an arbitrary volume. Additional charge can not be created from some sort of source inside the volume. The only way to change the charge is a flow of new charges from outside the volume, that would be our $\vec{j}$. So lets observe an infinitesimal surface element $d \vec{S}$ on the surface of our volume. The surface element is vector that is normal to the surface. So the infinitesimal charge $d q$ that passes through the surface element can be written as:

$$
    dq = \vec{j} \cdot d \vec{S} dt.
$$

We'll use this expression later. The total charge $q$ for the whole volume $V$ can be written as:

$$
    q = \int_V \rho dV,
$$

where earlier mentioned $\rho$ is the charge density for the volume. We can now use the conservation of charge and look at the change of the charge through time and flux of charge through the volume's surface $\partial V$. So we have two expressions, firstly:

$$
    \frac{dq}{dt} = \int_V \frac{\partial \rho}{\partial t} dV,
$$

which is just the derivative of the total charge for the volume. Secondly, we have:

$$
    \frac{dq}{dt} = - \int_{\partial V} \vec{j} \cdot d \vec{S},
$$

which uses the expression for the charge through a surface element. The reason for the minus sign is because we are calculating the inflow of charges but the surface element vector is pointing outwards from the volume and we have to counteract that. We are almost at our sought after expression but we'll have to use one more trick, the divergence theorem aka Guass's theorem which states:

$$
    \int_V \nabla \cdot \vec{F} dV = \int_{\partial V} \vec{F} \cdot d \vec{S},
$$

which fits perfectly with our charge flux expression:

$$
     - \int_{\partial V} \vec{j} \cdot d \vec{S} = - \int_V \nabla \cdot \vec{j} dV.
$$

So now we have:

$$
    \frac{dq}{dt} = \int_V \frac{\partial \rho}{\partial t} dV = - \int_V \nabla \cdot \vec{j} dV.
$$

As you can see we have two expression with a volume integral for the volume $V$. So we can remove the integral and we finally get our continuity equation for electric charges:

$$
    \frac{\partial \rho}{\partial t} = - \nabla \cdot \vec{j}.
$$

Now lets see how the continuity equation looks in different coordinate systems.

## Continuity Equation in Cartesian Coordinates

The continuity equation for charges in Cartesian coordinates is given by:

$$
    \frac{\partial \rho}{\partial t} + \frac{\partial j_x}{\partial x} + \frac{\partial j_y}{\partial y} + \frac{\partial j_z}{\partial z} = 0.
$$

## Continuity Equation in Cylindrical Coordinates

The continuity equation for charges in cylindrical coordinates is given by:

$$
    \frac{\partial \rho}{\partial t} + \frac{1}{r} \frac{\partial}{\partial r}(r j_r) + \frac{1}{r} \frac{\partial}{\partial \theta}(j_\theta) + \frac{\partial}{\partial z}(j_z) = 0.
$$

## Continuity Equation in Spherical Coordinates

The continuity equation for charges in spherical coordinates is given by:

$$
    \frac{\partial \rho}{\partial t} + \frac{1}{r^2} \frac{\partial}{\partial r}(r^2 j_r) + \frac{1}{r \sin{\theta}} \frac{\partial}{\partial \theta}(j_\theta \sin{\theta}) + \frac{1}{r \sin{\theta}} \frac{\partial}{\partial \varphi}(j_\varphi) = 0.
$$

Continuity Equation in Electromagnetism

The moment of inertia is dependent on where the axis of rotation is located. Lets calculate the moment of inertia for the different rotational axes for a rod and a disk.

**If you are just searching for the moment of inertia of a rod or a disk I will just give them to you. For a rod $I = \frac{L^2}{12}$ if the rotation axis is located at the middle the rod and $I = M \frac{L^2}{3}$ if the rotational axis is located at the end of the rod. For a homogeneous disk/cylinder $I = M \frac{R^2}{2}$ and for a disk/cylinder with just a shell $I = MR^2$. Now if you are interested in how we get these values feel free to read further.**

We know that the moment of a point of inertia of a point mass is

$$
    M R^2,
$$

where $M$ is the mass and $R$ is the distance between the mass and the rotational axis. We can use this to calculate the moment of inertia for objects that is not point masses. We do this by integrating over all infinitesmal mass elements $dm$ of the object. So our moment of inertia is now

$$
    I = \int_0^M r^2 dm.
$$

## The Moment of Inertia of a Rod

The moment of inertia is dependent on the location of the rotational axis, therefore the moment of inertia will vary for a rod depending on the axis of rotation. We have our general formula for the moment of inertia

$$
    I = \int_0^M r^2 dm.
$$

We now have to identify what $dm$ is. Since the rod can be looked at in one dimension

$$
    dm = \frac{M}{L} dr
$$

where $L$ is the length of the rod. $\frac{M}{L}$ can be looked at as an one dimensional density for the rod. Now lets calculate the moment of inertia if the rod is rotating at one of its ends. The moment of inertia then becomes

$$
    I = \int_0^L r^2 \frac{M}{L} dr = \frac{M}{L} \Bigg[ \frac{r^3}{3} \Bigg]_0^L = \frac{M}{L} \frac{L^3}{3} = M \frac{L^2}{3}.
$$

If the axis of rotation is in the middle of the rod the moment of inertia then becomes

$$
    I = \int_{-L/2}^{L/2} r^2 \frac{M}{L} dr = \frac{M}{L} \Bigg[ \frac{r^3}{3} \Bigg]_{-L/2}^{L/2} = \frac{M}{3L} \left[ \frac{L^3}{8} - \left(\frac{-L^3}{8} \right) \right] = M \frac{L^2}{12}.
$$

Lets create a more general formula. Lets use the variable $x$ for the point of rotation. Our moment of inertia integral then becomes

$$
    I = \int_{-x}^{L-x} r^2 \frac{M}{L} dr = \frac{M}{L} \Bigg[ \frac{r^3}{3} \Bigg]_{-x}^{L-x} = \frac{M}{3L} \left[ (L-x)^3 + x^3 \right].
$$

## The Moment of Inertia of a Disk

![Image of the moment of inertia for disks](https://i.imgur.com/lyfEGNK.png "Moment of inertia for a homogeneous disk and a hollow disk")

Lets move on to calculate the moment of inertia for a disk (this also applies to a cylinder which is kind of an elongated disk). The axis of rotation will be located in the center of the disk but we will get back to calculating the moment of inertia of a displaced rotational axis. Like usual we start of with the classic integral for the moment of inertia

$$
    I = \int_0^M r^2 dm.
$$

Defining $dm$ is a little bit more tricky than for the rod. In this problem we have to go two dimensional. Lets start with defining our two dimensional density which is

$$
    \frac{M}{\pi R^2}.
$$

The denominator is just the area of the disk. Our area element is 

$$
    r dr d\theta
$$

I won't go into to much detail here, you'll have to do some repetition on polar coordinates if don't understand how we get this area element. So our moment of inertia becomes

$$
    I = \frac{M}{\pi R^2} \int_0^{2 \pi} \int_0^R r^3 dr d\theta = \frac{2M}{R^2} \int_0^R r^3 dr = \frac{2M}{R^2} \Bigg[ \frac{r^4}{4} \Bigg]_0^R = \frac{2M}{R^2} \frac{R^4}{4} = M \frac{R^2}{2}.
$$

Now, what happens if the axis of rotation was shifted (but still parallel to the axis through the center of mass)? Well luckily for us we don't have to do a super complicated integral, we can just use the parallel axis thereom which says

$$
    I = I_\text{CM} + Md^2,
$$

where $I_\text{CM}$ is the moment of inertia through the center of mass, in this case $M \frac{R^2}{2}$. $d$ is the distance between the axis of rotation and the rotational axis throught the center of mass. You can test this formula by calculating the moment of inertia at an end of the rod. In that case we have $I_\text{CM} = M \frac{L^2}{12}$ and $d=\frac{L}{2}$.

So what if the disk/cylinder is hollow? Then we simply change the limits of integration and the density. If the disk is hollow up to $r=a$ and then has mass up to $r=b$ the 2D density changes to

$$
    \frac{M}{\pi (b^2-a^2)}.
$$

So we get the following moment of inertia

$$
    I = \frac{2M}{b^2-a^2} \int_a^b r^3 dr = \frac{2M}{b^2-a^2} \Bigg[ \frac{r^4}{4} \Bigg]_a^b = \frac{2M}{b^2-a^2} \frac{b^4-a^4}{4} = M \frac{b^2+a^2}{2}.
$$

Here we can see that if the disk/cylinder only consisted of a shell ($a \rightarrow b$) we get

$$
    I = M b^2.
$$

So now you know how to calculate the moment of inertia of a rod and a disk.

The Moment of Inertia of a Rod and a Disk

What makes the suspension wire of the Golden Gate Bridge curve like a parabola? Lets get into the mechanics of the shape of the suspension bridge's cables and what determines them to curve like a parabola or a catenary.

**What makes the suspension wire of the Golden Gate Bridge curve like a parabola? Lets get into the mechanics of the shape of the suspension bridge's cables and what determines them to curve like a parabola or a catenary.**

## The Assumptions about the Problem

Some assumptions about the problem have to be made first. The main wire that the bridge rests on is assumed to be inextensible, homogeneous and infinitely flexible. The bridge is also assumed to be homogeneous. The main wire is connected at each end of the bridge. The wires that ties together the bridge and the main wire is assumed to be in an infinite amount so the bridge's load function is constant. The origin of coordinate system is located at the bottom of the main wire.

## Equilibrium of an Infinitesimal Wire Element

Lets make an free body diagram of an infinitesmal wire element of the main wire:

![Image of free body diagram of an infinitesmal wire element](https://i.imgur.com/vGGtqIv.png "Free body diagram of an infinitesmal wire element")

Here $\vec{T}$ is the tension of the cable, $d\vec{T}$ the infinitesmal change in the tension, $ds$ the infinitesmal length of the wire element and $dx$ the length in the $x$-direction. $q_c$ and $q_b$ is respectively the cable's and bridge's linear mass density. The equilibrium of the forces is

$$
    \vec{T} + d\vec{T} - \vec{T} - (q_c ds + q_b dx)g \hat{y} = 0.
$$

With some cancellation we get

$$
    d\vec{T} = (q_c ds + q_b dx)g \hat{y}.
$$

Lastly, the Pythagorean theorem gives us 

$$
    ds = \sqrt{(dx)^2 + (dy)^2} = \sqrt{1 + (y')^2} dx.
$$

So we get

$$
    d\vec{T} = (q_c \sqrt{1 + (y')^2} dx + q_b dx)g \hat{y}.
$$

## The ODE which Dictates the Shape of The Wire

The tension at the bottom of the main wire has only a $x$-component which we'll call $T_0$. So $\vec{T}(0) = T_0 \hat{x}$. If we do some integration of the equilibrium equation from earlier we get

$$
    \vec{T} = T_0 \hat{x} + \left(q_c \int_0^x \sqrt{1 + y'(\tilde{x})^2} d\tilde{x} + q_b x \right)g \hat{y},
$$

here we are just using $\tilde{x}$ in the integral so the upper limit, $x$, of the integral is correct. In the $x$-direction we now have

$$
    T \cos{\theta} = T \frac{dx}{\sqrt{(dx)^2 + (dy)^2}} = T_0,
$$

where $\theta$ is the angle between $\vec{T}$ and the $x$-axis. So

$$
    T = T_0 \frac{\sqrt{(dx)^2 + (dy)^2}}{dx}.
$$

In the $y$-direction we have

$$
    T \sin{\theta} = T \frac{dy}{\sqrt{(dx)^2 + (dy)^2}} = T_0 \frac{dy}{dx} = \left(q_c \int_0^x \sqrt{1 + y'(\tilde{x})^2} d\tilde{x} + q_b x \right)g.
$$

We now have an expression for $y'$:

$$
    y' = \left(q_c \int_0^x \sqrt{1 + y'(\tilde{x})^2} d\tilde{x} + q_b x \right) \frac{g}{T_0}
$$

Lets do an additional derivation to get rid of the integral:

$$
    y'' = \frac{q_b g}{T_0} + \frac{q_c g \sqrt{1+(y')^2}}{T_0}.
$$

This is the differential equation that dictates the curve of suspension bridge.

## Is the Cable Shaped as a Parabola or Catenary?

So what shape is the cable? To answer this lets look at two extremes, neglecting either the cable's mass or the bridge's mass. If we neglect the cable's load function i.e. its mass we get the differential equation

$$
    y'' = \frac{q_b g}{T_0}.
$$

With some simple integration we get

$$
    \frac{q_b g}{2 T_0} x^2 + bx + c,
$$

where $b$ and $c$ are constants. The origin is loacted at the bottom of the cable so $c=0$ and the cable is symmetric so $b=0$. Here we have our parabola. Neglecting the cable's mass is a pretty reasonable approximation since the bridge's mass is likely much larger than the cable's mass. Can you guess what we get if we do the opposite? If we neglect the bridge's mass we get the following differential equation:

$$
    y'' = \frac{q_c g \sqrt{1+(y')^2}}{T_0}.
$$

Since we already have pretty good hunch that the solution is in the shape of a catenary, lets check if $y = a \cosh{\frac{x}{a}}$ works. Then the left side is

$$
    y'' = \frac{\cosh{\frac{x}{a}}}{a}
$$

and the right side is

$$
    \frac{q_c g \sqrt{1+(y')^2}}{T_0} = \frac{q_c g \sqrt{1+\sinh^2{\frac{x}{a}}}}{T_0} =  \frac{q_c g \cosh{\frac{x}{a}}}{T_0}.
$$

For this to work $a$ needs to be

$$
    a = \frac{T_0}{q_c g}.
$$

That is the solution if the bridge's mass can be neglected. I don't know what sort of bridge that is... But hey now you know. I'll leave the harder task of solving the ODE with both masses present for you.

Is the Suspension Bridge Shaped like a Parabola or a Catenary?

When studying transient processes, like when a switch is changed, in electrical circuits we have to make us of the Laplace transform. Here is a quick rundown on how to Laplace transform the resistor, capacitor and inductor.

**When studying transient processes, like when a switch is changed, in electrical circuits we have to make us of the Laplace transform. Here is a quick rundown on how to Laplace transform the resistor, capacitor and inductor.**

First we need to establish what happens when we Laplace transform the current and voltage. We can see the voltage and current as "unknown" functions dependent on time, therefore their Laplace transforms become:

$$
        \begin{aligned}
            & i(t) \Rightarrow I(s) \\
            & u(t) \Rightarrow U(s)
        \end{aligned}
        .
$$

## Laplace Transform of the Resistor

So lets start with the easy component: the resistor. We know from Ohm's law that:

$$
    u(t) = R i(t).
$$

If we Laplace transform this we get:

$$
    U(s) = R I(s),
$$

since the resistor is time independent there is not much happening when moving from the t-domain to the s-domain. The resistor remains the same in both domains.

## Laplace Transform of a Derivative

The inductor and capacitor is time dependent, so they get a little bit more trickier. First we need to see what happens when we apply the Laplace transform to a derivative of a function $f(t)$. The Laplace transform of an arbitrary function $f(t)$ is given by:

$$
    F(s) = \int_{0-}^\infty e^{-st} f(t) dt.
$$

The $0-$ can be looked at more or less as what happened exactly before significant change happened to the function $f(t)$ (we'll get back to this). So lets preform the Laplace transform on $f'(t)$:

$$
    \int_{0-}^\infty e^{-st} f'(t) dt.
$$

With a little partial integration we get:

$$
    \bigg[e^{-st} f(t) \bigg]_{0-}^\infty - \int_{0-}^\infty -se^{-st} f(t) dt = - f(0-) + sF(s) .
$$

## Laplace transform of the Inductor

With that out of the way lets preform the Laplace transform on the inductor first. We know that the potential over an inductor is given by:

$$
    u(t) = L \frac{d i(t)}{dt}.
$$

So with our knowledge of Laplace transforming a derivative we get the following:

$$
    U(s) = L (s I(s) - i(0-)).
$$

The inductor is said to be energized if there is an initial current $i(0-)$. If there is no initial current we get a simpler relation of:

$$
    U(s) = sL I(s).
$$

We can also rewrite the relation to get the current:

$$
    I(s) = \frac{U(s)}{sL} + \frac{i(0-)}{s}.
$$

We can make circuit diagrams with these relations:

![Image of Laplace transform of the inductor](https://i.imgur.com/IkgvQUe.png "Laplace transform of the inductor")

Note the Thevenin/Norton equivalents under "Initial current".

## Laplace transform of the Capacitor

Lets move on to the capacitor. The procedure of Laplace transforming the capacitor is the same as the inductor. The current through a capacitor is given by:

$$
    i(t) = C \frac{du(t)}{dt}.
$$

So the Laplace transform of the current is:

$$
    I(s) = C(sU(s) - u(0-)).
$$

We can see here that if there is no initial voltage over the capacitor we get the simpler relation:

$$
    I(s) = sC U(s).
$$

Lastly the the Laplace transformed voltage is:

$$
    U(s) = \frac{I(s)}{sC} + \frac{u(0-)}{s}.
$$

The circuit diagrams then become:

![Image of Laplace transform of the capacitor](https://i.imgur.com/za9EeTS.png "Laplace transform of the capacitor")

And that is it! You now know how to Laplace transform the resistor, inductor and capacitor!

Laplace Transform of DC Circuits

The Navier Stokes equation is famously known for being unsolvable if there is any turbulence in a fluid but there are certain scenarios where laminar flow exist and a solution for the velocity profile can be derived. One of these scenarios is flow past two plates either caused by one plate moving, Couette flow, or a pressure gradient.

**The Navier Stokes equation is famously known for being unsolvable if there is any turbulence in a fluid but there are certain scenarios where laminar flow exist and a solution for the velocity profile can be derived. One of these scenarios is flow past two plates either caused by one plate moving, Couette flow, or a pressure gradient.**

Before solving these two problems we need to make some assumptions. First we assume that our fluid is incompressible and Newtonian so we can use the Navier Stokes equation and the mass continuity equation. Secondly, the flow has to be stationary (so all time dependent terms in the Navier Stokes equation is zero) and in one direction, the $x$-direction, so the velocity vector only has one component ($v=w=0$). The plates also has to be considered infinite. Lastly, the flow has to be laminar which has already been mentioned and the impact from gravity has to be considered negligible ($g = 0$).

## Couette Flow - One Plate Moves

So lets start with the case where one of the plates is moving. From our assumptions that fluid is Newtonian, stationary and goes in the $x$-direction we get the following mass continuity equation:

$$
    \frac{\partial u}{\partial x} + \cancel{\frac{\partial v}{\partial y}} + \cancel{\frac{\partial w}{\partial z}} = \frac{\partial u}{\partial x} = 0,
$$

where $u$, $v$ and $w$ is the velocity in the $x$-, $y$- and $z$-direction. It is then safe to assume that the velocity profile $u$ is only dependent on $y$ because the plates' areas are infinite. Now lets look at the Navier Stoke equation in all directions and see which terms we can cancel out because of our assumptions and the mass contuinity equation:

$$
    \left\{
        \begin{aligned}
            & \rho \left( \cancel{\frac{\partial u}{\partial t}} + \cancel{u \frac{\partial u}{\partial x}} + \cancel{v \frac{\partial u}{\partial y}} + \cancel{w \frac{\partial u}{\partial z}} \right) = \cancel{- \frac{\partial p}{\partial x}} + \mu \left(\cancel{\frac{\partial^2 u}{\partial x^2}} + \frac{\partial^2 u}{\partial y^2} + \cancel{\frac{\partial^2 u}{\partial z^2}} \right) + \cancel{\rho g_x} \\
            & \rho \left( \cancel{\frac{\partial v}{\partial t}} + \cancel{u \frac{\partial v}{ \partial x}} + \cancel{v \frac{\partial v}{\partial y}} + \cancel{w \frac{\partial v}{\partial z}} \right) = \cancel{- \frac{\partial p}{\partial y}} + \mu \left(\cancel{\frac{\partial^2 v}{\partial x^2}} + \cancel{\frac{\partial^2 v}{\partial y^2}} + \cancel{\frac{\partial^2 v}{\partial z^2}} \right) + \cancel{\rho g_y} \\
            & \rho \left( \cancel{\frac{\partial w}{\partial t}} + \cancel{u \frac{\partial w}{ \partial x}} + \cancel{v \frac{\partial w}{\partial y}} + \cancel{w \frac{\partial w}{\partial z}} \right) = \cancel{- \frac{\partial p}{\partial z}} + \mu \left(\cancel{\frac{\partial^2 w}{\partial x^2}} + \cancel{\frac{\partial^2 w}{\partial y^2}} + \cancel{\frac{\partial^2 w}{\partial z^2}} \right) + \cancel{\rho g_z}
        \end{aligned}
    \right. .
$$

As you can see, a bunch of terms cancel out and all that remains is just one single expression:

$$
    \mu \frac{\partial^2 u}{\partial y^2} = 0 \Rightarrow \frac{\partial^2 u}{\partial y^2} = 0.
$$

So if we integrate this expression twice we get:

$$
    u(y) = C_1 y + C_2,
$$

where $C_1$ and $C_2$ are constants. To determine these constants we have to use the no-slip boundary conditions of the plates. If we assume that the origin is located exactly in between the two plates we know that:

$$
    \left\{
        \begin{aligned}
            & u(+h) = u_p \\
            & u(-h) = 0
        \end{aligned}
    \right. ,
$$

where $h$ is the shortest distance from the origin to a plate and $u_p$ is the speed of the moving plate. This gives the constants:

$$
    \left\{
        \begin{aligned}
            & C_1 = \frac{u_p}{2 h} \\
            & C_2 = \frac{u_p}{2}
        \end{aligned}
    \right. .
$$

So our velocity profile becomes:

$$
    u(y) = \frac{u_p}{2} \left(1 + \frac{y}{h} \right).
$$

This solution is called the Couette flow which is named after the french physicist Maurice Couette. 

## Laminar Flow Past Two Plates Caused by Pressure Gradient

So what solution do we get if the two plates are stationary and we instead have a constant pressure gradient? The assumptions mentioned earlier are still the same so we get the same mass continuity equation. So lets see what cancels out in our Navier Stokes equations:

$$
    \left\{
        \begin{aligned}
            & \rho \left( \cancel{\frac{\partial u}{\partial t}} + \cancel{u \frac{\partial u}{\partial x}} + \cancel{v \frac{\partial u}{\partial y}} + \cancel{w \frac{\partial u}{\partial z}} \right) = - \frac{\partial p}{\partial x} + \mu \left(\cancel{\frac{\partial^2 u}{\partial x^2}} + \frac{\partial^2 u}{\partial y^2} + \cancel{\frac{\partial^2 u}{\partial z^2}} \right) + \cancel{\rho g_x} \\
            & \rho \left( \cancel{\frac{\partial v}{\partial t}} + \cancel{u \frac{\partial v}{ \partial x}} + \cancel{v \frac{\partial v}{\partial y}} + \cancel{w \frac{\partial v}{\partial z}} \right) = - \frac{\partial p}{\partial y} + \mu \left(\cancel{\frac{\partial^2 v}{\partial x^2}} + \cancel{\frac{\partial^2 v}{\partial y^2}} + \cancel{\frac{\partial^2 v}{\partial z^2}} \right) + \cancel{\rho g_y} \\
            & \rho \left( \cancel{\frac{\partial w}{\partial t}} + \cancel{u \frac{\partial w}{ \partial x}} + \cancel{v \frac{\partial w}{\partial y}} + \cancel{w \frac{\partial w}{\partial z}} \right) = - \frac{\partial p}{\partial z} + \mu \left(\cancel{\frac{\partial^2 w}{\partial x^2}} + \cancel{\frac{\partial^2 w}{\partial y^2}} + \cancel{\frac{\partial^2 w}{\partial z^2}} \right) + \cancel{\rho g_z}
        \end{aligned}
    \right. .
$$

We can see that the pressure gradient is zero in the $y$- and $z$-direction. So we get the following differential equation:

$$
    \mu \frac{\partial^2 u}{\partial y^2} = \frac{\partial p}{\partial x}.
$$

With some integration we get:

$$
    u(y) = \frac{1}{\mu} \frac{\partial p}{\partial x} \frac{y^2}{2} + C_1 y + C_2,
$$

where $C_1$ and $C_2$ are constants. The velocity at both plates are now zero because of the no-slip condition, so:

$$
    \left\{
        \begin{aligned}
            & u(+h) = 0 \\
            & u(-h) = 0
        \end{aligned}
    \right. .
$$

This gives the constants:

$$
    \left\{
        \begin{aligned}
            & C_1 = 0 \\
            & C_2 = \frac{1}{\mu} \frac{\partial p}{\partial x} \frac{h^2}{2}
        \end{aligned}
    \right. .
$$

So our velocity profile becomes:

$$
    u(y) = \frac{1}{\mu} \frac{\partial p}{\partial x} \frac{h^2}{2} \left(1 + \left(\frac{y}{h} \right)^2 \right).
$$

The velocity profile is a parabola with its maximum at the middle of the plates. 

![Image of pressure gradient flow](https://i.imgur.com/xIiH6ge.png "Pressure Gradient Flow")

So that is it, that is how you calculate fluid flow past two plates!

Couette Flow - Laminar Flow Past Two Plates

Why is it that a falling chimney breaks at a third of its length? Lets find out why!

**The falling chimney problem is the perfect mix of mechanical dynamics and beam theory. By just using the fundamentals of those two fields we can accurately predict where the chimney cracks when falling, under some assumptions...**

So lets get the assumptions of the falling chimney problem out of the way first. The first assumption is that there is no air resistance. If air resistance had been present the top of the chimney would experience a greater stress from the air than the lower part which makes our calculations much more complex and this is not an article about fluid mechanics. The second assumption is that the chimney is homogeneous, meaning that it has the same diameter across its length and the same material. If it had been inhomogeneous certain parts would handle stress better than others.

## Angle Acceleration of the Chimney

Now when that is out of the way, lets make some calculations! The first thing we'll have to determine is the angle acceleration $\ddot{\theta}$. To do that we use a beam's, in this case a chimney's, moment of inertia when attached at one of its ends:

$$
    I = \frac{1}{3} m L^2.
$$

Here $m$ is the chimneys mass and $L$ is its length. The moment from the gravitation acting on the chimney must equal that of the angle acceleration:

$$
    mg \frac{L}{2} \sin{\theta} = I \ddot{\theta} = \frac{1}{3}mL^2 \ddot{\theta}.
$$

So now we have our angle acceleration $\ddot{\theta}$ and apparently it is dependent on $\theta$:

$$
    \ddot{\theta} = \frac{3g}{2L} \sin{\theta}.
$$

Is this plausible? Well $\ddot{\theta}$ increases with $g$ and decreases with $L$ which makes sense. At $\theta=0$ the angle acceleration is zero which also makes sense since at that angle the chimney is "standing up".

## The Chimney's Load Function

If the chimney was locked and stationary at an angle $\theta$ the load function would be:

$$
    w(d) = -\frac{mg}{L} \sin{\theta}.
$$

If we go to the other extreme and say that the chimney was in free fall it would lead to the load function being zero. So falling negates the load function and in our case the chimney is rotating. So $\ddot{\theta}$ has a similar effect to free fall. So the load function in our case will be the following:

$$
    w(d)=\frac{dm}{L} \ddot{\theta} - \frac{mg}{L} \sin{\theta} = \frac{3dmg}{2L^2} \sin{\theta} - \frac{mg}{L} \sin{\theta}.
$$

We can see here that the load function is zero at $d = \frac{2L}{3}$. This makes sense because $\ddot{\theta}$ is equal to $g \sin{\theta}$ at $d = \frac{2L}{3}$. In other words, this point is in a "virtual" free fall.

## The Chimney's Shear Force Function

The relation between the load function and the shear force function is:

$$
    -w(d)=V'(d).
$$

If we integrate that we get:

$$
    V(d) = -\frac{3d^2mg}{4L^2} \sin{\theta}+ \frac{dmg}{L} \sin{\theta}+ C,
$$

where $C$ is a constant. We can determine this constant with the boundary condition $V(L) = 0$. In other words: the shear force is zero at the end of the chimney which is logical. This means:

$$
    C = -\frac{mg}{4}.
$$

If we take the absolute value of the shear force function we can check where the force is the greatest:

$$
    |V(d)|=\left|mg\, \sin{\theta} \left(-\frac{3d^2}{4L^2} + \frac{d}{L} - \frac{1}{4}\right) \right|.
$$

Here we can see that the shear force is the greatest when $d \to 0$ on the interval $0< d < L$, so if the shear force was the cause of the crack in the chimney it would be at the base. But that is not the only place where the chimney can break... 

## The Chimney's Bending Moment Function

The final step is looking at the bending moment function which has the following relation with the shear force function:

$$
    M'(d)=V(d).
$$

So with a bit of integration we get:

$$
    M(d)=mg\, \sin{\theta} \left(-\frac{d^3}{4L^2} + \frac{d^2}{2L} - \frac{d}{4}\right) + K,
$$

where $K$ is a constant. We'll once again use the boundary condition at the end of the chimney, which in this case is also zero. This leads to $K=0$. So where is $|M(d)|$ the greatest value? At $d=0$ the bending moment is zero so it is not there. Well lets check where its derivative, the shear force function, is zero to get the maximums:

$$
    M'(d)=V(d)=0 \Longrightarrow -\frac{3d^2}{4L^2} + \frac{d}{L} - \frac{1}{4} = 0.
$$

From this we get the roots $d=L$ and $d=\frac{L}{3}$. $d=L$ does not make much sense since its the end of the chimney, you cant real chop of something that isn't there... So we get $d=\frac{L}{3}$ and that is why a falling chimney breaks at a third of its length from its base.

<div style="position: relative;
    padding-bottom: 56.25%; /* 16:9 */
    height: 0;">
<iframe style="position: absolute;
    top: 0;
    left: 0;
    width: 100%;
    height: 100%;" width="560" height="315" src="https://www.youtube.com/embed/efjegZ6aWKI?start=149" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>

</div>

The video above shows a demolition of a smokestack, as you can see it breaks close to a third of its length when it falls.

The Falling Chimney

The chemical potential has different characteristics depending on the temperature. Lets calculate and understand the approximations of the chemical potential in an electron doped semiconductor.

**The chemical potential has different characteristics depending on the temperature. Lets calculate and understand the approximations of the chemical potential in an electron doped semiconductor.**

## The Chemical Potential at Low Temperatures

Lets start with temperatures close to absolute zero. You have probably heard or read that the chemical potential must lie in the center between the donor level and the bottom of the conduction band because of charge neutrality. Lets do some math to find out why that is.

Lets use the variable $E_g$ for the bandgap energy, $E_d$ is the donor energy so the the donor level is situated at $E_g-E_d$. $\mu$ is our desired chemical potential. $n$ is the electron density in the conduction band, $n_d$ is the donor density and $n_d^+$ is the density of the ionized donors. From the Fermi-Dirac distribution we know that the probability of an electron at the energy $E_g-E_d$ and the temperature $T$ is:

$$
    f(E_g - E_d) = \frac{1}{e^{(E_g-E_d-\mu)/(k_B T)}+1}.
$$

So the donor electron density at the donor level is:

$$
    \frac{n_d}{e^{(E_g-E_d-\mu)/(k_B T)}+1}.
$$

Every donor electron not at the donor level is at the conduction band so:

$$
    n_d^+ = n_d (1 - \frac{1}{e^{(E_g-E_d-\mu)/(k_B T)}+1}) = \frac{n_d}{e^{-(E_g-E_d-\mu)/(k_B T)}+1}.
$$

Since we are at low temperatures $T \ll E_d/k_B$. This means $-(E_g-E_d-\mu)/(k_B T) \gg 1$ so we can make the following approximation (the 1 in the denominator is neglected):

$$
    n_d^+ \approx n_d e^{(E_g-E_d-\mu)/(k_B T)}.
$$

Having low temperatures also means that the electrons at the valence band don\'t have enough energy to make an excitation to the conduction band meaning $n \approx n_d^+$. But the electron density $n$ can also be expressed as:

$$
    n = \frac{1}{\sqrt{2}} (\frac{m_e^* k_B T}{\pi \hbar^2})^{3/2} e^{-(E_g-\mu)/(k_B T)} = n_0 e^{-(E_g-\mu)/(k_B T)}
$$

here is $m_e^*$ the effective mass of the electrons in the conduction band. Lets put the different expressions for $n$ at low temperature equal to each other to get the chemical potential:

$$
    n_d e^{(E_g-E_d-\mu)/(k_B T)} = n_0 e^{-(E_g-\mu)/(k_B T)}
$$

So now we get the chemical potential:

$$
    \mu = E_g - \frac{1}{2} E_d - \frac{1}{2} k_B T \ln{(\frac{n_0}{n_d})}.
$$

So when the temperature approaches zero, $T \to 0$, the chemical potential becomes:

$$
    \mu = E_g - \frac{1}{2}E_d,
$$

since polynomial functions dominates over logarithmic functions ($n_0$ contains $T^{3/2}$) when it comes to limits.

## The Chemical Potential at Medium Temperatures

So that is how the chemical potential behaves at low temperatures, how about "medium" temperatures like room temperature? Well here we have more thermal energy to work with so pretty much every donor atom is ionized ($n_d^+ \approx n_d$) but the thermal energy is still not big enough to excite a significant amount of electrons in the valence band so $n \approx n_d$. This makes the math easier, we now have:

$$
    n = n_0e^{-(E_g-\mu)/(k_B T)} = n_d,
$$

so:

$$
    \mu = E_g + k_B T \ln{(\frac{n_d}{n_0})},
$$

for "medium" temperatures.

## The Chemical Potential at High Temperatures

For high temperatures, like 1200 K, electron holes has manifested at the valence band and now plays the bigger role in the electron density in the conduction band. So now $n \approx p$, where $p$ is the electron hole density at the valence band. The semiconductor is said to be in an intrinsic state, meaning that it can be seen as undoped. So now we have the following relation:

$$
    \frac{1}{\sqrt{2}} (\frac{m_e^* k_B T}{\pi \hbar^2})^{3/2} e^{-(E_g-\mu)/(k_B T)} = \frac{1}{\sqrt{2}} (\frac{m_h^* k_B T}{\pi \hbar^2})^{3/2} e^{-\mu/(k_B T)},
$$

where $m_h^*$ is the effective mass of the electron holes. This leads to the chemical potential of:

$$
    \mu = \frac{1}{2} E_g + \frac{3}{4} k_B T \ln{(\frac{m_h^*}{m_e^*})}.
$$

The chemical potential is located midway between the conduction band and the valence band if $m_e^* = m_h^*$.

So there we have it! We have now calculated the chemical potential for different temperature regions for an n-type semiconductor, although we did it with some significant approximations but without those this would get a lot harder if not impossible.

The Chemical Potential of an n-type Semiconductor

Calculating the energy of the stationary state of a two electron system, like a helium atom or lithium ion, gets quiet a bit more tricky than for a single electron system. In fact, calculating the Schrödinger equation is impossible. Luckily we have methods like perturbation theory and the variational method to estimate the energy.

**Calculating the energy of the stationary state of a two electron system, like a helium atom or lithium ion, gets quiet a bit more tricky than for a single electron system. In fact, calculating the Schrödinger equation is impossible. Luckily we have methods like perturbation theory and the variational method to estimate the energy.**

First we need to establish the Hamiltonian for the system which is:
$$
    \hat{H} = -\frac{\hbar}{2 m_e} \nabla^2_1 - \frac{\hbar}{2 m_e} \nabla^2_2 - \frac{Ze^2}{4 \pi \epsilon_0 r_1 } - \frac{Ze^2}{4 \pi \epsilon_0 r_2 } + \frac{e^2}{4 \pi \epsilon_0 | \vec{r}_1 - \vec{r}_2 | }.
$$
$\frac{\hbar}{2 m_e} \nabla^2_1$ is the kinetic energy for the first electron, same goes for the second. $\frac{Ze^2}{4 \pi \epsilon_0 r_1 }$ is Coulomb's law between the atom core's charge $Ze$ and the first electron's charge $e$. Lastly we have $\frac{e^2}{4 \pi \epsilon_0 | \vec{r}_1 - \vec{r}_2 | }$ which is Coulomb's law between the two electrons, this term is positive because the electrons repel each other.

## Two Independent Electrons Estimation

The first estimation we'll do is having the two electrons not react with each other which means that the Coulomb's law term between them is zero. Then the combined wave function of the system is the following:

$$
    \psi_+ (\vec{r}_1, \vec{r}_2) = \psi_{100}(r_1) \psi_{100}(r_2),
$$

where

$$
    \psi_{100}(r) = \frac{Z^{3/2}}{\sqrt{\pi} a_0^{3/2}} e^{Zr/a_0}
$$

is the wave function of a single electron system, $a_0$ is the Bohr radius. So our wave function of the two electron system is

$$
    \psi_+ (r_1, r_2) = \frac{Z^{3}}{\sqrt{\pi} a_0^{3}} e^{Z(r_1+r_2)/a_0},
$$

which leads the stationary state energy of:

$$
    E = 2Z^2 \cdot (-13.6 \text{ eV}).
$$

## Perturbation Theory Estimation

But that is not so satisfactory nor a good estimate, which leads to our second estimation method: perturbation theory. Lets use the first order correction to estimate the impact that Coulomb's law between the two electrons has on the energy:

$$
    E^{(1)} = \langle \psi_+ | \frac{e^2}{4 \pi \epsilon_0} \frac{1}{|\vec{r}_1 - \vec{r}_2|} | \psi_+ \rangle = \frac{5}{4} Z \cdot 13.6 \text{ eV}.
$$

So now our estimated stationary state energy is:

$$
    E = 2Z^2 \cdot (-13.6 \text{ eV}) + \frac{5}{4} Z \cdot 13.6 \text{ eV}.
$$

## Variational Method Estimation

Lets try getting an even better estimation with the variational method. If we change the atom number $Z$ to a parameter $\alpha$ we get:

$$
    \psi_0 (r_1, r_2) = \frac{\alpha^{3}}{\sqrt{\pi} a_0^{3}} e^{\alpha (r_1+r_2)/a_0}.
$$

We now use this wave function as the base to calculate the expectation value of the Hamiltonian:

$$
    E(\alpha) = \langle \psi_0 (r_1, r_2) | \hat{H} | \psi_0 (r_1, r_2) \rangle = (-2 \alpha^2 + 4 Z \alpha - \frac{5}{4} \alpha) \cdot (-13.6 \text{ eV}).
$$

Here we can see that if we choose $\alpha = Z$ we get the same result as for the first order correction perturbation theory. The variational method says $E_\text{ground} \leq E(\alpha)$ so lets check where the minimum of $E(\alpha)$ is:

$$
    \frac{\partial E(\alpha)}{\partial \alpha} = (-4 \alpha + 4Z - \frac{5}{4}) = 0 \Longrightarrow \alpha = Z - \frac{5}{16}.
$$

So our estimated stationary state energy with the variational method is

$$
    E_\text{ground} = (Z-\frac{5}{16})^2 \cdot (-13.6 \text{ eV}).
$$

For example, the estimated stationary state energy with the variational method for helium ($Z=2$) is $-77.5$ eV which is close to the experimental value of $-78.975$ eV, in other words this method of estimation works very well! You now know a couple of methods to estimate the stationary state energy for electrons!

FILTERED POSTS

Time Series White Noise

Continuity Equation in Electromagnetism

The Moment of Inertia of a Rod and a Disk

Is the Suspension Bridge Shaped like a Parabola or a Catenary?

Laplace Transform of DC Circuits

RECENT POSTS