In [4]:
import numpy as np
import matplotlib.pyplot as plt
import random
import math

import ipywidgets
from ipywidgets import widgets
from IPython.display import display

Boundary Conditions

  • Dirichlet: $u(x,t)$ is given on the boundary
  • Neumann: $\frac{du}{dt}$ is given on the boundary
  • Robin: Linear combination of Dirichlet and Neumann (Ex: $\alpha u+\beta \frac{du}{dt}=c \,\text{on}\,\partial \Omega$)
  • Periodic: $u(x,t)$ is periodic in the spacial domain

Note: $\Delta$ is the Laplacian and $\Delta = \nabla^2$ and $\nabla = \begin{pmatrix}\frac{\partial}{\partial x}\\ \frac{\partial}{\partial y}\\ \frac{\partial}{\partial z}\end{pmatrix}$.

Brownian Motion and the Heat Equation

$$\frac{du}{dt}=D\Delta u$$ Models a Brownian motion. What is the mean squared displacement of a particle?
$$var = 2Dt$$

Normal Probability Distribution: $$f(x)=\frac{1}{2\sigma^2\pi}e^{\frac{-(x-\mu)^2}{2\sigma^2}}$$ Given a Heat Equation starting from $\delta_0$, $$u(x,t)=\frac{1}{\sqrt{4Dt\pi}}e^{\frac{-(x-0)^2}{4Dt}}$$

Einstein's Derivation of the Heat Equation

Numerical Solution to the Heat Equation

Approximate the second derivative in space with: $$\frac{\partial u}{\partial x} = \lim_{h\to 0} \frac{u(x+h,t)-u(x,t)}{h} = \lim_{h\to 0}\frac{u(x,t)-u(x-h,t)}{h}$$ We can apply these definitions onto each other to find, $$\frac{\partial^2 u}{\partial x^2} = \lim_{h\to 0} \frac{\frac{u(x+h,t)-u(x,t)}{h} - \frac{u(x,t)-u(x-h,t)}{h}}{h}$$ Choosing a small but nonzero $h$ we see: $$\frac{\partial^2 u}{\partial x^2} \approx \frac{u(x+h)-2u(x,t)+u(x-h,t)}{h^2}$$

Moving to a numerical solution, $$u(x,t)\approx \vec{u}(t) = \begin{pmatrix}u_1(t)\\u_2(t)\\\vdots\\u_n(t)\end{pmatrix}$$ $$\frac{d\vec{u}}{dt} = \frac{D}{h^2}\text{A}\vec{u}$$ With derivative matrix A: $$A = \begin{pmatrix} &&&\vec{b_0}\\ 1 & -2 & 1 & 0 & 0 & \dots & 0\\ 0 & 1 & -2 & 1 & 0 & \dots & 0\\ \vdots\\ 0 & 0 & \dots & 0 & 1 & -2 & 1\\ &&&\vec{b_1} \end{pmatrix}$$ The row vectors $\vec{b_0}$ and $\vec{b_1}$ control the boundary conditions. For Dirchlet, assume they are zero vectors. For Periodic, $\vec{b_0} = \begin{pmatrix} -2 & 1 & 0 & \dots & 1\end{pmatrix}$ and $\vec{b_1} = \begin{pmatrix} 1 & 0 & \dots & 1 & -2\end{pmatrix}$
So, what we're doing is solving $$\frac{\partial u}{\partial t} = \Delta u$$ $$\frac{u(x, t+h)-u(x,t)}{h} = Au(x,t)$$ $$u(x, t+h)= hAu(x,t) + u(x,t)$$ Which means we can calculate the spacial spread at time $t+h$ using the spacial second derivative matrix and the spacial spread at time $t$.

The Heat Equation

In [6]:
from enum import Enum
#size = 11
#boundary = np.concatenate([np.mat([-2,1]), np.zeros([1, size - 3])])
#print (boundary)

class BoundaryConditions(Enum):
    periodic = 1
    dirchlet = 0

def getLaplacianMatrix(boundary, size):  #may not be stable -- be careful, not sure I understand laplacian transform matrix
    A = np.zeros([size,size])
    A += -2*np.identity(size)
    A += np.diagonal(np.ones([size-1,1]),  1)
    A += np.diagonal(np.ones([size,1]), -1)
    if(boundary == BoundaryConditions.dirchlet):
        A[0] = np.zeros([size, 1])
        A[size-1] = np.zeros([size, 1])
    elif(boundary == BoundaryConditions.periodic):
        tmp = A[0]
        A[0] = A[size-1]
        A[size-1] = t
    return A

def run1DHeatSim(initx, dt, steps, boundary):
    r = [initx]
    A = dt*getLaplacianMatrix(boundary, np.size(initx))
    for i in range(0, steps):
    return r

def run(size, dt, steps, boundary):
    initx = np.zeros([size,1])
    initx[round(size/2)] = 1
    solution = run1DHeatSim(initx, dt, steps, boundary)

#run(5, 1, 20, BoundaryConditions.dirchlet)

Code up a 1D Random Walk with Drift

Lattice is $\in[-100, 100]\cup Z$
$p(\text{left}) = 0.45$
$p(\text{right}) = 0.55$
Walk for 80 steps

In [96]:

seeds = [874083802,29305632,94086608,-60502738,32466780,-8312131,72793714,35637054,53822097,-10185695,-94233408]
boundmin = -100
boundmax = 100

