From f224e2299b9906ed0073da2500974667a3853fd0 Mon Sep 17 00:00:00 2001
From: Adam Jozefiak <adam.jozefiak@alumni.ubc.ca>
Date: Mon, 25 Mar 2019 00:19:00 -0700
Subject: [PATCH 1/6] working copy

---
 docs/examples/FD_intro.jmd | 292 +++++++++++++++++++++++++++++++++++++
 1 file changed, 292 insertions(+)
 create mode 100644 docs/examples/FD_intro.jmd

diff --git a/docs/examples/FD_intro.jmd b/docs/examples/FD_intro.jmd
new file mode 100644
index 0000000..e3f76ee
--- /dev/null
+++ b/docs/examples/FD_intro.jmd
@@ -0,0 +1,292 @@
+
+```julia
+# Packages
+using Plots, LinearAlgebra, SimpleDifferentialOperators
+```
+
+# Discretizing Time-Independent Linear Operators without Control
+
+This is a tutorial that demonstrates how to discretize a linear operator of a time-independent diffusion process without any control. As an example a univariate diffusion process is utilizedin and the corresponding  Hamilton-Jacobi-Bellman equation is solved.
+
+Notation:
+
+- Define $\mathcal{A}$ to be an operator (or infinitesimal generator) associated with a particular stochastic process.
+
+- Define $q^{-} \equiv \min\{q,0\}$
+
+- Define $q^{+} \equiv \max\{q,0\}$
+
+- Derivatives are denoted by the operator $\partial$ and univariate derivatives are defined as $\partial_x v(x) = v'(x)$
+
+###### Grids:
+
+- For the univariate case, $x\in [\underline{x}, \overline{x}]$ we will discretize $x$ with $I$ uniformly spaced points by $\{x_i\}_{i=1}^I$ such that $x_1 = \underline{x}$ and $x_I=\overline{x}$. The distance between consecutive points is denoted by $\Delta = x_{i+1}-x_i$.
+
+### Time Independent Univariate Diffusion
+
+Assume that the time variable does not directly enter into anything important in the problem (e.g. drifts, payoff, stopping values, boundary values, or discount rates). This is also assuming that the drift, payoffs, stopping value are not a (direct) function of any $v_i$. Otherwise, the operator would be nonlinear. In the case of a linear $\mathcal{A}$, the resulting discretized operator will be a (heavily sparse) matrix $A\in\mathbb{R}^{I\times I}$, which is especially conveneint to solve.
+
+#### Stochastic Process and Boundary Values
+
+Take a diffusion process for $x_t$ according to the following stochastic difference equation,
+
+$$dx_t = \mu(x_t)dt + \sigma(x_t)d\mathbb{W}_t  \hspace{1cm}(1)$$
+
+where $\mathbb{W}_t$ is Brownian motion and $x\in [\underline{x}, \overline{x}]$ where $-\infty \leq \underline{x}<\overline{x}\leq \infty$. The functions $\mu(\cdot)$ and $\sigma(\cdot)$ are left general, but it is essential that they are time-homogenous.
+
+For instance, we can have the following functions $\mu(\cdot)$ and $\sigma(\cdot)$:
+
+- Define $\mu(x) = -0.01x$
+
+- Define $\sigma(x) = 0.1x$
+
+```julia
+function μ(x)
+    return -0.01*x
+end
+
+function σ2(x)
+    return (0.1*x)^2
+end
+```
+
+The partial differential operator (infinitesimal generator) associated with the stochastic process in (1) is
+
+$$\mathcal{A} \equiv \mu(x)\partial_x + \frac{\sigma(x)^2}{2}\partial_{xx}\hspace{1cm}  (2)$$
+
+Since the process is univariate, we require boundary conditions at both sides of $x$, and assume that the boundary and $\mu(\cdot)$ and $\sigma(\cdot)^2$ are independent of time. For this tutorial we will consider the boundary values of a reflecting barrier, for sides of $x$, for an operator on a generic function $v(\cdot)$.
+
+For reflecting barriers at both sides of $x$ we have the following boundary values:
+
+- Let $\partial_{x}v(\underline{x}) = v'(\underline{x}) = 0$
+
+- Let $\partial_{x}v(\overline{x}) = v'(\overline{x}) = 0$
+
+### Example Problem
+
+Using $\mu(\cdot)$ and $\sigma(\cdot)$ as defined earlier $\mu(x) = -0.01x$ and $\sigma(x) = 0.1x$ to define the stochastic difference equation
+
+$$dx_t = \mu(x_t)dt + \sigma(x_t)d\mathbb{W}_t$$
+
+we can consider the problem where:
+
+- Utility function: $u(x) = \log(x)$
+
+- Left Boundary: $\underline{x} = 0.01$
+
+- Right Boundary: $\overline{x} = 10$
+
+- Discount factor: $\rho = 0.05$
+
+Note that the choices of $\mu(\cdot)$, $\sigma(\cdot)$, and $u(\cdot)$ render this particular problem control-free. That is, the state variable $x_t$ determines the drift and variance of the stochastics process along with determining the utility without any control by the agent.
+
+We can then discretize $[\underline{x},\overline{x}]=[0.01,10]$ by increments of $0.001$.
+
+```julia
+x_min = 0.01
+x_max = 10
+
+grid = x_min:0.001:x_max
+
+function u(x)
+    return log(x)
+end
+
+ρ = 0.05
+```
+
+### Solving Control-Free Time-Independent Linear Hamilton-Jacobi Bellman Equation
+
+In the simple case of an uncontrolled linear stochastic process with no optimal stopping, if an agent has payoff $u(x)$ with discounting $\rho>0$, then the following Hamilton-Jacobi Bellman equation determines the value
+
+$$\rho v(x) = u(x) + \mathcal{A}v(x)$$
+
+subject to the appropriate boundary conditions. Discretizing the grid, let $u \equiv \{u(x_i)\}_{i=1}^I\in\mathbb{R}^I$, where $v\in\mathbb{R}^I$ is unknown, then the discretized version is,
+
+$$\rho v = u + Av$$
+
+where $A$ is a discretization of $\mathcal{A}$ on $\{x_i\}_{i=1}^I$. Note that $A$ necessarily must account for the boundary values of the grid.
+
+Next, we can solve algebraically for $v$ as the linear system:
+
+$$(\rho I -A)v = u$$
+
+We will construct the discretized operator $A$ using the SimpleDifferentialOperators package. Since the drift of the stochastic process (1) is always negative in our example on $[0.01,10.0]$ we will only utilize the downwind drift discretization $L_{1-}$. We note that the resulting matrix $A$ is a tridiagonal matrix, making it particularly convenient to solve numerically.
+
+```julia
+bc = (Reflecting(), Reflecting())
+
+A = μ.(grid).*L₁₋(grid, bc) + (σ2.(grid)/2).*L₂(grid, bc)
+
+v = (ρ*I - A)\ u.(grid)
+
+plt_u = plot(grid, u.(grid), label="u", xlabel="x", ylabel = "u(x)", lw=2, legend=:topleft)
+plt_v = plot(grid, v, label="v", xlabel="x", ylabel="v(x)", lw=2, legend=:topleft)
+plot(plt_u,plt_v)
+```
+
+### Appendix: Upwind Finite Difference Discretization of the Differential Operator with a Unifrom Grid
+
+Here, we will manually use finite difference to discretize $\mathcal{A}$ in order to construct $A$ as a means to solving the HJBE. We note that utilizing the SimpleDifferentialOpperators package avoids this (tedious) derivation.
+
+Note that
+
+$$\mathcal{A}v=\mu(x) v_x + \frac{\sigma^2(x)}{2}v_{xx}$$
+
+Next, we will proceed to derive an upwind finite difference scheme for $\mathcal{A}v$. In particular we will employ a uniformly spaced discretization of $[\underline{x},\overline{x}]$. Denote this discretization of $[\underline{x},\overline{x}]$ with $I$ uniformly spaced points as $x \equiv \{x_i\}_{i=1}^I$, where $x_1 = \underline{x}$ and $x_I=\overline{x}$. We will denote the distance between consecutive points by $\Delta = x_{i+1}-x_i$.
+
+Then for an interior $x_i$, that is $i\in[2,I-1]$, $\mathcal{A}v(x_i)$ can be discretized as such:
+
+$$\mathcal{A}v(x_i)\approx \mu_i^+(\frac{v_{i+1}-v_{i}}{\Delta}) + \mu_i^-(\frac{v_i-v_{i-1}}{\Delta}) +\frac{\sigma^2_i}{2}\frac{v_{i+1}-2v_i+v_{i-1}}{\Delta^2}$$
+
+such that the upwind scheme chooses forward or backwrd difference depending on the sign of the drift.
+
+We can collect the above terms such that:
+
+$$\mathcal{A}v(x_i)\approx (-\frac{\mu_i^-}{\Delta}+\frac{\sigma^2_i}{2\Delta^2})v_{i-1} + (\frac{\mu_i^-}{\Delta} - \frac{\mu_i^+}{\Delta}-\frac{\sigma^2_i}{\Delta^2})v_i + (\frac{\mu_i^+}{\Delta}+\frac{\sigma^2_i}{2\Delta^2})v_{i+1}$$
+
+For the boundary value at $\underline{x}$, consider the ghost node $x_0$. Then the approximation of $\mathcal{A}v(x_1)$ is:
+
+$$\mathcal{A}v(x_1)\approx (-\frac{\mu_1^-}{\Delta}+\frac{\sigma^2_1}{2\Delta^2})v_{0} + (\frac{\mu_1^-}{\Delta} - \frac{\mu_1^+}{\Delta}-\frac{\sigma^2_1}{\Delta^2})v_1 + (\frac{\mu_1^+}{\Delta}+\frac{\sigma^2_1}{2\Delta^2})v_{2}$$
+
+Since $v'(\underline{x}) = 0$ then $v'(\underline{x})$ can be approximated since $v'(\underline{x}) = \frac{v_1-v_0}{\Delta} = 0$, then $v_0 = v_1$. Then:
+
+$$\mathcal{A}v(x_1)\approx [(-\frac{\mu_1^-}{\Delta}+\frac{\sigma^2_1}{2\Delta^2}) + (\frac{\mu_1^-}{\Delta} - \frac{\mu_i^+}{\Delta}-\frac{\sigma^2_1}{\Delta^2})]v_1 + (\frac{\mu_1^+}{\Delta}+\frac{\sigma^2_1}{2\Delta^2})v_{2}$$
+
+Likewise, for the boundary value at $\overline{x}$, consider the ghost node $x_{I+1}$. Then the approximation of $\mathcal{A}v(x_I)$ is:
+
+$$\mathcal{A}v(x_I)\approx (-\frac{\mu_I^-}{\Delta}+\frac{\sigma^2_I}{2\Delta^2})v_{I-1} + (\frac{\mu_I^-}{\Delta} - \frac{\mu_I^+}{\Delta}-\frac{\sigma^2_I}{\Delta^2})v_I + (\frac{\mu_I^+}{\Delta}+\frac{\sigma^2_I}{2\Delta^2})v_{I+1}$$
+
+Since $v'(\overline{x}) = 0$ then $v'(\overline{x})$ can be approximated since $v'(\overline{x}) = \frac{v_{I+1}-v_I}{\Delta} = 0$, then $v_{I+1} = v_I$. Then:
+
+$$\mathcal{A}v(x_I)\approx (-\frac{\mu_I^-}{\Delta}+\frac{\sigma^2_I}{2\Delta^2})v_{I-1} + [(\frac{\mu_I^-}{\Delta} - \frac{\mu_I^+}{\Delta}-\frac{\sigma^2_I}{\Delta^2}) + (\frac{\mu_I^+}{\Delta}+\frac{\sigma^2_I}{2\Delta^2})]v_{I}$$
+
+Thus the formulas for the discretization $A$ of $\mathcal{A}$ on the uniform grid $\{x_i\}_{i=1}^I$ are:
+
+- For $i\in[2,I-1]$, $\mathcal{A}v(x_i)\approx (-\frac{\mu_i^-}{\Delta}+\frac{\sigma^2_i}{2\Delta^2})v_{i-1} + (\frac{\mu_i^-}{\Delta} - \frac{\mu_i^+}{\Delta}-\frac{\sigma^2_i}{\Delta^2})v_i + (\frac{\mu_i^+}{\Delta}+\frac{\sigma^2_i}{2\Delta^2})v_{i+1}$
+
+- For $i=1$, $\mathcal{A}v(x_1)\approx [(-\frac{\mu_1^-}{\Delta}+\frac{\sigma^2_1}{2\Delta^2}) + (\frac{\mu_1^-}{\Delta} - \frac{\mu_1^+}{\Delta}-\frac{\sigma^2_1}{\Delta^2})]v_1 + (\frac{\mu_1^+}{\Delta}+\frac{\sigma^2_1}{2\Delta^2})v_{2}$
+
+- For $i=I$, $\mathcal{A}v(x_I)\approx (-\frac{\mu_I^-}{\Delta}+\frac{\sigma^2_I}{2\Delta^2})v_{I-1} + [(\frac{\mu_I^-}{\Delta} - \frac{\mu_I^+}{\Delta}-\frac{\sigma^2_I}{\Delta^2}) + (\frac{\mu_I^+}{\Delta}+\frac{\sigma^2_I}{2\Delta^2})]v_{I}$
+
+To succinctly construct $A$ in our code,  the discretization of $\mathcal{A}$ on the discretized grid $\{x_i\}_{i=1}^I$, we can define $X,Y,$ and $Z$ such that:
+
+- Let $X = -\frac{\mu^-}{\Delta}+\frac{\sigma^2}{2\Delta^2}$
+
+- Let $Y = -\frac{\mu^+}{\Delta}+\frac{\mu^-}{\Delta^2}-\frac{\sigma^2}{\Delta^2}$
+
+- Let $Z = \frac{\mu^+}{\Delta}+\frac{\sigma^2}{2\Delta^2}$
+
+Where $\mu^- \equiv \{\mu(x_i)^-\}_{i=1}^I$, $\mu^+ \equiv \{\mu(x_i)^+\}_{i=1}^I$, and $\sigma^2 \equiv \{\sigma(x_i)^2\}_{i=1}^I$
+
+Then $A$ can be constucted as a tridiagonal matrix:
+
+$$A=\begin{bmatrix}
+Y_1+X_1 & Z_1 & 0 & \cdots & \cdots & \cdots & 0 \\
+X_2 & Y_2 & Z_2 & 0 & \ddots &  & \vdots \\
+0 & \ddots & \ddots & \ddots & \ddots & & \vdots \\
+\vdots & & \ddots & \ddots & \ddots & \ddots & \vdots \\
+\vdots & & & \ddots & X_{I-1} & Y_{I-1} & Z_{I-1} \\
+0 & \cdots & \cdots & \cdots & 0 & Y_{I} & Y_I+Z_{I} \\
+\end{bmatrix}$$
+
+It thus follows that $\mathcal{A}v(x) \approx Av(x)$ on the discretized grid $x\equiv\{x_{i=1}\}_{i=1}^I$
+
+
+### Example Problem
+
+Using $\mu(\cdot)$ and $\sigma(\cdot)$ as defined earlier $\mu(x) = -0.01x$ and $\sigma(x) = 0.1x$ to define the stochastic difference equation
+
+$$dx_t = \mu(x_t)dt + \sigma(x_t)d\mathbb{W}_t$$
+
+we can consider the problem where:
+
+- Utility function: $u(x) = \log(x)$
+
+- Left Boundary: $\underline{x} = 0.01$
+
+- Right Boundary: $\overline{x} = 10$
+
+- Discount factor: $\rho = 0.05$
+
+Note that the choices of $\mu(\cdot)$, $\sigma(\cdot)$, and $u(\cdot)$ render this particular problem control-free. That is, the state variable $x_t$ determines the drift and variance of the stochastics process along with determining the utility without any control by the agent.
+
+We can then discretize $[\underline{x},\overline{x}]=[0.01,10]$ by increments of $0.001$.
+
+```julia
+x_min = 0.01
+x_max = 10
+
+grid = x_min:0.001:x_max
+
+function μ(x)
+    return -0.01*x
+end
+
+function σ2(x)
+    return (0.1*x)^2
+end
+
+function u(x)
+    return log(x)
+end
+
+n = size(grid)[1]
+
+delta = grid[2]-grid[1] # Create the delta with respect to x grid
+delta_2 = delta^2
+
+function μ_p(x)
+    return max(μ(x),0)
+end
+
+function μ_m(x)
+    return min(μ(x),0)
+end
+
+X = (-μ_m.(grid)/delta) + (σ2.(grid)/(2*delta_2))
+Y = (-μ_p.(grid)/delta) + (μ_m.(grid)/delta) - (σ2.(grid)/delta_2)
+Z = (μ_p.(grid)/delta) + (σ2.(grid)/(2*delta_2))
+
+A_2 = Tridiagonal(X[2:n], Y, Z[1:n-1])
+A_2[1,1] = Y[1] + X[1]
+A_2[n,n] = Y[n] + Z[n]
+```
+
+### Solving Control-Free Time-Independent Linear Hamilton-Jacobi Bellman Equation
+
+In the simple case of an uncontrolled linear stochastic process with no optimal stopping, if an agent has payoff $u(x)$ with discounting $\rho>0$, then the following Hamilton-Jacobi Bellman equation determines the value
+
+$$\rho v(x) = u(x) + \mathcal{A}v(x)$$
+
+subject to the appropriate boundary conditions. Discretizing the grid, let $u \equiv \{u(x_i)\}_{i=1}^I\in\mathbb{R}^I$, where $v\in\mathbb{R}^I$ is unknown, then the discretized version is,
+
+$$\rho v = u + Av$$
+
+Which we can solve algebraically for $v$ as the linear system:
+
+$$(\rho I -A)v = u$$
+
+Note that $A$ is a tridiagonal matrix, making it particularly convenient to solve numerically.
+
+```julia
+ρ = 0.05
+
+v_2 = (ρ*I - A_2)\ u.(grid)
+
+plt_u = plot(grid, u.(grid), label="u", xlabel="x", ylabel = "u(x)", lw=2, legend=:topleft)
+plt_v2 = plot(grid, v_2, label="v_2", xlabel="x", ylabel="v(x)", lw=2, legend=:topleft)
+
+plot(plt_u,plt_v2)
+```
+
+```julia
+using Test
+
+# We can verify that the solution v obtained with the linear operator from the SimpleDifferentialOperators package
+# agrees with the solution v_2 obtained fromm manually discretizing the linear operator
+
+@test v[1] ≈ v_2[1]
+@test v[500] ≈ v_2[500]
+@test v[999] ≈ v_2[999]
+```