def walk1D(bias=0.55):
    pos = 0
    history = [pos]
    t = [0]
    steps = 100
    boundmin = -10
    boundmax = 10
    for i in range(0,steps):
        pos = step(pos, bias, 1)
    plt.plot(t, history, "#84B12F")

def roll(bias):
    return (random.random() < bias)

def step(ipos, bias, size):
    if(ipos == boundmin):
            return boundmax - size
            return ipos - size
    if(ipos == boundmax):
            return ipos + size
            return boundmin + size
        return ipos + size
    return ipos - 1

def update(x):

widgets.interact(update, x=0.55)
<function __main__.update>
In [41]:
def walk2D(xpos=0, ypos=0, steps=10000, seed=4, boundmin=-100, boundmax=100):
    xhistory = [xpos]
    yhistory = [ypos]

    for i in range(0,steps):
        stp = random.random()*2
        ypos = step(ypos, 0.5, stp)
        xpos = step(xpos, 0.5, stp)

    plt.plot(yhistory, xhistory, "#84B12F")
    plt.xlabel("x position")
    plt.ylabel("y position")

_steps = 1000000
_seed = 4
_boundmin = -100
_boundmax = 100

def func(btn):
    walk2D(steps=_steps, seed=_seed, boundmin=_boundmin, boundmax=_boundmax)

def fsteps(x):
    _steps = x

def fseed(x):
    _seed = x

def fmin(x):
    _boundmin = x

def fmax(x):
    _boundmax = x

widgets.interact(fsteps, x=1000000)
widgets.interact(fseed, x=4)
widgets.interact(fmin, x=-100)
widgets.interact(fmax, x=100)

btn = widgets.Button(description="Run")

Navier-Stokes Equation

In physics, the Navier–Stokes equations, named after Claude-Louis Navier and George Gabriel Stokes, describe the motion of viscous fluid substances. These balance equations arise from applying Newton's second law to fluid motion, together with the assumption that the stress in the fluid is the sum of a diffusing viscous term (proportional to the gradient of velocity) and a pressure term—hence describing viscous flow. The main difference between them and the simpler Euler equations for inviscid flow is that Navier–Stokes equations also in the Froude limit (no external field) are not conservation equations, but rather a dissipative system, in the sense that they cannot be put into the quasilinear homogeneous form:

$$y_t + A(y)y_x=0$$


$$\frac{\partial u}{\partial t} = D\frac{\partial^2 u}{\partial x^2} - a \frac{\partial u}{\partial x}$$ This can be interpreted as something spreading or diffusing in the second derivative term and something travelling in the first derivative term.


We begin with two assumptions, that mass and momentum are both conserved (that is: we are in an intertial frame of reference and our fluid body neither gains nor looses material). Consider a volume $\chi(\omega, t)$

From our assumption of conservation of mass: $$\frac{\partial m}{\partial t} = \frac{d}{dt}\int_{\chi(\omega, t)}\rho(\vec{x},t) = 0$$

And from conservation of momentum: $$p = \int_{\chi(\omega,t)}\rho\vec{u}d\vec{x}$$ $$\frac{dp}{dt} = \frac{d}{dt}\int_{\chi(\omega,t)}\rho\vec{u}d\vec{x} = $$ Which by the transport theorem is $$ = \int_{\chi(\omega,t)} \frac{\partial}{\partial t}(\rho\vec{u}) + \nabla\cdot(\vec{u}\otimes\vec{u}\rho)d\vec{x}$$

$$\int_{\chi(\omega,t)}\rho\left(\frac{\partial \vec{u}}{\partial t} +\vec{u}\cdot\nabla\vec{u}\right)d\vec{x} = \int_{\chi(\omega,t)} g(\omega, \vec{x}) d\vec{x}$$

Where $g(\omega, \vec{x})$ is the function of all forces acting on the fluid.

Continued Derivation

What are the forces that define $g(\omega, \vec{x})$? Consider the forces on the surface, $$\int_{\chi}\rho\frac{D\vec{u}}{Dt}d\vec{x}=\int_{\partial\chi}\sigma\cdot\hat{n}dS$$ With unit normal $\hat{n}$ and stress tensor $\sigma$.

(note: $\text{stress}=\frac{\text{force}}{\text{area}}$)

Using the divergence theorem, $$\int_{\chi}\rho\frac{D\vec{u}}{Dt}d\vec{x}=\int_{\chi}\nabla\cdot\sigma\vec{x}$$ The stress tensor can be built by $$\sigma = -PI+\tau$$ With $-P$ as a scalar for pressure and $\tau$ as a matrix of the viscous forces.

So $$\nabla\cdot\sigma = \nabla\cdot(-PI) + \nabla\cdot\tau$$ $$\nabla\cdot(-PI) = \begin{pmatrix}\frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z}\end{pmatrix} \cdot \begin{pmatrix}-P&0&0\\0&-P&0\\0&0&-P\end{pmatrix} = \begin{pmatrix}-\frac{\partial P}{\partial x}\\-\frac{\partial P}{\partial y}\\-\frac{-\partial P}{\partial z}\end{pmatrix}$$

For Newtonian Fluids, we can write the viscous forces as $$\tau=\mu\left(\nabla\vec{u}+(\nabla\vec{u})^T\right)$$

Extending the viscous forces, $$\tau=\mu\left(\nabla\vec{u}+(\nabla\vec{u})^T\right) = $$ $$ = \mu\begin{pmatrix} \frac{\partial \vec{u_1}}{\partial x} & \frac{\partial \vec{u_1}}{\partial y} & \frac{\partial \vec{u_1}}{\partial z}\\ \frac{\partial \vec{u_2}}{\partial x} & \frac{\partial \vec{u_2}}{\partial y} & \frac{\partial \vec{u_2}}{\partial z}\\ \frac{\partial \vec{u_3}}{\partial x} & \frac{\partial \vec{u_3}}{\partial y} & \frac{\partial \vec{u_3}}{\partial z} \end{pmatrix} +\mu \begin{pmatrix} \frac{\partial \vec{u_1}}{\partial x} & \frac{\partial \vec{u_2}}{\partial x} & \frac{\partial \vec{u_3}}{\partial x}\\ \frac{\partial \vec{u_1}}{\partial y} & \frac{\partial \vec{u_2}}{\partial y} & \frac{\partial \vec{u_3}}{\partial y}\\ \frac{\partial \vec{u_1}}{\partial z} & \frac{\partial \vec{u_2}}{\partial z} & \frac{\partial \vec{u_3}}{\partial z} \end{pmatrix} =$$ $$ = \begin{pmatrix} \mu\left(\frac{\partial \vec{u_1}}{\partial x}+\frac{\partial \vec{u_1}}{\partial x} \right) & \mu\left(\frac{\partial \vec{u_1}}{\partial y}+\frac{\partial \vec{u_2}}{\partial x} \right) & \mu\left(\frac{\partial \vec{u_1}}{\partial z}+\frac{\partial \vec{u_3}}{\partial x} \right) \\ \mu\left(\frac{\partial \vec{u_2}}{\partial x}+\frac{\partial \vec{u_1}}{\partial y} \right) & \mu\left( \frac{\partial \vec{u_2}}{\partial y}+\frac{\partial \vec{u_2}}{\partial y} \right) & \mu\left(\frac{\partial \vec{u_2}}{\partial z}+\frac{\partial \vec{u_3}}{\partial y} \right) \\ \mu\left(\frac{\partial \vec{u_1}}{\partial z}+\frac{\partial \vec{u_3}}{\partial x} \right) & \mu\left(\frac{\partial \vec{u_3}}{\partial y}+\frac{\partial \vec{u_2}}{\partial z} \right) & \mu\left(\frac{\partial \vec{u_3}}{\partial z}+\frac{\partial \vec{u_3}}{\partial z} \right) \end{pmatrix}$$

Writing this matrix as

$$\mu \begin{pmatrix} \tau_{xx} & \tau_{xy} & \tau_{xz} \\ \tau_{yx} & \tau_{yy} & \tau_{yz} \\ \tau_{zx} & \tau_{zy} & \tau_{zz} \end{pmatrix}$$

We can interpret $\tau_{nm}$ as the stress acting on the $n$ face from the $m$ direction.

Finally, taking the divergence of $\tau$ $$(\nabla\cdot\tau)= $$

So, finally, $$\rho\left(\frac{\partial\vec{u}}{\partial t} + \vec{u}\cdot\nabla\vec{u}\right)=-\nabla P + \mu\Delta\vec{u}$$

Definitions (from Wikipedia)

Nondimensionalization is the partial or full removal of units from an equation involving physical quantities by a suitable substitution of variables. This technique can simplify and parameterize problems where measured units are involved. It is closely related to dimensional analysis. In some physical systems, the term scaling is used interchangeably with nondimensionalization, in order to suggest that certain quantities are better measured relative to some appropriate unit. These units refer to quantities intrinsic to the system, rather than units such as SI units. Nondimensionalization is not the same as converting extensive quantities in an equation to intensive quantities, since the latter procedure results in variables that still carry units.

Nondimensionalization can also recover characteristic properties of a system. For example, if a system has an intrinsic resonance frequency, length, or time constant, nondimensionalization can recover these values. The technique is especially useful for systems that can be described by differential equations. One important use is in the analysis of control systems. One of the simplest characteristic units is the doubling time of a system experiencing exponential growth, or conversely the half-life of a system experiencing exponential decay; a more natural pair of characteristic units is mean age/mean lifetime, which correspond to base e rather than base 2.

Many illustrative examples of nondimensionalization originate from simplifying differential equations. This is because a large body of physical problems can be formulated in terms of differential equations.

Although nondimensionalization is well adapted for these problems, it is not restricted to them. An example of a non-differential-equation application is dimensional analysis; another example is normalization in statistics.

Measuring devices are practical examples of nondimensionalization occurring in everyday life. Measuring devices are calibrated relative to some known unit. Subsequent measurements are made relative to this standard. Then, the absolute value of the measurement is recovered by scaling with respect to the standard.

Reynolds Number

In fluid mechanics, the Reynolds number (Re) is a dimensionless quantity that is used to help predict similar flow patterns in different fluid flow situations. The concept was introduced by George Gabriel Stokes in 1851,[2] but the Reynolds number is named after Osborne Reynolds (1842–1912), who popularized its use in 1883.[3][4]

The Reynolds number is defined as the ratio of inertial forces to viscous forces and consequently quantifies the relative importance of these two types of forces for given flow conditions.[5] Reynolds numbers frequently arise when performing scaling of fluid dynamics problems, and as such can be used to determine dynamic similitude between two different cases of fluid flow. They are also used to characterize different flow regimes within a similar fluid, such as laminar or turbulent flow:

laminar flow occurs at low Reynolds numbers, where viscous forces are dominant, and is characterized by smooth, constant fluid motion; turbulent flow occurs at high Reynolds numbers and is dominated by inertial forces, which tend to produce chaotic eddies, vortices and other flow instabilities. In practice, matching the Reynolds number is not on its own sufficient to guarantee similitude. Fluid flow is generally chaotic, and very small changes to shape and surface roughness can result in very different flows. Nevertheless, Reynolds numbers are a very important guide and are widely used.

Reynolds number interpretation has been extended into the area of arbitrary complex systems as well: financial flows,[6] nonlinear networks,[7] etc. In the latter case an artificial viscosity is reduced to nonlinear mechanism of energy distribution in complex network media. Reynolds number then represents a basic control parameter which expresses a balance between injected and dissipated energy flows for open boundary system. It has been shown [7] that Reynolds critical regime separates two types of phase space motion: accelerator (attractor) and decelerator. High Reynolds number leads to a chaotic regime transition only in frame of strange attractor model.

The Reynolds number can be defined for several different situations where a fluid is in relative motion to a surface.[n 1] These definitions generally include the fluid properties of density and viscosity, plus a velocity and a characteristic length or characteristic dimension. This dimension is a matter of convention – for example radius and diameter are equally valid to describe spheres or circles, but one is chosen by convention. For aircraft or ships, the length or width can be used. For flow in a pipe or a sphere moving in a fluid the internal diameter is generally used today. Other shapes such as rectangular pipes or non-spherical objects have an equivalent diameter defined. For fluids of variable density such as compressible gases or fluids of variable viscosity such as non-Newtonian fluids, special rules apply. The velocity may also be a matter of convention in some circumstances, notably stirred vessels. The Reynolds number is defined as:

Nondimensionalization of a Mass-Spring System

$$m\frac{d^2x}{dt^2} + b\frac{dx}{dt} + kx = 0$$ x is in meters, set $L\chi=x$ such that $\chi$ is dimensionless.
t is in seconds, set $T\tau=t$ such that $\tau$ is dimensionless. $$f=f(t(\tau))$$ $$\frac{df}{dt} = \frac{d\tau}{dt}\frac{df}{d\tau}$$ So, $$\frac{d}{dt}=\frac{1}{T}\frac{d}{d\tau}$$ And, $$\frac{d^2}{dt^2}=\left(\frac{1}{T}\frac{d}{d\tau}\right)^2=\frac{1}{T^2}\frac{d^2}{d\tau^2}$$ With these, we can find the system $$m\frac{1}{T^2}\frac{d^2}{d\tau^2}(L\chi) +b\frac{1}{T}\frac{d}{d\tau}(L\chi)+k(L\chi)=0$$ Dividing by $\frac{mL}{T^2}$, $$\frac{d^2\chi}{d\tau}+\frac{bT}{m}\frac{d\chi}{d\tau}+\frac{k}{m}T^2\chi=0$$ Now we can set $T=\sqrt{\frac{m}{k}}$ so that $$\frac{d^2\chi}{d\tau^2}+\frac{b}{\sqrt{mk}}\frac{d\chi}{d\tau}+\chi=0$$ Finally, just to make the equation look nice, we can define $\zeta=\frac{b}{2\sqrt{mk}}$ so, $$\frac{d^2\chi}{d\tau^2}+2\zeta\frac{d\chi}{d\tau}+\chi=0$$ This equation can be solved as $$\left(\alpha^2+2\zeta\alpha+1\right)e^{\alpha\tau}=0$$ And, $$\alpha=-\zeta\pm\sqrt{\zeta^2-1}$$ So, $$\chi(\tau)=c_1e^{\left(-\zeta+\sqrt{\zeta^2-1}\right)\tau}+c_2e^{\left(-\zeta-\sqrt{\zeta^2-1}\right)\tau}$$

Nondimensionalization of Navier-Stokes

Examining again the incompressible Naiver-Stokes: $$\rho\left(\frac{d\vec{u}}{dt}+\vec{u}\cdot\nabla\vec{u}\right)=-\nabla P+\mu\Delta\vec{u}$$ As with the spring system, we first separate variables from their units: $$\vec{x}=L\vec{\chi}\quad\vec{u}=V\vec{U}\quad P=\rho V^2Q$$ Also, time can be defined with simply $$t=\frac{L}{V}\tau$$ Now we use the chain rule $$\frac{d}{dt}=\frac{d}{d\tau}\frac{d\tau}{dt}=\frac{V}{L}\frac{d}{d\tau}$$ And, $$\frac{\partial}{\partial z}=\frac{1}{L}\frac{\partial}{\partial Z}$$ $$\tilde{\nabla}=\left(\frac{\partial}{\partial X},\frac{\partial}{\partial Y}, \frac{\partial}{\partial Z}\right)$$ $$\nabla = \frac{1}{L}\tilde{\nabla}$$ $$\Delta=\frac{1}{L^2}\tilde{\nabla}^2$$ So, plugging these into the original equation: $$\rho\frac{V^2}{L}\left(\frac{\partial\vec{U}}{\partial\tau}+\vec{U}\cdot\tilde{\nabla}\vec{U}\right)=-\rho\frac{V^2}{L}\tilde{\nabla}^2Q+\frac{\mu V}{L^2}\tilde{\Delta}\vec{U}$$ $$\frac{\partial\vec{U}}{\partial\tau}+\vec{U}\cdot\tilde{\nabla}\vec{U}=-\tilde{\nabla}Q+\frac{\mu}{\rho VL}\tilde{\Delta}\vec{U}$$ The last coefficient on the right hand side is our nondimensionalized Reynolds number inverted so $$Re=\frac{\rho VL}{\mu}$$ Further, this can be broken down into the inertial term (the numerator, ?) and the viscous term (the viscosity, $Pa\cdot s$ or $\frac{kg}{s\cdot m}$)