From 8986053224a81af351e7e7c0c780a1d100a89e2d Mon Sep 17 00:00:00 2001
From: Adam Jozefiak <adam.jozefiak@alumni.ubc.ca>
Date: Mon, 1 Apr 2019 09:42:01 -0700
Subject: [PATCH 2/6] FD introduction first pass

---
 docs/examples/FD_intro.jmd | 512 +++++++++++++++++++------------------
 1 file changed, 259 insertions(+), 253 deletions(-)

diff --git a/docs/examples/FD_intro.jmd b/docs/examples/FD_intro.jmd
index e3f76ee..2f1b19f 100644
--- a/docs/examples/FD_intro.jmd
+++ b/docs/examples/FD_intro.jmd
@@ -1,292 +1,298 @@
-
 ```julia
-# Packages
-using Plots, LinearAlgebra, SimpleDifferentialOperators
+using SimpleDifferentialOperators,Plots, LinearAlgebra
 ```
 
-# Discretizing Time-Independent Linear Operators without Control
-
-This is a tutorial that demonstrates how to discretize a linear operator of a time-independent diffusion process without any control. As an example a univariate diffusion process is utilizedin and the corresponding  Hamilton-Jacobi-Bellman equation is solved.
-
-Notation:
-
-- Define $\mathcal{A}$ to be an operator (or infinitesimal generator) associated with a particular stochastic process.
-
-- Define $q^{-} \equiv \min\{q,0\}$
-
-- Define $q^{+} \equiv \max\{q,0\}$
-
-- Derivatives are denoted by the operator $\partial$ and univariate derivatives are defined as $\partial_x v(x) = v'(x)$
-
-###### Grids:
-
-- For the univariate case, $x\in [\underline{x}, \overline{x}]$ we will discretize $x$ with $I$ uniformly spaced points by $\{x_i\}_{i=1}^I$ such that $x_1 = \underline{x}$ and $x_I=\overline{x}$. The distance between consecutive points is denoted by $\Delta = x_{i+1}-x_i$.
-
-### Time Independent Univariate Diffusion
-
-Assume that the time variable does not directly enter into anything important in the problem (e.g. drifts, payoff, stopping values, boundary values, or discount rates). This is also assuming that the drift, payoffs, stopping value are not a (direct) function of any $v_i$. Otherwise, the operator would be nonlinear. In the case of a linear $\mathcal{A}$, the resulting discretized operator will be a (heavily sparse) matrix $A\in\mathbb{R}^{I\times I}$, which is especially conveneint to solve.
-
-#### Stochastic Process and Boundary Values
-
-Take a diffusion process for $x_t$ according to the following stochastic difference equation,
-
-$$dx_t = \mu(x_t)dt + \sigma(x_t)d\mathbb{W}_t  \hspace{1cm}(1)$$
-
-where $\mathbb{W}_t$ is Brownian motion and $x\in [\underline{x}, \overline{x}]$ where $-\infty \leq \underline{x}<\overline{x}\leq \infty$. The functions $\mu(\cdot)$ and $\sigma(\cdot)$ are left general, but it is essential that they are time-homogenous.
-
-For instance, we can have the following functions $\mu(\cdot)$ and $\sigma(\cdot)$:
-
-- Define $\mu(x) = -0.01x$
-
-- Define $\sigma(x) = 0.1x$
+## Gentle Introduction to Upwind Finite Differences
+
+This notebook is a brief tutorial on upwind finite difference method for solving
+linear ordinary differential equations.
+
+In particular, this tutorial will be motivated by the following problem. Consider payoffs
+$u(x)$ for a stochastic process $x\in [0.2,04]$ with infitesimal generate $\mathcal{L}_x$
+discounted at rate $\rho > 0$. Then the payoffs can be valued by the
+Bellman equation:
+
+$$
+\begin{equation}
+  \rho v(x) = u(x) + \mathcal{L}_x v(x)
+\end{equation}
+$$
+
+We will consider the example $\mathcal{L}_x \equiv \mu \partial_x + \frac{\sigma^2}{2}\partial_{xx}$
+where $\mu = -0.1$ and $\sigma = 0.1$. We will consider $u(x) = \log(x)$ and we will impose the boundary
+conditions that $v'(0.2) = v'(0.4) = 0$.
+
+We will begin by first discretizing the interval $[0.2,0.4]$. Define $x = \{0.2, 0.3, 0.4\}$. Additionally,
+we will require the extended grid $\bar{x} = \{0.1,0.2,0.3,0.4,0.5\}$ in order to compute $\partial_x v(x)$ and
+$\partial_{xx}v(x)$ at the original boundaries $0.2$ and $0.4$.
+
+First, we will discretize $v'(x_i)$ for each $x_i\in x$.
+
+Note that by forward difference, for an arbitrary $x$ and $\Delta$:
+
+$$
+\begin{equation}
+  v'(x) \approx \frac{v(x+\Delta)-v(x)}{\Delta}
+\end{equation}
+$$
+
+Additionally, by backward difference,
+
+$$
+\begin{equation}
+  v'(x) \approx \frac{v(x)-v(x-\Delta)}{\Delta}
+\end{equation}
+$$
+
+Since in the stochastic process for $x$ the drift $\mu = -0.1 < 0$, we will
+approximate $v'(x)$ by the backward difference. In our instance, the grid $x$ has
+uniformly spaced grid points, so $\Delta = 0.1$ for each $x_i\in x$.
+
+So it follows that:
+
+$$
+\begin{equation}
+  v'(0.2) = 0 \text{ by the reflecting boundary condition}
+\end{equation}
+$$
+
+$$
+\begin{equation}
+  v'(0.3) \approx \frac{v(0.3)-v(0.2)}{0.1}
+\end{equation}
+$$
+
+$$
+\begin{equation}
+  v'(0.4) \approx \frac{v(0.4)-v(0.3)}{0.1}
+\end{equation}
+$$
+
+The above equations give us a linear system and hence we can represent them as follows:
+
+$$
+\begin{equation}
+  \begin{bmatrix}
+    v'(0.2) \\
+    v'(0.3) \\
+    v'(0.4) \\
+  \end{bmatrix}
+  \approx
+  \frac{1}{0.1}
+  \begin{bmatrix}
+    0 & 0 & 0 \\
+    -1 & 1 & 0 \\
+    0 & -1 & 1 \\
+  \end{bmatrix}
+  \begin{bmatrix}
+    v(0.2) \\
+    v(0.3) \\
+    v(0.4) \\
+  \end{bmatrix}
+\end{equation}
+$$
+
+Additionally, we will define $L_{1-}$ to be:
+
+$$
+\begin{equation}
+  L_{1-} = \frac{1}{0.1}
+  \begin{bmatrix}
+    0 & 0 & 0 \\
+    -1 & 1 & 0 \\
+    0 & -1 & 1 \\
+  \end{bmatrix}
+\end{equation}
+$$
 
 ```julia
-function μ(x)
-    return -0.01*x
-end
-
-function σ2(x)
-    return (0.1*x)^2
-end
+L_1m = 10*[0 0 0; -1 1 0; 0 -1 1]
 ```
 