Interpretaion of Reynold's Number


Problem 4, ManyWalk2D

In [86]:
#function many_walk2D(N, M);

#Simulates many random walkers that have an equally likely chance of moving
#left, right, up or down.
#N is the number of time steps
#M is the number of walkers

N = 20
M = 3

x = np.zeros([N,1]) #initialize the vector to hold the x-positions
y = np.zeros([N,1]) #initialize the vector to hold the y-positions
x[0] = 0     #set initial x-position of walker
y[0] = 0     #set initial y-position of walker

#set parameters here
stepsize = 1 #stepsize of walker
bins = (N/80) #size of bins
binmax = (N/2) #maximum bin
binnum = round(binmax/bins) #number of bins

rad = np.zeros([M,1]) #final distance of walker
nn = np.arange(0, bins, binmax)
size1 = np.size(nn)

for j in range(1, M): #loop through all walkers
    for i in range(2, N): #loop through all of the time steps
        r = random.random() #pick a random number to determine up/down or right/left
        if r > 0.5: #go right/left
            r1 = random.random() #pick a random number to determine right or left
            if r1 > 0.5: #go right
                y[i] = y[i-1]
                x[i] = x[i-1] + stepsize    
            else: #go left
                y[i] = y[i-1]
                x[i] = x[i-1] - stepsize
        else: #go up/down
            r2 = random.random() #pick a random number to go up or down
            if r2 > 0.5: #go up
                x[i] = x[i-1]
                y[i] = y[i-1] + stepsize;    
            else: #go down
                x[i] = x[i-1]
                y[i] = y[i-1] - stepsize;

    rad[j] = math.sqrt(x[N-1]**2+y[N-1]**2) #This holds the final distance from the starting point of each walker
    #use these to set the axes
    xmax = max(abs(x))+1;
    ymax = max(abs(y))+1;

    plt.plot(x, y, 'r'); #plot the path
    plt.axis([-N/2, N/2, -N/2, N/2]);
    plt.title('Two-dimensional ' + str(N) + '-step random walk')
    #zoom on

#calculate theoretical distribution
x1 = np.zeros([100, 1])
pr = np.zeros([100, 1])
dx = nn[size1-1]/100

for k in range(0, 99):
    x1[k]= (k-1)*dx;
    pr[k]= (  2*x1[k] / ((stepsize^2)*N))* math.exp(-1*x1[k]*x1[k]/((stepsize^2)*N)  )

pr = pr * (bins*M);   #rescale theoretical distribution

#plot histogram and theoretical distribution
plt.hist(rad, nn)
plt.plot(x1, pr, 'r')
plt.ylabel('Number of walkers')

Problem 5, (chain.m translated to python by hand)

In [95]:
#  Thanks for making me do this conversion by hand...
#  chain.m simulates a Markov chain on {1, 2, ..., n} given an initial 
#  distribution and a transition matrix. 

# Program to simulate a Markov chain

# Below is a sample
k = 3 #number of states
v = np.mat([1/3, 1/3, 1/3]) # initial distribution
M = np.mat([[.5, 0.5, 0], [0.5, 0.5, 0], [0, 0.5, 0.5]]) # Markov transition matrix

n = 80 # number of time steps to take
x = v
t = [range(0, n)] # time indices


#Calculate v for n times steps.
for i in range (0,n):
    x = np.concatenate((x,(x[i]*M)), axis=0)

#print (x)
xt = np.transpose(x)
#plot the results
plt.plot(t, xt[0], 'ro')
plt.plot(t, xt[1], 'go')
plt.plot(t, xt[2], 'bo')
plt.xlabel('time step')
#plt.legend('Island A', 'Island B', 'Island C')

Problem 6; Chain Final

In [110]:
# Program to simulate a Markov chain

# Below is a sample
k = 3 #number of states
v = np.mat([0.25, 0.25, 0.25, 0.25]) # initial distribution
M = np.mat([[0.00, 0.33, 0.34, 0.33], 
            [0.34, 0.00, 0.33, 0.33], 
            [0.33, 0.33, 0.00, 0.34],
            [0.33, 0.34, 0.33, 0.00]]) # Markov transition matrix

n = 80 # number of time steps to take
x = v
t = [range(0, n)] # time indices


#Calculate v for n times steps.
for i in range (0,n):
    x = np.concatenate((x,(x[i]*M)), axis=0)