-The partial differential operator (infinitesimal generator) associated with the stochastic process in (1) is
-
-$$\mathcal{A} \equiv \mu(x)\partial_x + \frac{\sigma(x)^2}{2}\partial_{xx}\hspace{1cm}  (2)$$
-
-Since the process is univariate, we require boundary conditions at both sides of $x$, and assume that the boundary and $\mu(\cdot)$ and $\sigma(\cdot)^2$ are independent of time. For this tutorial we will consider the boundary values of a reflecting barrier, for sides of $x$, for an operator on a generic function $v(\cdot)$.
-
-For reflecting barriers at both sides of $x$ we have the following boundary values:
-
-- Let $\partial_{x}v(\underline{x}) = v'(\underline{x}) = 0$
-
-- Let $\partial_{x}v(\overline{x}) = v'(\overline{x}) = 0$
-
-### Example Problem
-
-Using $\mu(\cdot)$ and $\sigma(\cdot)$ as defined earlier $\mu(x) = -0.01x$ and $\sigma(x) = 0.1x$ to define the stochastic difference equation
-
-$$dx_t = \mu(x_t)dt + \sigma(x_t)d\mathbb{W}_t$$
-
-we can consider the problem where:
-
-- Utility function: $u(x) = \log(x)$
-
-- Left Boundary: $\underline{x} = 0.01$
-
-- Right Boundary: $\overline{x} = 10$
-
-- Discount factor: $\rho = 0.05$
-
-Note that the choices of $\mu(\cdot)$, $\sigma(\cdot)$, and $u(\cdot)$ render this particular problem control-free. That is, the state variable $x_t$ determines the drift and variance of the stochastics process along with determining the utility without any control by the agent.
-
-We can then discretize $[\underline{x},\overline{x}]=[0.01,10]$ by increments of $0.001$.
-
-```julia
-x_min = 0.01
-x_max = 10
-
-grid = x_min:0.001:x_max
-
-function u(x)
-    return log(x)
-end
-
-ρ = 0.05
-```
-
-### Solving Control-Free Time-Independent Linear Hamilton-Jacobi Bellman Equation
-
-In the simple case of an uncontrolled linear stochastic process with no optimal stopping, if an agent has payoff $u(x)$ with discounting $\rho>0$, then the following Hamilton-Jacobi Bellman equation determines the value
-
-$$\rho v(x) = u(x) + \mathcal{A}v(x)$$
-
-subject to the appropriate boundary conditions. Discretizing the grid, let $u \equiv \{u(x_i)\}_{i=1}^I\in\mathbb{R}^I$, where $v\in\mathbb{R}^I$ is unknown, then the discretized version is,
-
-$$\rho v = u + Av$$
-
-where $A$ is a discretization of $\mathcal{A}$ on $\{x_i\}_{i=1}^I$. Note that $A$ necessarily must account for the boundary values of the grid.
-
-Next, we can solve algebraically for $v$ as the linear system:
-
-$$(\rho I -A)v = u$$
-
-We will construct the discretized operator $A$ using the SimpleDifferentialOperators package. Since the drift of the stochastic process (1) is always negative in our example on $[0.01,10.0]$ we will only utilize the downwind drift discretization $L_{1-}$. We note that the resulting matrix $A$ is a tridiagonal matrix, making it particularly convenient to solve numerically.
+(Note, a similar derivation using the forward differences would yield an analagous
+matrix $L_{1+}$)
+
+In a similar manner, we will discretize $v''(x_i)$ for each $x_i \in x = \{0.2,0.3, 0.4\}$.
+
+Note, by central difference, for arbitrary $x$ and $\Delta$:
+
+$$
+\begin{equation}
+  v''(x) \approx \frac{v(x+\Delta)-2v(x) + v(x-\Delta)}{\Delta^2}
+\end{equation}
+$$
+
+So it follows that:
+
+$$
+\begin{equation}
+  v''(0.2) \approx \frac{v(0.3)-2v(0.2) + v(0.1)}{0.1^2}
+\end{equation}
+$$
+
+Note that by the boundary conditions and backward difference:
+
+$$
+\begin{equation}
+0 = v'(0.2) \approx \frac{v(0.2)-v(0.1)}{0.1}
+\end{equation}
+$$
+
+which implies that $v(0.2) = v(0.1)$.
+
+Therefore,
+
+$$
+\begin{equation}
+  v''(0.2) \approx \frac{v(0.3)-v(0.2)}{0.1^2}
+\end{equation}
+$$
+
+Then $0.3$:
+
+$$
+\begin{equation}
+  v''(0.3) \approx \frac{v(0.4)-2v(0.3) + v(0.2)}{0.1^2}
+\end{equation}
+$$
+
+For $0.4$:
+
+$$
+\begin{equation}
+  v''(0.4) \approx \frac{v(0.5)-2v(0.4) + v(0.3)}{0.1^2}
+\end{equation}
+$$
+
+Note that by the boundary conditions and forward difference:
+
+$$
+\begin{equation}
+0 = v'(0.4) \approx \frac{v(0.5)-v(0.4)}{0.1}
+\end{equation}
+$$
+
+which implies that $v(0.4) = v(0.5)$.
+
+Therefore,
+
+$$
+\begin{equation}
+  v''(0.4) \approx \frac{-v(0.4) + v(0.3)}{0.1^2}
+\end{equation}
+$$
+
+The above equations give us a linear system and hence we can represent them as follows:
+
+$$
+\begin{equation}
+  \begin{bmatrix}
+    v''(0.2) \\
+    v''(0.3) \\
+    v''(0.4) \\
+  \end{bmatrix}
+  \approx
+  \frac{1}{0.1^2}
+  \begin{bmatrix}
+    -1 & 1 & 0 \\
+    1 & -2 & 1 \\
+    0 & 1 & -1 \\
+  \end{bmatrix}
+  \begin{bmatrix}
+    v(0.2) \\
+    v(0.3) \\
+    v(0.4) \\
+  \end{bmatrix}
+\end{equation}
+$$
+
+Additionally, we will define $L_2$ to be:
+
+$$
+\begin{equation}
+  L_2 = \frac{1}{0.1^2}
+  \begin{bmatrix}
+    -1 & 1 & 0 \\
+    1 & -2 & 1 \\
+    0 & 1 & -1 \\
+  \end{bmatrix}
+\end{equation}
+$$
 
 ```julia
-bc = (Reflecting(), Reflecting())
-
-A = μ.(grid).*L₁₋(grid, bc) + (σ2.(grid)/2).*L₂(grid, bc)
-
-v = (ρ*I - A)\ u.(grid)
-
-plt_u = plot(grid, u.(grid), label="u", xlabel="x", ylabel = "u(x)", lw=2, legend=:topleft)
-plt_v = plot(grid, v, label="v", xlabel="x", ylabel="v(x)", lw=2, legend=:topleft)
-plot(plt_u,plt_v)
+L_2 = 100*[-1 1 0; 1 -2 1; 0 1 -1]
 ```
 
-### Appendix: Upwind Finite Difference Discretization of the Differential Operator with a Unifrom Grid
-
-Here, we will manually use finite difference to discretize $\mathcal{A}$ in order to construct $A$ as a means to solving the HJBE. We note that utilizing the SimpleDifferentialOpperators package avoids this (tedious) derivation.
-
-Note that
-
-$$\mathcal{A}v=\mu(x) v_x + \frac{\sigma^2(x)}{2}v_{xx}$$
-
-Next, we will proceed to derive an upwind finite difference scheme for $\mathcal{A}v$. In particular we will employ a uniformly spaced discretization of $[\underline{x},\overline{x}]$. Denote this discretization of $[\underline{x},\overline{x}]$ with $I$ uniformly spaced points as $x \equiv \{x_i\}_{i=1}^I$, where $x_1 = \underline{x}$ and $x_I=\overline{x}$. We will denote the distance between consecutive points by $\Delta = x_{i+1}-x_i$.
-
-Then for an interior $x_i$, that is $i\in[2,I-1]$, $\mathcal{A}v(x_i)$ can be discretized as such:
-
-$$\mathcal{A}v(x_i)\approx \mu_i^+(\frac{v_{i+1}-v_{i}}{\Delta}) + \mu_i^-(\frac{v_i-v_{i-1}}{\Delta}) +\frac{\sigma^2_i}{2}\frac{v_{i+1}-2v_i+v_{i-1}}{\Delta^2}$$
-
-such that the upwind scheme chooses forward or backwrd difference depending on the sign of the drift.
-
-We can collect the above terms such that:
-
-$$\mathcal{A}v(x_i)\approx (-\frac{\mu_i^-}{\Delta}+\frac{\sigma^2_i}{2\Delta^2})v_{i-1} + (\frac{\mu_i^-}{\Delta} - \frac{\mu_i^+}{\Delta}-\frac{\sigma^2_i}{\Delta^2})v_i + (\frac{\mu_i^+}{\Delta}+\frac{\sigma^2_i}{2\Delta^2})v_{i+1}$$
-
-For the boundary value at $\underline{x}$, consider the ghost node $x_0$. Then the approximation of $\mathcal{A}v(x_1)$ is:
-
-$$\mathcal{A}v(x_1)\approx (-\frac{\mu_1^-}{\Delta}+\frac{\sigma^2_1}{2\Delta^2})v_{0} + (\frac{\mu_1^-}{\Delta} - \frac{\mu_1^+}{\Delta}-\frac{\sigma^2_1}{\Delta^2})v_1 + (\frac{\mu_1^+}{\Delta}+\frac{\sigma^2_1}{2\Delta^2})v_{2}$$
-
-Since $v'(\underline{x}) = 0$ then $v'(\underline{x})$ can be approximated since $v'(\underline{x}) = \frac{v_1-v_0}{\Delta} = 0$, then $v_0 = v_1$. Then:
-
-$$\mathcal{A}v(x_1)\approx [(-\frac{\mu_1^-}{\Delta}+\frac{\sigma^2_1}{2\Delta^2}) + (\frac{\mu_1^-}{\Delta} - \frac{\mu_i^+}{\Delta}-\frac{\sigma^2_1}{\Delta^2})]v_1 + (\frac{\mu_1^+}{\Delta}+\frac{\sigma^2_1}{2\Delta^2})v_{2}$$
-
-Likewise, for the boundary value at $\overline{x}$, consider the ghost node $x_{I+1}$. Then the approximation of $\mathcal{A}v(x_I)$ is:
-
-$$\mathcal{A}v(x_I)\approx (-\frac{\mu_I^-}{\Delta}+\frac{\sigma^2_I}{2\Delta^2})v_{I-1} + (\frac{\mu_I^-}{\Delta} - \frac{\mu_I^+}{\Delta}-\frac{\sigma^2_I}{\Delta^2})v_I + (\frac{\mu_I^+}{\Delta}+\frac{\sigma^2_I}{2\Delta^2})v_{I+1}$$
-
-Since $v'(\overline{x}) = 0$ then $v'(\overline{x})$ can be approximated since $v'(\overline{x}) = \frac{v_{I+1}-v_I}{\Delta} = 0$, then $v_{I+1} = v_I$. Then:
-
-$$\mathcal{A}v(x_I)\approx (-\frac{\mu_I^-}{\Delta}+\frac{\sigma^2_I}{2\Delta^2})v_{I-1} + [(\frac{\mu_I^-}{\Delta} - \frac{\mu_I^+}{\Delta}-\frac{\sigma^2_I}{\Delta^2}) + (\frac{\mu_I^+}{\Delta}+\frac{\sigma^2_I}{2\Delta^2})]v_{I}$$
-
-Thus the formulas for the discretization $A$ of $\mathcal{A}$ on the uniform grid $\{x_i\}_{i=1}^I$ are:
-
-- For $i\in[2,I-1]$, $\mathcal{A}v(x_i)\approx (-\frac{\mu_i^-}{\Delta}+\frac{\sigma^2_i}{2\Delta^2})v_{i-1} + (\frac{\mu_i^-}{\Delta} - \frac{\mu_i^+}{\Delta}-\frac{\sigma^2_i}{\Delta^2})v_i + (\frac{\mu_i^+}{\Delta}+\frac{\sigma^2_i}{2\Delta^2})v_{i+1}$
-
-- For $i=1$, $\mathcal{A}v(x_1)\approx [(-\frac{\mu_1^-}{\Delta}+\frac{\sigma^2_1}{2\Delta^2}) + (\frac{\mu_1^-}{\Delta} - \frac{\mu_1^+}{\Delta}-\frac{\sigma^2_1}{\Delta^2})]v_1 + (\frac{\mu_1^+}{\Delta}+\frac{\sigma^2_1}{2\Delta^2})v_{2}$
-
-- For $i=I$, $\mathcal{A}v(x_I)\approx (-\frac{\mu_I^-}{\Delta}+\frac{\sigma^2_I}{2\Delta^2})v_{I-1} + [(\frac{\mu_I^-}{\Delta} - \frac{\mu_I^+}{\Delta}-\frac{\sigma^2_I}{\Delta^2}) + (\frac{\mu_I^+}{\Delta}+\frac{\sigma^2_I}{2\Delta^2})]v_{I}$
-
-To succinctly construct $A$ in our code,  the discretization of $\mathcal{A}$ on the discretized grid $\{x_i\}_{i=1}^I$, we can define $X,Y,$ and $Z$ such that:
-
-- Let $X = -\frac{\mu^-}{\Delta}+\frac{\sigma^2}{2\Delta^2}$
-
-- Let $Y = -\frac{\mu^+}{\Delta}+\frac{\mu^-}{\Delta^2}-\frac{\sigma^2}{\Delta^2}$
-
-- Let $Z = \frac{\mu^+}{\Delta}+\frac{\sigma^2}{2\Delta^2}$
-
-Where $\mu^- \equiv \{\mu(x_i)^-\}_{i=1}^I$, $\mu^+ \equiv \{\mu(x_i)^+\}_{i=1}^I$, and $\sigma^2 \equiv \{\sigma(x_i)^2\}_{i=1}^I$
-
-Then $A$ can be constucted as a tridiagonal matrix:
+Treating $x= \{0.2,0.3,0.4\}$ as an array $x = [0.2, 0.3, 0.4], then
+considering the original Bellman equation, we obtain:
 
-$$A=\begin{bmatrix}
-Y_1+X_1 & Z_1 & 0 & \cdots & \cdots & \cdots & 0 \\
-X_2 & Y_2 & Z_2 & 0 & \ddots &  & \vdots \\
-0 & \ddots & \ddots & \ddots & \ddots & & \vdots \\
-\vdots & & \ddots & \ddots & \ddots & \ddots & \vdots \\
-\vdots & & & \ddots & X_{I-1} & Y_{I-1} & Z_{I-1} \\
-0 & \cdots & \cdots & \cdots & 0 & Y_{I} & Y_I+Z_{I} \\
-\end{bmatrix}$$
+$$
+\begin{equation}
+  \rho v(x) \approx u(x) +\mu L_{1-}v(x) +\frac{\sigma^2}{2}L_2v(x)
+\end{equation}
+$$
 
-It thus follows that $\mathcal{A}v(x) \approx Av(x)$ on the discretized grid $x\equiv\{x_{i=1}\}_{i=1}^I$
+Which further implies:
 
+$$
+\begin{equation}
+(\rho I - \mu L_{1-} - \frac{\sigma^2}{2}L_2)v(x) \approx u(x)
+\end{equation}
+$$
 
-### Example Problem
+Then setting:
 
-Using $\mu(\cdot)$ and $\sigma(\cdot)$ as defined earlier $\mu(x) = -0.01x$ and $\sigma(x) = 0.1x$ to define the stochastic difference equation
-
-$$dx_t = \mu(x_t)dt + \sigma(x_t)d\mathbb{W}_t$$
-
-we can consider the problem where:
-
-- Utility function: $u(x) = \log(x)$
-
-- Left Boundary: $\underline{x} = 0.01$
-
-- Right Boundary: $\overline{x} = 10$
-
-- Discount factor: $\rho = 0.05$
-
-Note that the choices of $\mu(\cdot)$, $\sigma(\cdot)$, and $u(\cdot)$ render this particular problem control-free. That is, the state variable $x_t$ determines the drift and variance of the stochastics process along with determining the utility without any control by the agent.
-
-We can then discretize $[\underline{x},\overline{x}]=[0.01,10]$ by increments of $0.001$.
+$$
+\begin{equation}
+L = \rho I - \mu\frac{1}{0.1}L_{1-} - \frac{\sigma^2}{2}L_2
+\end{equation}
+$$
 
 ```julia
-x_min = 0.01
-x_max = 10
-
-grid = x_min:0.001:x_max
-
-function μ(x)
-    return -0.01*x
-end
-
-function σ2(x)
-    return (0.1*x)^2
-end
+x = [0.2, 0.3, 0.4]
 
 function u(x)
     return log(x)
 end
 
-n = size(grid)[1]
-
-delta = grid[2]-grid[1] # Create the delta with respect to x grid
-delta_2 = delta^2
-
-function μ_p(x)
-    return max(μ(x),0)
-end
-
-function μ_m(x)
-    return min(μ(x),0)
-end
-
-X = (-μ_m.(grid)/delta) + (σ2.(grid)/(2*delta_2))
-Y = (-μ_p.(grid)/delta) + (μ_m.(grid)/delta) - (σ2.(grid)/delta_2)
-Z = (μ_p.(grid)/delta) + (σ2.(grid)/(2*delta_2))
+ρ = 0.05
+μ = -0.1
+σ = 0.1
 
-A_2 = Tridiagonal(X[2:n], Y, Z[1:n-1])
-A_2[1,1] = Y[1] + X[1]
-A_2[n,n] = Y[n] + Z[n]
+L = I*ρ - μ*L_1m - σ^2 / 2 * L_2
 ```
 
-### Solving Control-Free Time-Independent Linear Hamilton-Jacobi Bellman Equation
-
-In the simple case of an uncontrolled linear stochastic process with no optimal stopping, if an agent has payoff $u(x)$ with discounting $\rho>0$, then the following Hamilton-Jacobi Bellman equation determines the value
+We obtain the linear system:
 
-$$\rho v(x) = u(x) + \mathcal{A}v(x)$$
+$$
+\begin{equation}
+Lv(x) \approx u(x)
+\end{equation}
+$$
 
-subject to the appropriate boundary conditions. Discretizing the grid, let $u \equiv \{u(x_i)\}_{i=1}^I\in\mathbb{R}^I$, where $v\in\mathbb{R}^I$ is unknown, then the discretized version is,
+And thus we obtain:
 
-$$\rho v = u + Av$$
+$$
+\begin{equation}
+v(x) \approx L^{-1}u(x)
+\end{equation}
+$$
 
-Which we can solve algebraically for $v$ as the linear system:
-
-$$(\rho I -A)v = u$$
-
-Note that $A$ is a tridiagonal matrix, making it particularly convenient to solve numerically.
 
 ```julia
-ρ = 0.05
-
-v_2 = (ρ*I - A_2)\ u.(grid)
+v = L \ u.(x)
+```
 
-plt_u = plot(grid, u.(grid), label="u", xlabel="x", ylabel = "u(x)", lw=2, legend=:topleft)
-plt_v2 = plot(grid, v_2, label="v_2", xlabel="x", ylabel="v(x)", lw=2, legend=:topleft)
+### Using the SDO Package
 
-plot(plt_u,plt_v2)
-```
+We can solve the same problem using the simple differential operators package:
 
 ```julia
-using Test
+x̄ = range(0.1, 0.5, length = 5)
+x = interiornodes(x̄)
 
-# We can verify that the solution v obtained with the linear operator from the SimpleDifferentialOperators package
-# agrees with the solution v_2 obtained fromm manually discretizing the linear operator
+bc = (Reflecting(), Reflecting())
+L = I * ρ - μ*L₁₋(x̄, bc) - σ^2 / 2 * L₂(x̄, bc)
 
-@test v[1] ≈ v_2[1]
-@test v[500] ≈ v_2[500]
-@test v[999] ≈ v_2[999]
+v_SDO = L \ u.(x)
 ```

From aad0630870540f96aa1714a8634f000947d3fd77 Mon Sep 17 00:00:00 2001
From: Adam Jozefiak <adam.jozefiak@alumni.ubc.ca>
Date: Sat, 20 Apr 2019 15:26:37 -0700
Subject: [PATCH 3/6] FD_intro first revision

---
 docs/examples/FD_intro.jmd     | 458 +++++++++++++++++++++++++++++++++
 docs/examples/FD_intro_old.jmd | 298 +++++++++++++++++++++
 2 files changed, 756 insertions(+)
 create mode 100644 docs/examples/FD_intro_old.jmd

diff --git a/docs/examples/FD_intro.jmd b/docs/examples/FD_intro.jmd
index 2f1b19f..d891a8f 100644
--- a/docs/examples/FD_intro.jmd
+++ b/docs/examples/FD_intro.jmd
@@ -106,7 +106,9 @@ $$
 $$
 
 ```julia
+
 L_1m = 10*[0 0 0; -1 1 0; 0 -1 1]
+
 ```
 
 (Note, a similar derivation using the forward differences would yield an analagous
@@ -296,3 +298,459 @@ L = I * ρ - μ*L₁₋(x̄, bc) - σ^2 / 2 * L₂(x̄, bc)
 
 v_SDO = L \ u.(x)
 ```