#print (x)
xt = np.transpose(x)
#plot the results
plt.plot(t, xt[0], 'r.', label="Island A")
plt.plot(t, xt[1], 'g.', label="Island B")
plt.plot(t, xt[2], 'b.', label="Island C")
plt.plot(t, xt[3], 'y.', label="Island D")
plt.xlabel('time step')
plt.ylabel('proportion of population')
plt.legend(['Island A', 'Island B', 'Island C', 'Island D'])
plt.title("Question 6 part B; No one stays")
In [42]:
#oh joy another matlab file to transcribe.  One of the two of us is really stubborn.  
#(It might be the one writing pithy comments no living human but himself will ever read)

def run(numx, periodic):
    numx = 101   #number of grid points in x
    numt = 2000  #number of time steps to be iterated
    dx = 1/(numx)
    dt = 0.00005
    pi = math.pi
    x = np.arange(0,1,dx)   #vector of x values, to be used for plotting
    periodic = False

    C = np.zeros([numx,numt])   #initialize everything to zero

    #specify initial conditions
    t = [0]      #t=0

    for i in range(0,numx-1):
        C[i] = math.cos(2*pi*x[i]-pi)+1

    #iterate difference equations
        for j in range(0, numt-1):
            t.append(t[j] + dt)  #update time
            for i in range(0, numx-1):  #loop through the middle spatial steps
                C[i][j+1] = C[i][j] + (dt/dx**2)*(C[i+1][j] - 2*C[i][j] + C[i-1][j]) 
            tmp = C[0][j+1]
            C[0][j+1] = C[numx-1][j+1]
            C[numx-1][j+1] = tmp
        for j in range(0, numt-1):
            t.append(t[j] + dt)  #update time
            for i in range(1, numx-2):  #loop through the middle spatial steps
                C[i][j+1] = C[i][j] + (dt/dx**2)*(C[i+1][j] - 2*C[i][j] + C[i-1][j]) 
            C[0][j+1] = 0       #Dirichlet boundary conditions
            C[numx-1][j+1] = 0  #Dirichlet boundary conditions

    D = np.transpose(C)

    plt.plot(x,D[1],    '-')
    plt.plot(x,D[9],    'b-')
    plt.plot(x,D[99],   'r--')
    plt.plot(x,D[999],  'y-.')
    plt.plot(x,D[1999], 'g:')
numx = 0
periodic = False

Hydrostatics Sheet (to Navier-Stokes)

Examining the incompressible form $$\left(\frac{\partial u(\vec{x}, t)}{\partial t} + u(x,t)\cdot\nabla u(\vec{x},t)\right)\rho=-\nabla P(\vec{x},t)+\mu \Delta u(\vec{x},t)+\rho \vec{g}$$ Using the helpful assumption that u = 0 at all locations in the fluid and at all times we can simplify: $$0 = -\nabla P(\vec{x},t)+\rho \vec{g}$$ $$\rho\vec{g}=\nabla P(\vec{x},t)$$ $$\begin{pmatrix}0\\0\\-\rho g\end{pmatrix} = \begin{pmatrix}\frac{\partial P}{\partial x}\\\frac{\partial P}{\partial y}\\\frac{\partial P}{\partial z}\end{pmatrix}$$ $$-\rho g = \frac{\partial P}{\partial z}$$ $$-\rho g\int_{\Omega}dz = P$$ $$P = -\rho gz + P_0$$ This is an encouraging result, as it looks very similar to the pressure equation $P=\rho gh$.

Problem 1, Flow Through a Channel (Homework 8)

In this problem you will consider two different fluid models and find the fluid flow profile within a rectangular channel, as in the figure below

In the first part, we consider the inviscid, incompressible Navier-Stokes equations. As the name suggests, the inviscid models assumes the fluid does not undergo any shearing-typeforces. This model is also known as the Euler Equations. The governing equations for this situation are given by the following equations

$$\left(\frac{\partial u(\vec{x}, t)}{\partial t} + u(\vec{x},t)\cdot\nabla u(\vec{x},t)\right)\rho=-\nabla P(\vec{x},t)$$$$\nabla\cdot\vec{u}=0$$

Find the velocity of the fluid through the channel.

Now in this part of the problem, we will consider the viscous (Newtonian), incompressible Navier-Stokes equations to model the situation in Figure 1. Again, as the name suggests here, we will reintroduce viscous forces back into the model. The governing equations for this situation are given by the following equations

$$\left(\frac{\partial u(\vec{x}, t)}{\partial t} + u(\vec{x},t)\cdot\nabla u(\vec{x},t)\right)\rho=-\nabla P(\vec{x},t)+\mu\Delta\vec{u}$$$$\nabla\cdot\vec{u}=0$$

a) $$u(\vec{x},t)=\begin{pmatrix}u_x(\vec{x},t)\\u_y(\vec{x},t)\\u_z(\vec{x},t)\end{pmatrix}$$ $$ \rho \left( \begin{pmatrix} \frac{\partial u_x(\vec{x},t)}{\partial t} \\ \frac{\partial u_y(\vec{x},t)}{\partial t} \\ \frac{\partial u_z(\vec{x},t)}{\partial t} \end{pmatrix} + \begin{pmatrix} \frac{\partial u_x(\vec{x},t)}{\partial x} & \frac{\partial u_x(\vec{x},t)}{\partial y} & \frac{\partial u_x(\vec{x},t)}{\partial z}\\ \frac{\partial u_y(\vec{x},t)}{\partial x} & \frac{\partial u_y(\vec{x},t)}{\partial y} & \frac{\partial u_y(\vec{x},t)}{\partial z}\\ \frac{\partial u_z(\vec{x},t)}{\partial x} & \frac{\partial u_z(\vec{x},t)}{\partial y} & \frac{\partial u_z(\vec{x},t)}{\partial z} \end{pmatrix} \begin{pmatrix}u_x(\vec{x},t)\\u_y(\vec{x},t)\\u_z(\vec{x},t)\end{pmatrix} \right)= -\begin{pmatrix} \frac{\partial P(\vec{x},t)}{\partial x} \\ \frac{\partial P(\vec{x},t)}{\partial y} \\ \frac{\partial P(\vec{x},t)}{\partial z}\end{pmatrix}$$ We are concerned here with the velocity in the x direction so we can ignore the $y$ and $z$ dimensions,