+
+# Motivating Example
+
+This notebook is a gentle introduction to finite difference along with demonstrating the application of SimpleDifferntialOperators.jl package.
+
+As a simple example, let
+
+- $x\in [0.2, 0.4]$ be a state variable on a domain following the SDE
+
+$$
+\begin{equation}
+  dx_t = \mu dx_t + \sigma W_t
+\end{equation}
+$$
+
+where the variable $x_t$ is reflected at $x_{\text{min}} = 0.2$ and $x_{\text{max}}=0.4$. Let $\mu = -0.1$ and $\sigma = 0.01$.
+
+- The payoffs for state $x$ are a function $\widetilde{f}(x)$ defined on the domain $[0.2, 0.4]$. For this particular example, let $\widetilde{f}(x) = log(x)$
+
+- $\widetilde{v}(x)$ is the value of the stream of payoffs discounted at rate $r=0.05$
+
+```julia
+
+µ = -0.1
+σ = 0.01
+r = 0.05
+
+function f(x)
+    return log(x)
+end
+
+```
+
+Then, through standard arguments, the stationary Bellman equation along with boundary conditions is
+
+$$
+\begin{align}
+  r\widetilde{v}(x) &= \widetilde{f}(x) + \mu\partial_x\widetilde{v}(x) + \frac{\sigma^2}{2}\partial_{xx}\widetilde{v}(x) \\
+  \partial_x\tilde{v}(0.2) &= 0 \\
+  \partial_x\tilde{v}(0.4) &= 0
+\end{align}
+$$
+
+Mapping to the notation of (derivation document ...)
+
+$$
+\begin{align}
+\tilde{L} &\equiv r - \mu\partial_x - \frac{\sigma^2}{2}\partial_{xx} \\
+\tilde{B} & \equiv \begin{bmatrix}
+  \partial_x\lvert_{x=0.2} \\
+  \partial_x\lvert_{x=0.4} \\
+  \end{bmatrix} \\
+b &\equiv \begin{bmatrix}
+0 \\
+0 \\
+\end{bmatrix}
+\end{align}
+$$
+
+We will solve for $\widetilde{v}(x)$ on the domain by the use of finite difference. We will discretize the domain $[0.2, 0.4]$ using a finite size of $\Delta = 0.1$, and hence let $x = [0.2, 0.3, 0.4]^T$ will employ $\bar{x} = [0.1, 0.2, 0.3, 0.4, 0.5]^T$. Likewise,
+we will solve for $v \equiv \{\widetilde{v}(x_i)\}_{i=1}^3$ and will employ $\bar{v} \equiv \{\widetilde{v}(x_i)\}_{i=0}^4$. $x_0 = 0.1, x_1 = 0.2, \ldots, x_4 = 0.5$.
+
+```julia
+
+x̄ = range(0.1, 0.5, length = (5)) # extended grid
+x = interiornodes(x̄) # interior grid
+
+```
+
+Next, we will derive the stencils for the discretized $3\times 5$ operators. $L_{1-}$ (due to downdrift), $I$, $L_2$.
+
+# First Derivative Operators
+
+To discretize the $\widetilde{L}_1$ operators we can use a backward difference approximation
+
+$$
+\begin{equation}
+\widetilde{L}_1\widetilde{v}(x_i) = \partial_x\widetilde{v}(x_i)\approx \frac{\bar{v}_i-\bar{v}_{i-1}}{\Delta}, \text{ for } i=1,\ldots, 3
+\end{equation}
+$$
+
+And with forward differences
+
+$$
+\begin{equation}
+\widetilde{L}_1\widetilde{v}(x_i) = \partial_x\widetilde{v}(x_i)\approx \frac{\bar{v}_{i+1}-\bar{v}_{i}}{\Delta}, \text{ for } i=1,\ldots, 3
+\end{equation}
+$$
+
+In order to calculate the derivaters for all $i=1,\ldots, 3$ (i.e. in the interior) we can stack these up and apply to the extension $\bar{v}$
+
+$$
+\begin{equation}
+\{\partial_x\widetilde{v}(x_i)\}_{i=1}^3 \approx L_{1-}\cdot \bar{v}
+\end{equation}
+$$
+
+or,
+
+$$
+\begin{equation}
+\{\partial_x\widetilde{v}(x_i)\}_{i=1}^3 \approx L_{1+}\cdot \bar{v}
+\end{equation}
+$$
+
+Where we define $L_{1-}$ from applying to the $\bar{v}$ vectors for all $i=1,\ldots, 3$
+
+$$
+\begin{equation}
+L_{1-} = \frac{1}{\Delta}\begin{bmatrix}
+  -1 & 1 & 0 & 0 & 0 \\
+  0 & -1 & 1 & 0 & 0 \\
+  0 & 0 & -1 & 1 & 0
+  \end{bmatrix}
+  \end{equation}
+  $$
+
+And similarly define $L_{1+}$ from applying to the $\bar{v}$ vector for all $i=1,\ldots, 3$
+
+$$
+\begin{equation}
+L_{1+} = \frac{1}{\Delta}\begin{bmatrix}
+  0 & -1 & 1 & 0 & 0 \\
+  0 & 0 & -1 & 1 & 0 \\
+  0 & 0 & 0 & -1 & 1
+  \end{bmatrix}
+  \end{equation}
+  $$
+
+  Since $\mu < 0$, we will only utilize the downdrift and hence only $L_{1-}$
+
+```julia
+
+L_1m = 1/0.1*[-1 1 0 0 0; 0 -1 1 0 0; 0 0 -1 1 0]
+
+```
+
+# Second Derivative Operators
+
+To discretize the $\widetilde{L}_2$ second order operator, we can use central differences
+
+$$
+\begin{equation}
+\widetilde{L}_2\widetilde{v}(x_i) \equiv \partial_{xx}\widetilde{v}(x_i)\approx \frac{\bar{v}_{i+1}-2\bar{v}_i\bar{v}_{i-1}}{\Delta^2} \text{ for } i=1,\ldots, 3
+\end{equation}
+$$
+
+In order to calculate the derivatives for all $i=1, \ldots, 3$ (i.e the interior) we can stack these up and aplly to the extension $\bar{v}$
+
+$$
+\begin{equation}
+\{\partial_{xx}\widetilde{v}(x_i)\}_{i=1}^3 \approx L_{2}\cdot \bar{v}
+\end{equation}
+$$
+
+Where we define $L_2$ from applying to the $\bar{v}$ vector for all $i=1,\ldots 3$,
+
+$$
+\begin{equation}
+L_{2} = \frac{1}{\Delta^2}\begin{bmatrix}
+  1 & -2 & 1 & 0 & 0 \\
+  0 & 1 & -2 & 1 & 0 \\
+  0 & 0 & 1 & -2 & 1
+  \end{bmatrix}
+  \end{equation}
+  $$
+
+```julia
+
+L_2 = 1/0.1^2*[1 -2 1 0 0 ; 0 1 -2 1 0; 0 0 1 -2 1 ]
+
+```
+
+  Additionally, we will require the identity operator for compositions. For simplicity, define the identity operator as the 0-th
+  order operator $\widetilde{L}_0\equiv I$ so that $\widetilde{L}_0\widetilde{v}(x) = \widetilde{v}(x)$.
+
+  With this, the operator applied to the $\bar{v}$ vector is trivial
+
+  $$
+  \begin{equation}
+  \widetilde{L}_0\widetilde{v}(x_i) \equiv \widetilde{v}(x_i) = \bar{v}_i, \text{ for } i = 1,\ldots, 3
+  \end{equation}
+  $$
+
+  And stacking it up for all $i=1,\ldots 3$
+
+  $$
+  \begin{equation}
+    \{\widetilde{v}(x_i)\}_{i=1}^3 \approx L_0\cdot \bar{v}
+  \end{equation}
+  $$
+
+  Where
+
+  $$
+  \begin{equation}
+  L_{0} = \begin{bmatrix}
+    0 & 1 & 0 & 0 & 0 \\
+    0 & 0 & 1 & 0 & 0 \\
+    0 & 0 & 0 & 1 & 0
+    \end{bmatrix}
+    \end{equation}
+    $$
+
+```julia
+
+L_0 = [0 1 0 0 0 ; 0 0 1 0 0; 0 0 0 1 0]
+
+```
+
+  Now that we have the operators, we can compose them together to create the complete differential operator for dthe diffusion process $L$
+
+$$
+\begin{equation}
+  L = rL_0 - \mu L_{1-} - \frac{σ^2}{2}L_2
+\end{equation}
+$$
+
+```julia
+
+L = r*L_0 - µ*L_1m - σ^2/2*L_2
+
+```
+
+    $$
+    \begin{equation}
+    L = \begin{bmatrix}
+      0 & 1 & 0 & 0 & 0 \\
+      0 & 0 & 1 & 0 & 0 \\
+      0 & 0 & 0 & 1 & 0
+      \end{bmatrix}
+      \end{equation}
+      $$
+
+
+
+# Boundary Conditions
+
+We will additionally derive the stencil for the boundary operator $B$
+
+Given the reflecting boundary conditions, using forward difference and backward differences
+
+$$
+\begin{equation}
+0 = \partial_xv(x_4) \approx \frac{v(x_4)-v(x_3)}{\Delta}
+\end{equation}
+$$
+
+and
+
+$$
+\begin{equation}
+0 = \partial_xv(x_1) \approx \frac{v(0.1)-v(0.0)}{\Delta}
+\end{equation}
+$$
+
+which imply that $v(x_1) = v(x_0)$ and that $v(x_4) = v(x_3)$
+
+Stacking the two conditions above, it follows that
+
+$$
+\begin{equation}
+B = \begin{bmatrix}
+  -1 & 1 & 0 & 0 & 0 \\
+  0 & 0 & 0 & -1 & 1
+  \end{bmatrix}
+\end{equation}
+$$
+
+```julia
+
+B = [-1 1 0 0 0; 0 0 0 -1 1]
+
+```
+
+#Putting it Together
+
+Then we can stack the operators and affine terms and obtain the system
+
+$$
+\begin{equation}
+LB\bar{v} =
+\begin{bmatrix}
+L \\
+B
+\end{bmatrix}\bar{v} =
+\begin{bmatrix}
+f \\
+b
+\end{bmatrix}
+\end{equation}
+$$
+
+Solving the above system as an inverse problem for $\bar{v}$, we can then extract the interior nodes $\bar{v}$ in order to obtain $\{\widetilde{v}(x_i)\}_{i=1}^3.
+
+```julia
+
+LB = vcat(L, B)
+
+b = [0; 0]
+
+ab = vcat(f.(x), b)
+
+v̄ = LB \ ab
+
+v = interiornodes(v̄)
+
+```
+
+# Applying Boundary Conditions to Operators, the SimpleDifferentialOperators.jl Implementation
+
+Additionally, we could apply the boundary conditions directly to the $L$ operators and solve for the $v$ vector directly as an alternative way of solving the system of $M+2$ equations where we solved for $\bar{v}$ (where in the above example $M=3$).
+In fact, this is the exact same implementation that the SimpleDifferentialOperators.jl package utilizes and we will thus utilize the SimpleDifferentialOperators.jl package throughout this derivation.e
+
+Subtracting the first row of $L_{1-}$ by the first of $B$ multiplied by $\Delta^{-1}$ gives, with the corresponding row operation matrix R,
+
+$$
+\begin{equation}
+RL_{1-} = \frac{1}{\Delta}
+\begin{bmatrix}
+0 & 0 & 0 & 0 & 0 \\
+0 & -1 & 1 & 0  & 0 \\
+0 & 0 & -1 & 1 & 0
+\end{bmatrix}
+\end{equation}
+$$
+
+Note that there is no non-zero element in the first and last column for nodes on the boundaries. Hence, solving the corresponding extended system,
+
+$$
+\begin{equation}
+\begin{bmatrix}
+L_{1-} \\
+B
+\end{bmatrix}
+\bar{v} =
+\begin{bmatrix}
+f \\
+b
+\end{bmatrix}
+\end{equation}
+$$
+
+is identical as solving the following system
+
+$$
+\begin{align}
+R\begin{bmatrix}
+L_{1-} \\
+B
+\end{bmatrix}
+\bar{v} &=
+R\begin{bmatrix}
+f \\
+b
+\end{bmatrix} \\
+&= \begin{bmatrix}
+f \\
+b
+\end{bmatrix}
+\end{align}
+$$
+
+as $b$ is a zero vector so that the row operations $R$ do not change anything on the right hand side. Furthermore,
+ limited to the interior, solving $b$ in the system above is identical as solving the following system with an operator $L^B_{1-}$ on interior nodes:
+
+$$
+\begin{equation}
+L^{B}_{1-}v=f
+\end{equation}
+$$
+
+ where we have
+
+$$
+\begin{equation}
+L^B_{1-} = \frac{1}{\Delta}\begin{bmatrix}
+0 & 0 & 0 \\
+-1 & 1 & 0 \\
+0 & -1 & 1
+\end{bmatrix}
+\end{equation}
+$$
+
+instead of solving the full system without boundary conditions.
+
+```julia
+
+bc = (Reflecting(), Reflecting())
+L₁₋(x̄, bc)
+
+```
+
+For the $L_2$ operator, we can add the first row $L_2$ by the first row $B$ multiplied by $\Delta^{-2}$ and subtract the
+last row $L_2$ by the second row of $B$ multiplied by $\Delta^{-1}$, we have the following differential operator with the boundary condition $B$ applied for $L_2$
+
+$$
+\begin{equation}
+L^B_2 = \frac{1}{\Delta^2}
+\begin{bmatrix}
+-1 & 1 & 0 \\
+1 & -2 & 1 \\
+0 & -1 & 1
+\end{bmatrix}
+\end{equation}
+$$
+
+```julia
+
+L₂(x̄, bc)
+
+```
+
+Note that for $L_0$ there is no non-zero element in the first and last column for nodes on the boundaries and hence restricting $L_0$ to the interior
+columns we obtain the identity matrix $I$.
+
+Thus, setting
+
+$$
+\begin{equation}
+L^B = rI - \mu L^B_{1-} - \frac{\sigma^2}{2}L^B_2
+\end{equation}
+$$
+
+we note that solving
+
+$$
+\begin{equation}
+\begin{bmatrix}
+L \\
+B
+\end{bmatrix}\bar{v} =
+\begin{bmatrix}
+f \\
+b
+\end{bmatrix}
+\end{equation}
+$$
+
+is identical to solving
+
+$$
+\begin{equation}
+L^Bv = f
+\end{equation}
+$$
+
+```julia
+
+L = I * ρ - μ*L₁₋(x̄, bc) - σ^2 / 2 * L₂(x̄, bc)
+
+v_SDO = L \ u.(x)
+
+```
+
+We observe that this result is identical as the result obtained from solving the stacked system $LB$ above.
diff --git a/docs/examples/FD_intro_old.jmd b/docs/examples/FD_intro_old.jmd
new file mode 100644
index 0000000..2f1b19f
--- /dev/null
+++ b/docs/examples/FD_intro_old.jmd
@@ -0,0 +1,298 @@
+```julia
+using SimpleDifferentialOperators,Plots, LinearAlgebra
+```
+
+## Gentle Introduction to Upwind Finite Differences
+
+This notebook is a brief tutorial on upwind finite difference method for solving
+linear ordinary differential equations.
+
+In particular, this tutorial will be motivated by the following problem. Consider payoffs
+$u(x)$ for a stochastic process $x\in [0.2,04]$ with infitesimal generate $\mathcal{L}_x$
+discounted at rate $\rho > 0$. Then the payoffs can be valued by the
+Bellman equation:
+
+$$
+\begin{equation}
+  \rho v(x) = u(x) + \mathcal{L}_x v(x)
+\end{equation}
+$$
+
+We will consider the example $\mathcal{L}_x \equiv \mu \partial_x + \frac{\sigma^2}{2}\partial_{xx}$
+where $\mu = -0.1$ and $\sigma = 0.1$. We will consider $u(x) = \log(x)$ and we will impose the boundary
+conditions that $v'(0.2) = v'(0.4) = 0$.
+
+We will begin by first discretizing the interval $[0.2,0.4]$. Define $x = \{0.2, 0.3, 0.4\}$. Additionally,
+we will require the extended grid $\bar{x} = \{0.1,0.2,0.3,0.4,0.5\}$ in order to compute $\partial_x v(x)$ and
+$\partial_{xx}v(x)$ at the original boundaries $0.2$ and $0.4$.
+
+First, we will discretize $v'(x_i)$ for each $x_i\in x$.
+
+Note that by forward difference, for an arbitrary $x$ and $\Delta$:
+
+$$
+\begin{equation}
+  v'(x) \approx \frac{v(x+\Delta)-v(x)}{\Delta}
+\end{equation}
+$$
+
+Additionally, by backward difference,
+
+$$
+\begin{equation}
+  v'(x) \approx \frac{v(x)-v(x-\Delta)}{\Delta}
+\end{equation}
+$$
+
+Since in the stochastic process for $x$ the drift $\mu = -0.1 < 0$, we will
+approximate $v'(x)$ by the backward difference. In our instance, the grid $x$ has
+uniformly spaced grid points, so $\Delta = 0.1$ for each $x_i\in x$.
+
+So it follows that:
+
+$$
+\begin{equation}
+  v'(0.2) = 0 \text{ by the reflecting boundary condition}
+\end{equation}
+$$
+
+$$
+\begin{equation}
+  v'(0.3) \approx \frac{v(0.3)-v(0.2)}{0.1}
+\end{equation}
+$$
+
+$$
+\begin{equation}
+  v'(0.4) \approx \frac{v(0.4)-v(0.3)}{0.1}
+\end{equation}
+$$
+
+The above equations give us a linear system and hence we can represent them as follows:
+
+$$
+\begin{equation}
+  \begin{bmatrix}
+    v'(0.2) \\
+    v'(0.3) \\
+    v'(0.4) \\
+  \end{bmatrix}
+  \approx
+  \frac{1}{0.1}
+  \begin{bmatrix}
+    0 & 0 & 0 \\
+    -1 & 1 & 0 \\
+    0 & -1 & 1 \\
+  \end{bmatrix}
+  \begin{bmatrix}
+    v(0.2) \\
+    v(0.3) \\
+    v(0.4) \\
+  \end{bmatrix}
+\end{equation}
+$$
+
+Additionally, we will define $L_{1-}$ to be:
+
+$$
+\begin{equation}
+  L_{1-} = \frac{1}{0.1}
+  \begin{bmatrix}
+    0 & 0 & 0 \\
+    -1 & 1 & 0 \\
+    0 & -1 & 1 \\
+  \end{bmatrix}
+\end{equation}
+$$
+
+```julia
+L_1m = 10*[0 0 0; -1 1 0; 0 -1 1]
+```
+
+(Note, a similar derivation using the forward differences would yield an analagous
+matrix $L_{1+}$)
+
+In a similar manner, we will discretize $v''(x_i)$ for each $x_i \in x = \{0.2,0.3, 0.4\}$.
+
+Note, by central difference, for arbitrary $x$ and $\Delta$:
+
+$$
+\begin{equation}
+  v''(x) \approx \frac{v(x+\Delta)-2v(x) + v(x-\Delta)}{\Delta^2}
+\end{equation}
+$$
+
+So it follows that:
+
+$$
+\begin{equation}
+  v''(0.2) \approx \frac{v(0.3)-2v(0.2) + v(0.1)}{0.1^2}
+\end{equation}
+$$
+
+Note that by the boundary conditions and backward difference:
+
+$$
+\begin{equation}
+0 = v'(0.2) \approx \frac{v(0.2)-v(0.1)}{0.1}
+\end{equation}
+$$
+
+which implies that $v(0.2) = v(0.1)$.
+
+Therefore,
+
+$$
+\begin{equation}
+  v''(0.2) \approx \frac{v(0.3)-v(0.2)}{0.1^2}
+\end{equation}
+$$
+
+Then $0.3$:
+
+$$
+\begin{equation}
+  v''(0.3) \approx \frac{v(0.4)-2v(0.3) + v(0.2)}{0.1^2}
+\end{equation}
+$$
+
+For $0.4$:
+
+$$
+\begin{equation}
+  v''(0.4) \approx \frac{v(0.5)-2v(0.4) + v(0.3)}{0.1^2}
+\end{equation}
+$$
+
+Note that by the boundary conditions and forward difference:
+
+$$
+\begin{equation}
+0 = v'(0.4) \approx \frac{v(0.5)-v(0.4)}{0.1}
+\end{equation}
+$$
+
+which implies that $v(0.4) = v(0.5)$.
+
+Therefore,
+
+$$
+\begin{equation}
+  v''(0.4) \approx \frac{-v(0.4) + v(0.3)}{0.1^2}
+\end{equation}
+$$
+
+The above equations give us a linear system and hence we can represent them as follows:
+
+$$
+\begin{equation}
+  \begin{bmatrix}
+    v''(0.2) \\
+    v''(0.3) \\
+    v''(0.4) \\
+  \end{bmatrix}
+  \approx
+  \frac{1}{0.1^2}
+  \begin{bmatrix}
+    -1 & 1 & 0 \\
+    1 & -2 & 1 \\
+    0 & 1 & -1 \\
+  \end{bmatrix}
+  \begin{bmatrix}
+    v(0.2) \\
+    v(0.3) \\
+    v(0.4) \\
+  \end{bmatrix}
+\end{equation}
+$$
+
+Additionally, we will define $L_2$ to be:
+
+$$
+\begin{equation}
+  L_2 = \frac{1}{0.1^2}
+  \begin{bmatrix}
+    -1 & 1 & 0 \\
+    1 & -2 & 1 \\
+    0 & 1 & -1 \\
+  \end{bmatrix}
+\end{equation}
+$$
+
+```julia
+L_2 = 100*[-1 1 0; 1 -2 1; 0 1 -1]
+```
+
+Treating $x= \{0.2,0.3,0.4\}$ as an array $x = [0.2, 0.3, 0.4], then
+considering the original Bellman equation, we obtain:
+
+$$
+\begin{equation}
+  \rho v(x) \approx u(x) +\mu L_{1-}v(x) +\frac{\sigma^2}{2}L_2v(x)
+\end{equation}
+$$
+
+Which further implies:
+
+$$
+\begin{equation}
+(\rho I - \mu L_{1-} - \frac{\sigma^2}{2}L_2)v(x) \approx u(x)
+\end{equation}
+$$
+
+Then setting:
+
+$$
+\begin{equation}
+L = \rho I - \mu\frac{1}{0.1}L_{1-} - \frac{\sigma^2}{2}L_2
+\end{equation}
+$$
+
+```julia
+x = [0.2, 0.3, 0.4]
+
+function u(x)
+    return log(x)
+end
+
+ρ = 0.05
+μ = -0.1
+σ = 0.1
+
+L = I*ρ - μ*L_1m - σ^2 / 2 * L_2
+```
+
+We obtain the linear system:
+
+$$
+\begin{equation}
+Lv(x) \approx u(x)
+\end{equation}
+$$
+
+And thus we obtain:
+
+$$
+\begin{equation}
+v(x) \approx L^{-1}u(x)
+\end{equation}
+$$
+
+
+```julia
+v = L \ u.(x)
+```
+
+### Using the SDO Package
+
+We can solve the same problem using the simple differential operators package:
+
+```julia
+x̄ = range(0.1, 0.5, length = 5)
+x = interiornodes(x̄)
+
+bc = (Reflecting(), Reflecting())
+L = I * ρ - μ*L₁₋(x̄, bc) - σ^2 / 2 * L₂(x̄, bc)
+
+v_SDO = L \ u.(x)
+```