$$u_x(\vec{x},t) =-\left(\frac{\partial P(\vec{x},t)}{\partial x}\frac{1}{\rho} + \frac{\partial u_x(\vec{x},t)}{\partial t} + \frac{\partial u_x(\vec{x},t)}{\partial y}\cdot u_y(\vec{x},t) + \frac{\partial u_z(\vec{x},t)}{\partial z}\cdot u_z(\vec{x},t)\right) \left(\frac{\partial u_x(\vec{x},t)}{\partial x}\right)^{-1}$$

We will approximate the change in pressure over the length $L$ as $\frac{p_2-p_1}{L}$ so we have $$u_x(\vec{x},t) =-\left(\frac{p_2-p_1}{\rho L} + \frac{\partial u_x(\vec{x},t)}{\partial t} + \frac{\partial u_x(\vec{x},t)}{\partial y}\cdot u_y(\vec{x},t) + \frac{\partial u_z(\vec{x},t)}{\partial z}\cdot u_z(\vec{x},t)\right) \left(\frac{\partial u_x(\vec{x},t)}{\partial x}\right)^{-1}$$

Now we add the viscous forces to the left hand side $$ \rho \left( \begin{pmatrix} \frac{\partial u_x(\vec{x},t)}{\partial t} \\ \frac{\partial u_y(\vec{x},t)}{\partial t} \\ \frac{\partial u_z(\vec{x},t)}{\partial t} \end{pmatrix} + \begin{pmatrix} \frac{\partial u_x(\vec{x},t)}{\partial x} & \frac{\partial u_x(\vec{x},t)}{\partial y} & \frac{\partial u_x(\vec{x},t)}{\partial z}\\ \frac{\partial u_y(\vec{x},t)}{\partial x} & \frac{\partial u_y(\vec{x},t)}{\partial y} & \frac{\partial u_y(\vec{x},t)}{\partial z}\\ \frac{\partial u_z(\vec{x},t)}{\partial x} & \frac{\partial u_z(\vec{x},t)}{\partial y} & \frac{\partial u_z(\vec{x},t)}{\partial z} \end{pmatrix} \begin{pmatrix}u_x(\vec{x},t)\\u_y(\vec{x},t)\\u_z(\vec{x},t)\end{pmatrix} \right)= -\begin{pmatrix} \frac{\partial P(\vec{x},t)}{\partial x} \\ \frac{\partial P(\vec{x},t)}{\partial y} \\ \frac{\partial P(\vec{x},t)}{\partial z}\end{pmatrix} + \begin{pmatrix} \frac{\partial^2u_x(\vec{x},t)}{\partial x^2} + \frac{\partial^2u_x(\vec{x},t)}{\partial y^2} + \frac{\partial^2u_x(\vec{x},t)}{\partial z^2} \\ \frac{\partial^2u_y(\vec{x},t)}{\partial x^2} + \frac{\partial^2u_y(\vec{x},t)}{\partial y^2} + \frac{\partial^2u_y(\vec{x},t)}{\partial z^2} \\ \frac{\partial^2u_z(\vec{x},t)}{\partial x^2} + \frac{\partial^2u_z(\vec{x},t)}{\partial y^2} + \frac{\partial^2u_z(\vec{x},t)}{\partial z^2} \end{pmatrix} $$

$$u_x(\vec{x},t) =-\left( \left(\frac{\partial P(\vec{x},t)}{\partial x} + \frac{\partial^2u_x(\vec{x},t)}{\partial x^2} + \frac{\partial^2u_x(\vec{x},t)}{\partial y^2} + \frac{\partial^2u_x(\vec{x},t)}{\partial z^2} \right) \frac{1}{\rho} + \frac{\partial u_x(\vec{x},t)}{\partial t} + \frac{\partial u_x(\vec{x},t)}{\partial y}\cdot u_y(\vec{x},t) + \frac{\partial u_z(\vec{x},t)}{\partial z}\cdot u_z(\vec{x},t)\right) \left(\frac{\partial u_x(\vec{x},t)}{\partial x}\right)^{-1}$$

Problem 2, Reynold's Number Calculation (Homework 8)

In this problem you are asked to calculate the $Re$ for a few biological systems. Recall that the Reynolds Number, $Re$, is

$$Re=\frac{\rho LV}{\mu}$$

First we consider a human embryonic heart around 3 weeks post fertilization. At this stage of development, the heart is a valveless tube, that pumps blood using peristalsis. The diameter of the heart tube is roughly $0.5\cdot10^{-4}m$ and the peak velocity of the blood is approximately $1\frac{mm}{s}$. Although physical properties of embryonic blood are not explicitly known, one tends to approximate the density as $1025\frac{kg}{m^3}$ and kinematic viscosity of $2\cdot10^{-3}Pa\cdot s$. What is the $Re$ for this cardiac system?

Once the heart has fully matured to a multi-chambered valvular pumping system, it has a different $Re$. However, rather than compute the $Re$ of the adult heart, you will calculate the maximum velocity of the blood running through the aorta, e.g., the main artery of the body that is responsible for supplying oxygenated blood throughout the entire circulatory system. The adult human ascending aorta is roughly $2.3 cm$ in diameter on average. The properties of adult human blood are better documented with a density of approximately $1060 \frac{kg}{m^3}$ and kinematic viscosity of roughly $3\cdot10^{-3}Pa\cdot s$. If the $Re$ of the adult aorta is 10 000, what is the maximum velocity of the blood running through it?

In the above two cases, you have calculated the either the $Re$ or parameters involved for blood flowing through different vessels. Having learned what the $Re$ was from part (a) and the maximum velocity of blood flow in part (b), what can you say about the human heart as it develops?

a) $$Re=\frac{\rho LV}{\mu} = \frac{\left(1025\frac{kg}{m^3}\right)\left(0.5\cdot10^{-4}m\right)\left(1\cdot10^{-3}\frac{m}{s}\right)}{2\cdot10^{-3}\frac{kg}{s\cdot m}}=0.025625$$

b) $$Re=\frac{\rho LV}{\mu} = \frac{\left(1060\frac{kg}{m^3}\right)\left(2.3\cdot10^{-3}m\right)\left(V\right)}{3\cdot10^{-3}\frac{kg}{s\cdot m}}=10000$$ $$V=12.3052\frac{m}{s}$$

c) The vorticity of the embryonic heart is far less than that off a fully developed one and much of this comes from the massive increase in peak velocity. The embryonic heart seems to deal with a much simpler, less turbulent fluid system.

Problem 3, Fluid Nondimensionalization (Homework 8)

For this problem, find a non-dimensional number in fluid dynamics (other than the Reynolds Number), and write out what parameters it depends on. You do not have to derive it. Briefly discuss what the non-dimensional number describes, and find two fluids applications it is used in (hint: cough, cough - Google “dimensionless fluids numbers”).

One interesting case is the Stanton Number, which governs forced convection. That is, it is used to characterize heat flow from a boundary to a fluid. It is constructed from three other nondimensional numbers: the Nusselt Number, the Prandtl Number, and the Reynolds Number.

The Nusselt Number is the ratio of convective to conductive heat transfer from a boundary to the surface of a fluid. Essentially it seems to describe how temperature differences propagate from the surface of a fluid to that fluid's body. The number is defined as $$Nu=\frac{hL}{k}=$$ $$=\frac{\text{heat transfer coefficient}\left( \frac{kg}{s^3\cdot K} \right)\text{characteristic length}\left( m \right)}{\text{thermal conductivity}\left( \frac{kg\cdot m}{s^3\cdot K} \right)}$$

The Prandtl Number is a ratio of viscous diffusion to thermal diffusion. Thus, it seems to define how much heat diffusion in a material is due to standard conduction versus how much is due to convection currents within the material. It is simply $$Pr=\frac{\nu}{\alpha}=\frac{\text{momentum diffusivity}\left( \frac{m^2}{s} \right)}{\text{thermal diffusivity}\left( \frac{m^2}{s} \right)}$$

The last number, the Reynolds Number, has already been covered exhaustively enough that I will avoid examining it in the same way.

From these three, we can construct the Stanton Number as $$St=\frac{Nu}{Pr\cdot Re}$$ What does this tell us? We seem to be dealing with a ratio between the convection of heat into a fluid and the movement of heat within that fluid (as the Prandtl Number describes the movement of heat by conduction and by current and the Reynolds number describes those currents). Taken together, this number seems to wholly capture the process by which temperature differences propagate from a boundary to within the fluid until the system reaches equilibrium.

I can see many possibly applications for this number. One case would be the pot of warming (though not boiling) water; the hot bottom and warm sides along with the cool top would create convection currents while the heat also conducts itself through the water statically. I could even imagine a more complete model that took into account the diffusion of the pot (via the heat equation) and modeled the air as a separate, large system.

The Stanton Number might also help with situations in which materials must be heated or cooled very precisely (telescope refractor mirrors come to mind; as do silicone wafers for microchip production). How best should heating elements be placed to ensure even heating? How should air currents be directed to best control cooling? How does a change in humidity or air composition affect outcomes? These are questions the Stanton Number could help answer (although exact numerical modelling is probably required, especially for those first two).

Problem 4, Rubber Band Simulation (Homework 8)

Here you are asked to run a fluid-structure interaction simulation using the immersed boundary method, and

  1. Take a screenshot of the simulation

  2. Describe what is happening within the simulation

This is the simulation at t=60:

From what I can gather, this example models a 2 dimensional rubber band in a fluid which is pulled by springs at the Lagrangian points of the system. I am not entirely sure what Lagrangian points are in the context of fluids, but in gravitational systems they are points of stability between two bodies, so perhaps they are still points in the system? In either case, this fluid seems to be under periodic and exponentially decaying forcing which creates waves of pressure; at first these are simple with only four main waves ossilating back and forth from t=0 to around t=40. This doubles to eight main waves after t=40 and twenty-four distinct waves can be seen at t=60.