From 6097cc96fbc9e41d46be4de8832742931e2a22db Mon Sep 17 00:00:00 2001
From: Adam Jozefiak <adam.jozefiak@alumni.ubc.ca>
Date: Sat, 20 Apr 2019 16:29:03 -0700
Subject: [PATCH 4/6] April 20 edits

---
 docs/examples/FD_intro.jmd | 404 +++++--------------------------------
 1 file changed, 47 insertions(+), 357 deletions(-)

diff --git a/docs/examples/FD_intro.jmd b/docs/examples/FD_intro.jmd
index d891a8f..3ee3aef 100644
--- a/docs/examples/FD_intro.jmd
+++ b/docs/examples/FD_intro.jmd
@@ -2,310 +2,16 @@
 using SimpleDifferentialOperators,Plots, LinearAlgebra
 ```
 
-## Gentle Introduction to Upwind Finite Differences
+# Gentle Introduction to Upwind Finite Differences
 
-This notebook is a brief tutorial on upwind finite difference method for solving
-linear ordinary differential equations.
-
-In particular, this tutorial will be motivated by the following problem. Consider payoffs
-$u(x)$ for a stochastic process $x\in [0.2,04]$ with infitesimal generate $\mathcal{L}_x$
-discounted at rate $\rho > 0$. Then the payoffs can be valued by the
-Bellman equation:
-
-$$
-\begin{equation}
-  \rho v(x) = u(x) + \mathcal{L}_x v(x)
-\end{equation}
-$$
-
-We will consider the example $\mathcal{L}_x \equiv \mu \partial_x + \frac{\sigma^2}{2}\partial_{xx}$
-where $\mu = -0.1$ and $\sigma = 0.1$. We will consider $u(x) = \log(x)$ and we will impose the boundary
-conditions that $v'(0.2) = v'(0.4) = 0$.
-
-We will begin by first discretizing the interval $[0.2,0.4]$. Define $x = \{0.2, 0.3, 0.4\}$. Additionally,
-we will require the extended grid $\bar{x} = \{0.1,0.2,0.3,0.4,0.5\}$ in order to compute $\partial_x v(x)$ and
-$\partial_{xx}v(x)$ at the original boundaries $0.2$ and $0.4$.
-
-First, we will discretize $v'(x_i)$ for each $x_i\in x$.
-
-Note that by forward difference, for an arbitrary $x$ and $\Delta$:
-
-$$
-\begin{equation}
-  v'(x) \approx \frac{v(x+\Delta)-v(x)}{\Delta}
-\end{equation}
-$$
-
-Additionally, by backward difference,
-
-$$
-\begin{equation}
-  v'(x) \approx \frac{v(x)-v(x-\Delta)}{\Delta}
-\end{equation}
-$$
-
-Since in the stochastic process for $x$ the drift $\mu = -0.1 < 0$, we will
-approximate $v'(x)$ by the backward difference. In our instance, the grid $x$ has
-uniformly spaced grid points, so $\Delta = 0.1$ for each $x_i\in x$.
-
-So it follows that:
-
-$$
-\begin{equation}
-  v'(0.2) = 0 \text{ by the reflecting boundary condition}
-\end{equation}
-$$
-
-$$
-\begin{equation}
-  v'(0.3) \approx \frac{v(0.3)-v(0.2)}{0.1}
-\end{equation}
-$$
-
-$$
-\begin{equation}
-  v'(0.4) \approx \frac{v(0.4)-v(0.3)}{0.1}
-\end{equation}
-$$
-
-The above equations give us a linear system and hence we can represent them as follows:
-
-$$
-\begin{equation}
-  \begin{bmatrix}
-    v'(0.2) \\
-    v'(0.3) \\
-    v'(0.4) \\
-  \end{bmatrix}
-  \approx
-  \frac{1}{0.1}
-  \begin{bmatrix}
-    0 & 0 & 0 \\
-    -1 & 1 & 0 \\
-    0 & -1 & 1 \\
-  \end{bmatrix}
-  \begin{bmatrix}
-    v(0.2) \\
-    v(0.3) \\
-    v(0.4) \\
-  \end{bmatrix}
-\end{equation}
-$$
-
-Additionally, we will define $L_{1-}$ to be:
-
-$$
-\begin{equation}
-  L_{1-} = \frac{1}{0.1}
-  \begin{bmatrix}
-    0 & 0 & 0 \\
-    -1 & 1 & 0 \\
-    0 & -1 & 1 \\
-  \end{bmatrix}
-\end{equation}
-$$
-
-```julia
-
-L_1m = 10*[0 0 0; -1 1 0; 0 -1 1]
-
-```
-
-(Note, a similar derivation using the forward differences would yield an analagous
-matrix $L_{1+}$)
-
-In a similar manner, we will discretize $v''(x_i)$ for each $x_i \in x = \{0.2,0.3, 0.4\}$.
-
-Note, by central difference, for arbitrary $x$ and $\Delta$:
-
-$$
-\begin{equation}
-  v''(x) \approx \frac{v(x+\Delta)-2v(x) + v(x-\Delta)}{\Delta^2}
-\end{equation}
-$$
-
-So it follows that:
-
-$$
-\begin{equation}
-  v''(0.2) \approx \frac{v(0.3)-2v(0.2) + v(0.1)}{0.1^2}
-\end{equation}
-$$
-
-Note that by the boundary conditions and backward difference:
-
-$$
-\begin{equation}
-0 = v'(0.2) \approx \frac{v(0.2)-v(0.1)}{0.1}
-\end{equation}
-$$
-
-which implies that $v(0.2) = v(0.1)$.
-
-Therefore,
-
-$$
-\begin{equation}
-  v''(0.2) \approx \frac{v(0.3)-v(0.2)}{0.1^2}
-\end{equation}
-$$
-
-Then $0.3$:
-
-$$
-\begin{equation}
-  v''(0.3) \approx \frac{v(0.4)-2v(0.3) + v(0.2)}{0.1^2}
-\end{equation}
-$$
-
-For $0.4$:
-
-$$
-\begin{equation}
-  v''(0.4) \approx \frac{v(0.5)-2v(0.4) + v(0.3)}{0.1^2}
-\end{equation}
-$$
-
-Note that by the boundary conditions and forward difference:
-
-$$
-\begin{equation}
-0 = v'(0.4) \approx \frac{v(0.5)-v(0.4)}{0.1}
-\end{equation}
-$$
-
-which implies that $v(0.4) = v(0.5)$.
-
-Therefore,
-
-$$
-\begin{equation}
-  v''(0.4) \approx \frac{-v(0.4) + v(0.3)}{0.1^2}
-\end{equation}
-$$
-
-The above equations give us a linear system and hence we can represent them as follows:
-
-$$
-\begin{equation}
-  \begin{bmatrix}
-    v''(0.2) \\
-    v''(0.3) \\
-    v''(0.4) \\
-  \end{bmatrix}
-  \approx
-  \frac{1}{0.1^2}
-  \begin{bmatrix}
-    -1 & 1 & 0 \\
-    1 & -2 & 1 \\
-    0 & 1 & -1 \\
-  \end{bmatrix}
-  \begin{bmatrix}
-    v(0.2) \\
-    v(0.3) \\
-    v(0.4) \\
-  \end{bmatrix}
-\end{equation}
-$$
-
-Additionally, we will define $L_2$ to be:
-
-$$
-\begin{equation}
-  L_2 = \frac{1}{0.1^2}
-  \begin{bmatrix}
-    -1 & 1 & 0 \\
-    1 & -2 & 1 \\
-    0 & 1 & -1 \\
-  \end{bmatrix}
-\end{equation}
-$$
-
-```julia
-L_2 = 100*[-1 1 0; 1 -2 1; 0 1 -1]
-```
-
-Treating $x= \{0.2,0.3,0.4\}$ as an array $x = [0.2, 0.3, 0.4], then
-considering the original Bellman equation, we obtain:
-
-$$
-\begin{equation}
-  \rho v(x) \approx u(x) +\mu L_{1-}v(x) +\frac{\sigma^2}{2}L_2v(x)
-\end{equation}
-$$
-
-Which further implies:
-
-$$
-\begin{equation}
-(\rho I - \mu L_{1-} - \frac{\sigma^2}{2}L_2)v(x) \approx u(x)
-\end{equation}
-$$
-
-Then setting:
-
-$$
-\begin{equation}
-L = \rho I - \mu\frac{1}{0.1}L_{1-} - \frac{\sigma^2}{2}L_2
-\end{equation}
-$$
-
-```julia
-x = [0.2, 0.3, 0.4]
-
-function u(x)
-    return log(x)
-end
-
-ρ = 0.05
-μ = -0.1
-σ = 0.1
-
-L = I*ρ - μ*L_1m - σ^2 / 2 * L_2
-```
-
-We obtain the linear system:
-
-$$
-\begin{equation}
-Lv(x) \approx u(x)
-\end{equation}
-$$
-
-And thus we obtain:
-
-$$
-\begin{equation}
-v(x) \approx L^{-1}u(x)
-\end{equation}
-$$
-
-
-```julia
-v = L \ u.(x)
-```
-
-### Using the SDO Package
-
-We can solve the same problem using the simple differential operators package:
-
-```julia
-x̄ = range(0.1, 0.5, length = 5)
-x = interiornodes(x̄)
-
-bc = (Reflecting(), Reflecting())
-L = I * ρ - μ*L₁₋(x̄, bc) - σ^2 / 2 * L₂(x̄, bc)
-
-v_SDO = L \ u.(x)
-```
+This notebook is a gentle introduction to upwind finite differences along with demonstrating the application of the SimpleDifferntialOperators.jl package.
+This notebook follows the notation of the SimpleDifferntialOperators.jl package and serves as a motivating example for the package.
 
 # Motivating Example
 
-This notebook is a gentle introduction to finite difference along with demonstrating the application of SimpleDifferntialOperators.jl package.
-
 As a simple example, let
 
-- $x\in [0.2, 0.4]$ be a state variable on a domain following the SDE
+* $x\in [0.2, 0.4]$ be a state variable on a domain following the stochastic differential equation
 
 $$
 \begin{equation}
@@ -315,9 +21,9 @@ $$
 
 where the variable $x_t$ is reflected at $x_{\text{min}} = 0.2$ and $x_{\text{max}}=0.4$. Let $\mu = -0.1$ and $\sigma = 0.01$.
 
-- The payoffs for state $x$ are a function $\widetilde{f}(x)$ defined on the domain $[0.2, 0.4]$. For this particular example, let $\widetilde{f}(x) = log(x)$
+* The payoffs for state $x$ are a function $\widetilde{f}(x)$ defined on the domain $[0.2, 0.4]$. For this particular example, let $\widetilde{f}(x) = log(x)$
 
-- $\widetilde{v}(x)$ is the value of the stream of payoffs discounted at rate $r=0.05$
+* $\widetilde{v}(x)$ is the value of the stream of payoffs discounted at rate $r=0.05$
 
 ```julia
 
@@ -336,17 +42,17 @@ Then, through standard arguments, the stationary Bellman equation along with bou
 $$
 \begin{align}
   r\widetilde{v}(x) &= \widetilde{f}(x) + \mu\partial_x\widetilde{v}(x) + \frac{\sigma^2}{2}\partial_{xx}\widetilde{v}(x) \\
-  \partial_x\tilde{v}(0.2) &= 0 \\
-  \partial_x\tilde{v}(0.4) &= 0
+  \partial_x\widetilde{v}(0.2) &= 0 \\
+  \partial_x\widetilde{v}(0.4) &= 0
 \end{align}
 $$
 
-Mapping to the notation of (derivation document ...)
+Our goal is to solve for $\widetilde{v}(x)$ on the domain $[0.2,0.4]$. Then, mapping to the notation of (derivation document ...)
 
 $$
 \begin{align}
-\tilde{L} &\equiv r - \mu\partial_x - \frac{\sigma^2}{2}\partial_{xx} \\
-\tilde{B} & \equiv \begin{bmatrix}
+\widetilde{L} &\equiv r - \mu\partial_x - \frac{\sigma^2}{2}\partial_{xx} \\
+\widetilde{B} & \equiv \begin{bmatrix}
   \partial_x\lvert_{x=0.2} \\
   \partial_x\lvert_{x=0.4} \\
   \end{bmatrix} \\
@@ -357,8 +63,9 @@ b &\equiv \begin{bmatrix}
 \end{align}
 $$
 
-We will solve for $\widetilde{v}(x)$ on the domain by the use of finite difference. We will discretize the domain $[0.2, 0.4]$ using a finite size of $\Delta = 0.1$, and hence let $x = [0.2, 0.3, 0.4]^T$ will employ $\bar{x} = [0.1, 0.2, 0.3, 0.4, 0.5]^T$. Likewise,
-we will solve for $v \equiv \{\widetilde{v}(x_i)\}_{i=1}^3$ and will employ $\bar{v} \equiv \{\widetilde{v}(x_i)\}_{i=0}^4$. $x_0 = 0.1, x_1 = 0.2, \ldots, x_4 = 0.5$.
+We will solve for $\widetilde{v}(x)$ on the domain $[0.2, 0.4]$ by discretizing this domain (and in turn the solution $\widetilde{v}(x)$) and then utilizing an upwind finite difference method.
+We will discretize the domain $[0.2, 0.4]$ using a constant grid size of $\Delta = 0.1$, and hence defining $x = [0.2, 0.3, 0.4]^T$ the discretized grid, and defining $\bar{x} = [0.1, 0.2, 0.3, 0.4, 0.5]^T$, the extended discretized grid.
+Likewise, we will solve for $v \equiv \{\widetilde{v}(x_i)\}_{i=1}^3$ by directly solving the extended-grid solution $\bar{v} \equiv \{\widetilde{v}(x_i)\}_{i=0}^4$. Note, we label the grid points as $x_0 = 0.1, x_1 = 0.2, \ldots, x_4 = 0.5$.
 
 ```julia
 
@@ -367,7 +74,7 @@ x = interiornodes(x̄) # interior grid
 
 ```
 
-Next, we will derive the stencils for the discretized $3\times 5$ operators. $L_{1-}$ (due to downdrift), $I$, $L_2$.
+Next, we will derive the stencils for the discretized $3\times 5$ operators: $L_{1-}$ (due to downdrift), $I$, $L_2$.
 
 # First Derivative Operators
 
@@ -387,7 +94,7 @@ $$
 \end{equation}
 $$
 
-In order to calculate the derivaters for all $i=1,\ldots, 3$ (i.e. in the interior) we can stack these up and apply to the extension $\bar{v}$
+In order to calculate the derivatives for all $i=1,\ldots, 3$ (i.e. in the interior) we can stack these up and apply to the extension $\bar{v}$
 
 $$
 \begin{equation}
@@ -411,9 +118,9 @@ L_{1-} = \frac{1}{\Delta}\begin{bmatrix}
   -1 & 1 & 0 & 0 & 0 \\
   0 & -1 & 1 & 0 & 0 \\
   0 & 0 & -1 & 1 & 0
-  \end{bmatrix}
-  \end{equation}
-  $$
+\end{bmatrix}
+\end{equation}
+$$
 
 And similarly define $L_{1+}$ from applying to the $\bar{v}$ vector for all $i=1,\ldots, 3$
 
@@ -423,11 +130,11 @@ L_{1+} = \frac{1}{\Delta}\begin{bmatrix}
   0 & -1 & 1 & 0 & 0 \\
   0 & 0 & -1 & 1 & 0 \\
   0 & 0 & 0 & -1 & 1
-  \end{bmatrix}
-  \end{equation}
-  $$
+\end{bmatrix}
+\end{equation}
+$$
 
-  Since $\mu < 0$, we will only utilize the downdrift and hence only $L_{1-}$
+  Since $\mu < 0$ throughout the domain, we will only utilize the downdrift and hence only $L_{1-}$
 
 ```julia
 
@@ -445,7 +152,7 @@ $$
 \end{equation}
 $$
 
-In order to calculate the derivatives for all $i=1, \ldots, 3$ (i.e the interior) we can stack these up and aplly to the extension $\bar{v}$
+In order to calculate the derivatives for all $i=1, \ldots, 3$ (i.e the interior) we can stack these up and apply to the extension $\bar{v}$
 
 $$
 \begin{equation}
@@ -453,7 +160,7 @@ $$
 \end{equation}
 $$
 
-Where we define $L_2$ from applying to the $\bar{v}$ vector for all $i=1,\ldots 3$,
+Where we define $L_2$ from applying to the $\bar{v}$ vector for all $i=1,\ldots, 3$,
 
 $$
 \begin{equation}
@@ -462,14 +169,15 @@ L_{2} = \frac{1}{\Delta^2}\begin{bmatrix}
   0 & 1 & -2 & 1 & 0 \\
   0 & 0 & 1 & -2 & 1
   \end{bmatrix}
-  \end{equation}
-  $$
+\end{equation}
+$$
 
 ```julia
 
 L_2 = 1/0.1^2*[1 -2 1 0 0 ; 0 1 -2 1 0; 0 0 1 -2 1 ]
 
 ```
+# Identity Operator
 
   Additionally, we will require the identity operator for compositions. For simplicity, define the identity operator as the 0-th
   order operator $\widetilde{L}_0\equiv I$ so that $\widetilde{L}_0\widetilde{v}(x) = \widetilde{v}(x)$.
@@ -482,7 +190,7 @@ L_2 = 1/0.1^2*[1 -2 1 0 0 ; 0 1 -2 1 0; 0 0 1 -2 1 ]
   \end{equation}
   $$
 
-  And stacking it up for all $i=1,\ldots 3$
+  And stacking it up for all $i=1,\ldots, 3$
 
   $$
   \begin{equation}
@@ -498,17 +206,17 @@ L_2 = 1/0.1^2*[1 -2 1 0 0 ; 0 1 -2 1 0; 0 0 1 -2 1 ]
     0 & 1 & 0 & 0 & 0 \\
     0 & 0 & 1 & 0 & 0 \\
     0 & 0 & 0 & 1 & 0
-    \end{bmatrix}
-    \end{equation}
-    $$
+\end{bmatrix}
+\end{equation}
+$$
 
 ```julia
 
-L_0 = [0 1 0 0 0 ; 0 0 1 0 0; 0 0 0 1 0]
+L_0 = [0 1 0 0 0; 0 0 1 0 0; 0 0 0 1 0]
 
 ```
 
-  Now that we have the operators, we can compose them together to create the complete differential operator for dthe diffusion process $L$
+  Now that we have the $L_{1-}$, $L_2$, and $L_0$ operators, we can compose them together to construct the discretized differential operator $L$ for the diffusion process
 
 $$
 \begin{equation}
@@ -522,18 +230,6 @@ L = r*L_0 - µ*L_1m - σ^2/2*L_2
 
 ```
 
-    $$
-    \begin{equation}
-    L = \begin{bmatrix}
-      0 & 1 & 0 & 0 & 0 \\
-      0 & 0 & 1 & 0 & 0 \\
-      0 & 0 & 0 & 1 & 0
-      \end{bmatrix}
-      \end{equation}
-      $$
-
-
-
 # Boundary Conditions
 
 We will additionally derive the stencil for the boundary operator $B$
@@ -563,7 +259,7 @@ $$
 B = \begin{bmatrix}
   -1 & 1 & 0 & 0 & 0 \\
   0 & 0 & 0 & -1 & 1
-  \end{bmatrix}
+\end{bmatrix}
 \end{equation}
 $$
 
@@ -573,13 +269,12 @@ B = [-1 1 0 0 0; 0 0 0 -1 1]
 
 ```
 
-#Putting it Together
+#Solving the System
 
-Then we can stack the operators and affine terms and obtain the system
+Now we can stack the operators and the affine terms and obtain the following $5\times 5$ system
 
 $$
 \begin{equation}
-LB\bar{v} =
 \begin{bmatrix}
 L \\
 B
@@ -591,17 +286,12 @@ b
 \end{equation}
 $$
 
-Solving the above system as an inverse problem for $\bar{v}$, we can then extract the interior nodes $\bar{v}$ in order to obtain $\{\widetilde{v}(x_i)\}_{i=1}^3.
+Solving the above system as an inverse problem for $\bar{v}$, we can then extract the interior nodes $\bar{v}$ in order to obtain $v=\{\widetilde{v}(x_i)\}_{i=1}^3$, the
+discretization of $\widetilde{v}$ on $[0.2,0.4]$.
 
 ```julia
 
-LB = vcat(L, B)
-
-b = [0; 0]
-
-ab = vcat(f.(x), b)
-
-v̄ = LB \ ab
+v̄ = [L; B] \ [f.(x); 0 ; 0]
 
 v = interiornodes(v̄)
 
@@ -609,8 +299,8 @@ v = interiornodes(v̄)
 
 # Applying Boundary Conditions to Operators, the SimpleDifferentialOperators.jl Implementation
 
-Additionally, we could apply the boundary conditions directly to the $L$ operators and solve for the $v$ vector directly as an alternative way of solving the system of $M+2$ equations where we solved for $\bar{v}$ (where in the above example $M=3$).
-In fact, this is the exact same implementation that the SimpleDifferentialOperators.jl package utilizes and we will thus utilize the SimpleDifferentialOperators.jl package throughout this derivation.e
+Additionally, we can apply the boundary conditions directly to the $L$ operators and solve for the $v$ vector directly as an alternative way of solving the system of $5\times 5$ system which directly gave us $\bar{v}$.
+In fact, this is the exact same implementation that the SimpleDifferentialOperators.jl package utilizes and we will thus utilize the SimpleDifferentialOperators.jl package throughout this derivation.
 
 Subtracting the first row of $L_{1-}$ by the first of $B$ multiplied by $\Delta^{-1}$ gives, with the corresponding row operation matrix R,
 
@@ -661,7 +351,7 @@ b
 \end{align}
 $$
 
-as $b$ is a zero vector so that the row operations $R$ do not change anything on the right hand side. Furthermore,
+as $b$ is a zero vector and hence the row operations $R$ do not change anything on the right hand side. Furthermore,
  limited to the interior, solving $b$ in the system above is identical as solving the following system with an operator $L^B_{1-}$ on interior nodes:
 
 $$
@@ -692,7 +382,7 @@ L₁₋(x̄, bc)
 ```
 
 For the $L_2$ operator, we can add the first row $L_2$ by the first row $B$ multiplied by $\Delta^{-2}$ and subtract the
-last row $L_2$ by the second row of $B$ multiplied by $\Delta^{-1}$, we have the following differential operator with the boundary condition $B$ applied for $L_2$
+last row $L_2$ by the second row of $B$ multiplied by $\Delta^{-2}$, we have the following differential operator with the boundary condition $B$ applied for $L_2$
 
 $$
 \begin{equation}
@@ -747,10 +437,10 @@ $$
 
 ```julia
 
-L = I * ρ - μ*L₁₋(x̄, bc) - σ^2 / 2 * L₂(x̄, bc)
+L = r*I - μ*L₁₋(x̄, bc) - σ^2 / 2 * L₂(x̄, bc)
 
-v_SDO = L \ u.(x)
+v_SDO = L \ f.(x)
 
 ```
 
-We observe that this result is identical as the result obtained from solving the stacked system $LB$ above.
+We observe that this result is identical to the result obtained from solving the stacked system above.

From 85bcd9bb4b2e7fc0f73bbb81c5d1653bdaed59c9 Mon Sep 17 00:00:00 2001
From: Adam Jozefiak <adam.jozefiak@alumni.ubc.ca>
Date: Fri, 26 Apr 2019 13:36:50 -0700
Subject: [PATCH 5/6] April 26

---
 docs/examples/FD_intro.jmd | 42 ++++++++++++++++++++------------------
 1 file changed, 22 insertions(+), 20 deletions(-)

diff --git a/docs/examples/FD_intro.jmd b/docs/examples/FD_intro.jmd
index 3ee3aef..1a7ce95 100644
--- a/docs/examples/FD_intro.jmd
+++ b/docs/examples/FD_intro.jmd
@@ -1,10 +1,10 @@
 ```julia
-using SimpleDifferentialOperators,Plots, LinearAlgebra
+using SimpleDifferentialOperators, Plots, LinearAlgebra
 ```
 
 # Gentle Introduction to Upwind Finite Differences
 
-This notebook is a gentle introduction to upwind finite differences along with demonstrating the application of the SimpleDifferntialOperators.jl package.
+This notebook is a gentle introduction to upwind finite difference methods along with demonstrating the application of the SimpleDifferntialOperators.jl package.
 This notebook follows the notation of the SimpleDifferntialOperators.jl package and serves as a motivating example for the package.
 
 # Motivating Example
@@ -47,7 +47,7 @@ $$
 \end{align}
 $$
 
-Our goal is to solve for $\widetilde{v}(x)$ on the domain $[0.2,0.4]$. Then, mapping to the notation of (derivation document ...)
+Our goal is to solve for $\widetilde{v}(x)$ on the domain $[0.2,0.4]$. Then, mapping to the PDE notation in the Derivation and Applications for SimpleDifferentialOperators.jl document.
 
 $$
 \begin{align}
@@ -74,7 +74,7 @@ x = interiornodes(x̄) # interior grid
 
 ```
 
-Next, we will derive the stencils for the discretized $3\times 5$ operators: $L_{1-}$ (due to downdrift), $I$, $L_2$.
+Next, we will derive the stencils for the discretized $3\times 5$ operators: $L_{1-}$, $L_{1+}$, $I$, and $L_2$.
 
 # First Derivative Operators
 
@@ -110,7 +110,7 @@ $$
 \end{equation}
 $$
 
-Where we define $L_{1-}$ from applying to the $\bar{v}$ vectors for all $i=1,\ldots, 3$
+Where we define $L_{1-}$ from applying the backward difference approximation to the $\bar{v}$ vectors for all $i=1,\ldots, 3$
 
 $$
 \begin{equation}
@@ -122,7 +122,7 @@ L_{1-} = \frac{1}{\Delta}\begin{bmatrix}
 \end{equation}
 $$
 
-And similarly define $L_{1+}$ from applying to the $\bar{v}$ vector for all $i=1,\ldots, 3$
+And similarly define $L_{1+}$ from applying the forward difference approximation to the $\bar{v}$ vector for all $i=1,\ldots, 3$
 
 $$
 \begin{equation}
@@ -134,7 +134,7 @@ L_{1+} = \frac{1}{\Delta}\begin{bmatrix}
 \end{equation}
 $$
 
-  Since $\mu < 0$ throughout the domain, we will only utilize the downdrift and hence only $L_{1-}$
+  Since $\mu < 0$ throughout the domain, by the upwind scheme, we will utilize $L_{1-}$
 
 ```julia
 
@@ -160,7 +160,7 @@ $$
 \end{equation}
 $$
 
-Where we define $L_2$ from applying to the $\bar{v}$ vector for all $i=1,\ldots, 3$,
+Where we define $L_2$ from applying the central difference approximation to the $\bar{v}$ vector for all $i=1,\ldots, 3$,
 
 $$
 \begin{equation}
@@ -174,7 +174,7 @@ $$
 
 ```julia
 
-L_2 = 1/0.1^2*[1 -2 1 0 0 ; 0 1 -2 1 0; 0 0 1 -2 1 ]
+L_2 = 1/0.1^2*[1 -2 1 0 0 ; 0 1 -2 1 0; 0 0 1 -2 1]
 
 ```
 # Identity Operator
@@ -190,7 +190,7 @@ L_2 = 1/0.1^2*[1 -2 1 0 0 ; 0 1 -2 1 0; 0 0 1 -2 1 ]
   \end{equation}
   $$
 
-  And stacking it up for all $i=1,\ldots, 3$
+  And stacking these up for all $i=1,\ldots, 3$
 
   $$
   \begin{equation}
@@ -232,7 +232,7 @@ L = r*L_0 - µ*L_1m - σ^2/2*L_2
 
 # Boundary Conditions
 
-We will additionally derive the stencil for the boundary operator $B$
+Next, we will derive the stencil for the boundary operator $B$
 
 Given the reflecting boundary conditions, using forward difference and backward differences
 
@@ -269,7 +269,7 @@ B = [-1 1 0 0 0; 0 0 0 -1 1]
 
 ```
 
-#Solving the System
+# Solving the System
 
 Now we can stack the operators and the affine terms and obtain the following $5\times 5$ system
 
@@ -286,7 +286,7 @@ b
 \end{equation}
 $$
 
-Solving the above system as an inverse problem for $\bar{v}$, we can then extract the interior nodes $\bar{v}$ in order to obtain $v=\{\widetilde{v}(x_i)\}_{i=1}^3$, the
+Solving the above system as an inverse problem for $\bar{v}$, we can then extract the interior nodes of $\bar{v}$ in order to obtain $v=\{\widetilde{v}(x_i)\}_{i=1}^3$, the
 discretization of $\widetilde{v}$ on $[0.2,0.4]$.
 
 ```julia
@@ -299,10 +299,10 @@ v = interiornodes(v̄)
 
 # Applying Boundary Conditions to Operators, the SimpleDifferentialOperators.jl Implementation
 
-Additionally, we can apply the boundary conditions directly to the $L$ operators and solve for the $v$ vector directly as an alternative way of solving the system of $5\times 5$ system which directly gave us $\bar{v}$.
-In fact, this is the exact same implementation that the SimpleDifferentialOperators.jl package utilizes and we will thus utilize the SimpleDifferentialOperators.jl package throughout this derivation.
+Additionally, we can apply the boundary conditions directly to the operators that compose $L$ and solve for the $v$ vector directly as opposed to solving for $\bar{v}$ directly with the stacked $5\times 5$ system.
+In fact, this is the implementation of the SimpleDifferentialOperators.jl package and we will thus utilize the SimpleDifferentialOperators.jl package throughout this derivation.
 
-Subtracting the first row of $L_{1-}$ by the first of $B$ multiplied by $\Delta^{-1}$ gives, with the corresponding row operation matrix R,
+We note that subtracting the first row of $L_{1-}$ by the first of $B$ multiplied by $\Delta^{-1}$ gives, with the corresponding row operation matrix R,
 
 $$
 \begin{equation}
@@ -352,7 +352,7 @@ b
 $$
 
 as $b$ is a zero vector and hence the row operations $R$ do not change anything on the right hand side. Furthermore,
- limited to the interior, solving $b$ in the system above is identical as solving the following system with an operator $L^B_{1-}$ on interior nodes:
+ limited to the interior, solving $v$ in the system above is identical to solving the following system with an operator $L^B_{1-}$ on interior nodes:
 
 $$
 \begin{equation}
@@ -372,7 +372,7 @@ L^B_{1-} = \frac{1}{\Delta}\begin{bmatrix}
 \end{equation}
 $$
 
-instead of solving the full system without boundary conditions.
+instead of solving the full system with boundary conditions.
 
 ```julia
 
@@ -381,8 +381,9 @@ L₁₋(x̄, bc)
 
 ```
 
-For the $L_2$ operator, we can add the first row $L_2$ by the first row $B$ multiplied by $\Delta^{-2}$ and subtract the
-last row $L_2$ by the second row of $B$ multiplied by $\Delta^{-2}$, we have the following differential operator with the boundary condition $B$ applied for $L_2$
+Likewise, for the $L_2$ operator, we can obtain the differential operator $L_2^B$ on the interior nodes with the boundary condition $B$ applied to $L_2^B$ by the following procedure.
+Add the first row $L_2$ by the first row $B$ multiplied by $\Delta^{-2}$ and subtract the
+last row $L_2$ by the second row of $B$ multiplied by $\Delta^{-2}$, and restricting to the interior columns we obtain
 
 $$
 \begin{equation}
@@ -435,6 +436,7 @@ L^Bv = f
 \end{equation}
 $$
 
+
 ```julia
 
 L = r*I - μ*L₁₋(x̄, bc) - σ^2 / 2 * L₂(x̄, bc)

From 415aab97c3d7b42144cbbfd1d6267d54ca8d463f Mon Sep 17 00:00:00 2001
From: Adam Jozefiak <adam.jozefiak@alumni.ubc.ca>
Date: Fri, 26 Apr 2019 13:43:19 -0700
Subject: [PATCH 6/6] remove old notebook

---
 docs/examples/FD_intro_old.jmd | 298 ---------------------------------
 1 file changed, 298 deletions(-)
 delete mode 100644 docs/examples/FD_intro_old.jmd

diff --git a/docs/examples/FD_intro_old.jmd b/docs/examples/FD_intro_old.jmd
deleted file mode 100644
index 2f1b19f..0000000
--- a/docs/examples/FD_intro_old.jmd
+++ /dev/null
@@ -1,298 +0,0 @@
-```julia
-using SimpleDifferentialOperators,Plots, LinearAlgebra
-```
-
-## Gentle Introduction to Upwind Finite Differences
-
-This notebook is a brief tutorial on upwind finite difference method for solving
-linear ordinary differential equations.
-
-In particular, this tutorial will be motivated by the following problem. Consider payoffs
-$u(x)$ for a stochastic process $x\in [0.2,04]$ with infitesimal generate $\mathcal{L}_x$
-discounted at rate $\rho > 0$. Then the payoffs can be valued by the
-Bellman equation:
-
-$$
-\begin{equation}
-  \rho v(x) = u(x) + \mathcal{L}_x v(x)
-\end{equation}
-$$
-
-We will consider the example $\mathcal{L}_x \equiv \mu \partial_x + \frac{\sigma^2}{2}\partial_{xx}$
-where $\mu = -0.1$ and $\sigma = 0.1$. We will consider $u(x) = \log(x)$ and we will impose the boundary
-conditions that $v'(0.2) = v'(0.4) = 0$.
-
-We will begin by first discretizing the interval $[0.2,0.4]$. Define $x = \{0.2, 0.3, 0.4\}$. Additionally,
-we will require the extended grid $\bar{x} = \{0.1,0.2,0.3,0.4,0.5\}$ in order to compute $\partial_x v(x)$ and
-$\partial_{xx}v(x)$ at the original boundaries $0.2$ and $0.4$.
-
-First, we will discretize $v'(x_i)$ for each $x_i\in x$.
-
-Note that by forward difference, for an arbitrary $x$ and $\Delta$:
-
-$$
-\begin{equation}
-  v'(x) \approx \frac{v(x+\Delta)-v(x)}{\Delta}
-\end{equation}
-$$
-
-Additionally, by backward difference,
-
-$$
-\begin{equation}
-  v'(x) \approx \frac{v(x)-v(x-\Delta)}{\Delta}
-\end{equation}
-$$
-
-Since in the stochastic process for $x$ the drift $\mu = -0.1 < 0$, we will
-approximate $v'(x)$ by the backward difference. In our instance, the grid $x$ has
-uniformly spaced grid points, so $\Delta = 0.1$ for each $x_i\in x$.
-
-So it follows that:
-
-$$
-\begin{equation}
-  v'(0.2) = 0 \text{ by the reflecting boundary condition}
-\end{equation}
-$$
-
-$$
-\begin{equation}
-  v'(0.3) \approx \frac{v(0.3)-v(0.2)}{0.1}
-\end{equation}
-$$
-
-$$
-\begin{equation}
-  v'(0.4) \approx \frac{v(0.4)-v(0.3)}{0.1}
-\end{equation}
-$$
-
-The above equations give us a linear system and hence we can represent them as follows:
-
-$$
-\begin{equation}
-  \begin{bmatrix}
-    v'(0.2) \\
-    v'(0.3) \\
-    v'(0.4) \\
-  \end{bmatrix}
-  \approx
-  \frac{1}{0.1}
-  \begin{bmatrix}
-    0 & 0 & 0 \\
-    -1 & 1 & 0 \\
-    0 & -1 & 1 \\
-  \end{bmatrix}
-  \begin{bmatrix}
-    v(0.2) \\
-    v(0.3) \\
-    v(0.4) \\
-  \end{bmatrix}
-\end{equation}
-$$
-
-Additionally, we will define $L_{1-}$ to be:
-
-$$
-\begin{equation}
-  L_{1-} = \frac{1}{0.1}
-  \begin{bmatrix}
-    0 & 0 & 0 \\
-    -1 & 1 & 0 \\
-    0 & -1 & 1 \\
-  \end{bmatrix}
-\end{equation}
-$$
-
-```julia
-L_1m = 10*[0 0 0; -1 1 0; 0 -1 1]
-```
-
-(Note, a similar derivation using the forward differences would yield an analagous
-matrix $L_{1+}$)
-
-In a similar manner, we will discretize $v''(x_i)$ for each $x_i \in x = \{0.2,0.3, 0.4\}$.
-
-Note, by central difference, for arbitrary $x$ and $\Delta$:
-
-$$
-\begin{equation}
-  v''(x) \approx \frac{v(x+\Delta)-2v(x) + v(x-\Delta)}{\Delta^2}
-\end{equation}
-$$
-
-So it follows that:
-
-$$
-\begin{equation}
-  v''(0.2) \approx \frac{v(0.3)-2v(0.2) + v(0.1)}{0.1^2}
-\end{equation}
-$$
-
-Note that by the boundary conditions and backward difference:
-
-$$
-\begin{equation}
-0 = v'(0.2) \approx \frac{v(0.2)-v(0.1)}{0.1}
-\end{equation}
-$$
-
-which implies that $v(0.2) = v(0.1)$.
-
-Therefore,
-
-$$
-\begin{equation}
-  v''(0.2) \approx \frac{v(0.3)-v(0.2)}{0.1^2}
-\end{equation}
-$$
-
-Then $0.3$:
-
-$$
-\begin{equation}
-  v''(0.3) \approx \frac{v(0.4)-2v(0.3) + v(0.2)}{0.1^2}
-\end{equation}
-$$
-
-For $0.4$:
-
-$$
-\begin{equation}
-  v''(0.4) \approx \frac{v(0.5)-2v(0.4) + v(0.3)}{0.1^2}
-\end{equation}
-$$
-
-Note that by the boundary conditions and forward difference:
-
-$$
-\begin{equation}
-0 = v'(0.4) \approx \frac{v(0.5)-v(0.4)}{0.1}
-\end{equation}
-$$
-
-which implies that $v(0.4) = v(0.5)$.
-
-Therefore,
-
-$$
-\begin{equation}
-  v''(0.4) \approx \frac{-v(0.4) + v(0.3)}{0.1^2}
-\end{equation}
-$$
-
-The above equations give us a linear system and hence we can represent them as follows:
-
-$$
-\begin{equation}
-  \begin{bmatrix}
-    v''(0.2) \\
-    v''(0.3) \\
-    v''(0.4) \\
-  \end{bmatrix}
-  \approx
-  \frac{1}{0.1^2}
-  \begin{bmatrix}
-    -1 & 1 & 0 \\
-    1 & -2 & 1 \\
-    0 & 1 & -1 \\
-  \end{bmatrix}
-  \begin{bmatrix}
-    v(0.2) \\
-    v(0.3) \\
-    v(0.4) \\
-  \end{bmatrix}
-\end{equation}
-$$
-
-Additionally, we will define $L_2$ to be:
-
-$$
-\begin{equation}
-  L_2 = \frac{1}{0.1^2}
-  \begin{bmatrix}
-    -1 & 1 & 0 \\
-    1 & -2 & 1 \\
-    0 & 1 & -1 \\
-  \end{bmatrix}
-\end{equation}
-$$
-
-```julia
-L_2 = 100*[-1 1 0; 1 -2 1; 0 1 -1]
-```
-
-Treating $x= \{0.2,0.3,0.4\}$ as an array $x = [0.2, 0.3, 0.4], then
-considering the original Bellman equation, we obtain:
-
-$$
-\begin{equation}
-  \rho v(x) \approx u(x) +\mu L_{1-}v(x) +\frac{\sigma^2}{2}L_2v(x)
-\end{equation}
-$$
-
-Which further implies:
-
-$$
-\begin{equation}
-(\rho I - \mu L_{1-} - \frac{\sigma^2}{2}L_2)v(x) \approx u(x)
-\end{equation}
-$$
-
-Then setting:
-
-$$
-\begin{equation}
-L = \rho I - \mu\frac{1}{0.1}L_{1-} - \frac{\sigma^2}{2}L_2
-\end{equation}
-$$
-
-```julia
-x = [0.2, 0.3, 0.4]
-
-function u(x)
-    return log(x)
-end
-
-ρ = 0.05
-μ = -0.1
-σ = 0.1
-
-L = I*ρ - μ*L_1m - σ^2 / 2 * L_2
-```
-
-We obtain the linear system:
-
-$$
-\begin{equation}
-Lv(x) \approx u(x)
-\end{equation}
-$$
-
-And thus we obtain:
-
-$$
-\begin{equation}
-v(x) \approx L^{-1}u(x)
-\end{equation}
-$$
-
-
-```julia
-v = L \ u.(x)
-```
-
-### Using the SDO Package
-
-We can solve the same problem using the simple differential operators package:
-
-```julia
-x̄ = range(0.1, 0.5, length = 5)
-x = interiornodes(x̄)
-
-bc = (Reflecting(), Reflecting())
-L = I * ρ - μ*L₁₋(x̄, bc) - σ^2 / 2 * L₂(x̄, bc)
-
-v_SDO = L \ u.(x)
-```