<jupyter_start><jupyter_text>Announcements- Problem Set 5 posted on D2L, due Oct 20. Stay on at the end of class for last-minute clarifications and discussion.- Problem Set 6 will be on finite difference methods. Due either Oct 27 or Nov 2 - so, *ask questions during today's class*!- __Outlook__: Finite difference methods for BVP, application to ODEs and PDEs. Then on to Monte Carlo Markov Chain methods.- Conference for Undergraduate Women in Physics: online event in 2021, [applications accepted until 10/25](https://www.aps.org/programs/women/cuwip/)- TIMESTEP on 10/21 at 5pm: Grad School Application Process [details + zoom registration](https://lavinia.as.arizona.edu/~timestep/)  Example BVP Problem (Lecture 17)We want to solve an ODE (PDE) that instead of having initial conditions is contained to an interval and has values at the edges of the interval.  This naturally comes about when we consider spatial problems.  One of the simplest cases for this is the Poisson problem in one-dimension$$   u''(x) = \frac{\text{d}^2 u}{\text{d} x^2} = f(x).$$Because there are two derivatives in the equation we need two boundary conditions to solve the PDE (really an ODE in this case) uniquely.  To start, we consider the following basic problem with __Dirichlet boundary conditions__$$\begin{aligned}    u''(x) = f(x) \quad \Omega = [a, b] \\    u(a) = \alpha \quad u(b) = \beta.\end{aligned}$$In Lecture 17, we solved this problem using a shooting method (iterating the initial values), effectively turning the boundary value problem into an initial value problem. Linear System Approach FormulationThe second approach we will consider involves the formation of a system of equations to solve based on finite difference approximations.  Again let's consider an example problem where$$    u_{xx} = f(x)$$with the initial conditions $u(a) = u_a$ and $u(b) = u_b$.  We will show in the finite difference discussion below that the second order centered difference approximation for the second derivative for a function $u(x)$ is$$    u_{xx} \approx \frac{u(x_{i-1}) - 2 u(x_i) + u(x_{i+1})}{\Delta x^2}.$$If we discretize the domain of the original BVP into $N$ points (not including the boundaries) such that$$    x_i = a + \frac{b - a}{N+1} \cdot i ~~~ \text{where} ~~~ i = 1, \ldots, N$$we can then write the finite difference approximation as a system of linear equations!  If for instance we take $N = 5$ then$$\begin{aligned}    (U_{xx})_1 &\approx \frac{U_a - 2 U_1 + U_2}{\Delta x^2} \\    (U_{xx})_2 &\approx \frac{U_1 - 2 U_2 + U_3}{\Delta x^2} \\    (U_{xx})_3 &\approx \frac{U_2 - 2 U_3 + U_4}{\Delta x^2} \\    (U_{xx})_4 &\approx \frac{U_3 - 2 U_4 + U_5}{\Delta x^2} \\    (U_{xx})_5 &\approx \frac{U_4 - 2 U_5 + U_b}{\Delta x^2} \\\end{aligned}$$where we have used $U_a = u(a)$ and $U_b = u(b)$ as the boundary conditions.  Using these approximations to the derivatives we can then write the ODE as$$    \frac{1}{\Delta x^2}\begin{bmatrix}    -2 &  1 &    &    &    \\     1 & -2 &  1 &    &    \\       &  1 & -2 &  1 &    \\       &    &  1 & -2 &  1 \\       &    &    &  1 & -2 \\    \end{bmatrix} \begin{bmatrix}        U_1 \\ U_2 \\ U_3 \\ U_4 \\ U_5    \end{bmatrix} =     \begin{bmatrix}        f(x_1) \\ f(x_2) \\ f(x_3) \\ f(x_4) \\ f(x_5) \\    \end{bmatrix}.$$ Boundary ConditionsThis does not include the boundary conditions though.  We can add these values easily for Dirichlet boundary conditions by sending the values we know to the $b$ vector:$$\begin{aligned}    \frac{U_a - 2 U_1 + U_2}{\Delta x^2} = f(x_1) &\Rightarrow& \frac{- 2 U_1 + U_2}{\Delta x^2} = f(x_1) - \frac{U_a}{\Delta x^2} \\    \frac{U_4 - 2 U_5 + U_b}{\Delta x^2} = f(x_1) &\Rightarrow& \frac{U_4 - 2 U_5}{\Delta x^2} = f(x_5) - \frac{U_b}{\Delta x^2}\end{aligned}$$so that final system looks like$$    \frac{1}{\Delta x^2} \begin{bmatrix}    -2 &  1 &    &    &    \\     1 & -2 &  1 &    &    \\       &  1 & -2 &  1 &    \\       &    &  1 & -2 &  1 \\       &    &    &  1 & -2 \\    \end{bmatrix} \begin{bmatrix}        U_1 \\ U_2 \\ U_3 \\ U_4 \\ U_5    \end{bmatrix} =     \begin{bmatrix}        f(x_1) - \frac{U_a}{\Delta x^2} \\ f(x_2) \\ f(x_3) \\ f(x_4) \\ f(x_5) - \frac{U_b}{\Delta x^2} \\    \end{bmatrix}.$$ ExampleWant to solve the BVP$$    u_{xx} = e^x, \quad x \in [0, 1] \quad \text{with} \quad u(0) = 0.0, \text{ and } u(1) = 3$$via the construction of a linear system of equations.\begin{align*}    u_{xx} &= e^x \\    u_x &= A + e^x \\    u &= Ax + B + e^x\\    u(0) &= B + 1 = 0 \Rightarrow B = -1 \\    u(1) &= A - 1 + e^{1} = 3 \Rightarrow A = 4 - e\\     ~\\    u(x) &= (4 - e) x - 1 + e^x\end{align*}<jupyter_code># Problem setup
a = 0.0
b = 1.0
u_a = 0.0
u_b = 3.0
f = lambda x: numpy.exp(x)
u_true = lambda x: (4.0 - numpy.exp(1.0)) * x - 1.0 + numpy.exp(x)

# Discretization
N = 10
x_bc = numpy.linspace(a, b, N + 2)
x = x_bc[1:-1]
delta_x = (b - a) / (N + 1)

# Construct matrix A
A = numpy.zeros((N, N))
diagonal = numpy.ones(N) / delta_x**2
A += numpy.diag(diagonal * -2.0, 0)
A += numpy.diag(diagonal[:-1], 1)
A += numpy.diag(diagonal[:-1], -1)

# Construct RHS
b = f(x)
b[0] -= u_a / delta_x**2
b[-1] -= u_b / delta_x**2

# Solve system
U = numpy.empty(N + 2)
U[0] = u_a
U[-1] = u_b
U[1:-1] = numpy.linalg.solve(A, b)

# Plot result
fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)
axes.plot(x_bc, U, 'o', label="Computed")
axes.plot(x_bc, u_true(x_bc), 'k', label="True")
axes.set_title("Solution to $u_{xx} = e^x$")
axes.set_xlabel("x")
axes.set_ylabel("u(x)")
axes.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Error AnalysisA natural question to ask given our approximation $U_i$ is how close this is to the true solution $u(x)$ at the grid points $x_i$.  To address this we will define the error $E$ as$$    E = U - \widehat{U}$$where $U$ is the vector of the approximate solution and $\widehat{U}$ is the vector composed of the $u(x_i)$ such that$$    \widehat{U_i} = u(x_i)$$This leaves $E$ as a vector still so often we ask the question how does the norm of $E$ behave given a particular $\Delta x$.  For the $\infty$-norm we would have$$    ||E||_\infty = \max_{1 \leq i \leq m} |E_i| = \max_{1 \leq i \leq m} |U_i - u(x_i)|$$If we can show that $||E||_\infty$ goes to zero as $\Delta x \rightarrow 0$ we can then claim that the approximate solution $U_i$ at any of the grid points $E_i \rightarrow 0$.  If we would like to use other norms we often define slightly modified versions of the norms that also contain the grid width $\Delta x$ where$$\begin{aligned}    ||E||_1 &= \Delta x \sum^m_{i=1} |E_i| \\    ||E||_2 &= \left( \Delta x \sum^m_{i=1} |E_i|^2 \right )^{1/2}\end{aligned}$$These are referred to as *grid function norms*.The $E$ defined above is known as the *global error*.  We will now discuss how $E$ behaves given other factors as we define later. Local Truncation ErrorThe *local truncation error* (LTE) can be defined by replacing the approximate solution $U_i$ by the approximate solution $u(x_i)$.  Since the algebraic equations are an approximation to the original BVP, we do not expect that the true solution will exactly satisfy these equations, this resulting difference is the LTE.For our one-dimensional finite difference approximation from above we have$$    \frac{1}{\Delta x^2} (U_{i+1} - 2 U_i + U_{i-1}) = f(x_i).$$Replacing $U_i$ with $u(x_i)$ in this equation leads to$$    \tau_i = \frac{1}{\Delta x^2} (u(x_{i+1}) - 2 u(x_i) + u(x_{i-1})) - f(x_i).$$In this form the LTE is not as useful but if we assume $u(x)$ is smooth we can repalce the $u(x_i)$ with their Taylor series counterparts, similar to what we do for finite differences (below).  The relevant Taylor series are$$    u(x_{i \pm 1}) = u(x_i) \pm u'(x_i) \Delta x + \frac{1}{2} u''(x_i) \Delta x^2 \pm \frac{1}{6} u'''(x_i) \Delta x^3 + \frac{1}{24} u^{(4)}(x_i) \Delta x^4 + \mathcal{O}(\Delta x^5)$$This leads to an expression for $\tau_i$ of$$\begin{aligned}    \tau_i &= \frac{1}{\Delta x^2} \left [u''(x_i) \Delta x^2 + \frac{1}{12} u^{(4)}(x_i) \Delta x^4 + \mathcal{O}(\Delta x^5) \right ] - f(x_i) \\    &= u''(x_i) + \frac{1}{12} u^{(4)}(x_i) \Delta x^2 + \mathcal{O}(\Delta x^4) - f(x_i) \\    &= \frac{1}{12} u^{(4)}(x_i) \Delta x^2 + \mathcal{O}(\Delta x^4)\end{aligned}$$where we note that the true solution would satisfy $u''(x) = f(x)$.As long as $ u^{(4)}(x_i) $ remains finite (smooth) we know that $\tau_i \rightarrow 0$ as $\Delta x \rightarrow 0$We can also write the vector of LTEs as$$    \tau = A \widehat{U} - F$$which implies$$    A\widehat{U} = F + \tau.$$Note that here it is the exact solution evaluated stencil, not the approximated function. Global ErrorWhat we really want to bound is the global error $E$.  To relate the global error and LTE we can substitute $E = U - \widehat{U}$ into our expression for the LTE to find$$    A E = -\tau.$$This means that the global error is the solution to the system of equations we defined for the approximation except with $\tau$ as the forcing function rather than $F$!This also implies that the global error $E$ can be thought of as an approximation to similar BVP as we started with where$$    e''(x) = -\tau(x) \quad \Omega = [0, 1] \\    e(0) = 0 \quad e(1) = 0.$$We can solve this ODE directly by integrating twice since to find to leading order$$\begin{aligned}    e(x) &\approx -\frac{1}{12} \Delta x^2 u''(x) + \frac{1}{12} \Delta x^2 (u''(0) + x (u''(1) - u''(0))) \\    &= \mathcal{O}(\Delta x^2) \rightarrow 0 \quad \text{as} \quad \Delta x \rightarrow 0.\end{aligned}$$ StabilityWe showed that the continuous analog to $E$, $e(x)$, does in fact go to zero as $\Delta x \rightarrow 0$ but what about $E$?  Instead of showing something based on $e(x)$ let's look back at the original system of equations for the global error$$    A_{\Delta x} E_{\Delta x} = - \tau_{\Delta x}$$where we now denote a particular realization of the system by the corresponding grid spacing $\Delta x$.  If we could invert $A_{\Delta x}$ we could compute $E_{\Delta x}$ directly.  Assuming that we can and taking an appropriate norm we find$$\begin{aligned}    E_{\Delta x} &= (A_{\Delta x})^{-1} \tau_{\Delta x} \\    ||E_{\Delta x}|| &= ||(A_{\Delta x})^{-1} \tau_{\Delta x}|| \\    & \leq ||(A_{\Delta x})^{-1} ||~|| \tau_{\Delta x}||\end{aligned}$$We know that $\tau_{\Delta x} \rightarrow 0$ as $\Delta x \rightarrow 0$ already for our example so if we can bound the norm of the matrix $(A_{\Delta x})^{-1}$ by some constant $C$ for sufficiently small $\Delta x$ we can then write a bound on the global error of$$    ||E_{\Delta x}|| \leq C ||\tau_{\Delta x}||$$demonstrating that $E_{\Delta x} \rightarrow 0 $ at least as fast as $\tau_{\Delta x} \rightarrow 0$.We can generalize this observation to all linear BVP problems by supposing that we have a finite difference approximation to a linear BVP of the form$$    A_{\Delta x} U_{\Delta x} = F_{\Delta x},$$where $\Delta x$ is the grid spacing.  We say the approximation is *stable* if $(A_{\Delta x})^{-1}$ exists $\forall \Delta x < \Delta x_0$ and there is a constant $C$ such that$$    ||(A_{\Delta x})^{-1}|| \leq C \quad \forall \Delta x < \Delta x_0.$$ ConsistencyA related and important idea for the discretization of any PDE is that it be consistent with the equation we are approximating.  If$$    ||\tau_{\Delta x}|| \rightarrow 0 \quad \text{as}\quad  \Delta x \rightarrow 0$$then we say an approximation is *consistent* with the differential equation. ConvergenceWe now have all the pieces to say something about the global error $E$.  A method is said to be *convergent* if$$    ||E_{\Delta x}|| \rightarrow 0 \quad \text{as} \quad \Delta x \rightarrow 0.$$If an approximation is both consistent ($||\tau_{\Delta x}|| \rightarrow 0$ as $\Delta x \rightarrow 0$) and stable ($||E_{\Delta x}|| \leq C ||\tau_{\Delta x}||$) then the approximation is convergent.We have only derived this in the case of linear BVPs but in fact these criteria for convergence are often found to be true for any finite difference approximation (and beyond for that matter).  This statement of convergence can also often be strengthened to say$$    \mathcal{O}(\Delta x^p) \text{ LTE } + \text{ stability } \Rightarrow \mathcal{O}(\Delta x^p) \text{ global error}.$$_It turns out the most difficult part of this process is usually the statement regarding stability. We will now discuss how to prove stability in the 2-norm for our simple example next._ Stability in the 2-NormRecalling our definition of stability, we need to show that for our previously defined $A$ that$$    (A_{\Delta x})^{-1}$$exists and$$    ||(A_{\Delta x})^{-1}|| \leq C \quad \forall \Delta x < \Delta x_0$$for some $C$.  How do we know that $(A_{\Delta x})^{-1}$ exists?We can assume now that $A$ is in fact invertible but can we bound the norm of the inverse?  Recall that the 2-norm of a symmetric matrix is equal to its spectral radius,$$    ||A||_2 = \rho(A) = \max_{1\leq p \leq m} |\lambda_p|.$$Since the inverse of $A$ is also symmetric the eigenvalues of $A^{-1}$ are the inverses of the eigenvalues of $A$ implying that$$    ||A^{-1}||_2 = \rho(A^{-1}) = \max_{1\leq p \leq m} \left| \frac{1}{\lambda_p} \right| = \frac{1}{\max_{1\leq p \leq m} \left| \lambda_p \right|}.$$If none of the $\lambda_p$ of $A$ are zero for sufficiently small $\Delta x$ and the rest are finite as $\Delta x \rightarrow 0$ we have shown the stability of the approximation.The eigenvalues of the matrix $A$ from above can be written as$$    \lambda_p = \frac{2}{\Delta x^2} (\cos(p \pi \Delta x) - 1)$$with the corresponding eigenvectors $v^p$ $$    v^p_j = \sin(p \pi j \Delta x)$$as the $j$th component with $j = 1, \ldots, m$. Compute the smallest eigenvalueIf we can show that the eigenvalues are away from the origin then we know $||A||_2$ will be bounded.  In this case the eigenvalues are negative so we need to show that they are always strictly less than zero.$$    \lambda_p = \frac{2}{\Delta x^2} (\cos(p \pi \Delta x) - 1)$$Use a Taylor series to get an idea of how this behaves with respect to $\Delta x$From these expressions we know that smallest eigenvalue is$$\begin{aligned}    \lambda_1 &= \frac{2}{\Delta x^2} (\cos(p \pi \Delta x) - 1) \\    &= \frac{2}{\Delta x^2} \left (-\frac{1}{2} p^2 \pi^2 \Delta x^2 + \frac{1}{24} p^4 \pi^4 \Delta x^4 + \mathcal{O}(\Delta^6) \right ) \\    &= -p^2 \pi^2 + \mathcal{O}(\Delta x^2).\end{aligned}$$Note that this also gives us an error bound as this eigenvalue also will also lead to the largest eigenvalue of the inverse matrix.  We can therefore say$$    ||E^{\Delta x}||_2 \leq ||(A^{\Delta x})^{-1}||_2 ||\tau^{\Delta x}||_2 \approx \frac{1}{\pi^2} ||\tau^{\Delta x}||_2.$$ Stability in the $\infty$-NormThe straight forward approach to show that $||E||_\infty \rightarrow 0$ as $\Delta x \rightarrow 0$ would be to use the matrix bound$$    ||E||_\infty \leq \frac{1}{\sqrt{\Delta x}} ||E||_2.$$For our example problem we showed that $||E||_2 = \mathcal{O}(\Delta x^2)$ so this implies that we at least know that $||E||_\infty = \mathcal{O}(\Delta x^{3/2})$.  This is unfortunate as we expect $||E||_\infty = \mathcal{O}(\Delta x^{2})$ due to the discretization.  In order to alleviate this problem let's go back and consider our definition of stability but this time consider the $\infty$-norm.  Finite DifferencesFinite differences are expressions that approximate derivatives of a function evaluated at a set of points, often called a *stencil*.  These expressions can come in many different flavors including types of stencils, order of accuracy, and order of derivatives.  In this lecture we will review the process of derivation, error analysis and application of finite differences. Derivation of Finite DifferencesConsider three different ways to define a derivative at a point $x_i$$$    u'(x_i) = \lim_{\Delta x \rightarrow 0} \left \{ \begin{aligned}         &\frac{u(x_i + \Delta x) - u(x_i)}{\Delta x} & \equiv D_+ u(x_i)\\        &\frac{u(x_i + \Delta x) - u(x_i - \Delta_x)}{2 \Delta x} & \equiv D_0 u(x_i)\\        &\frac{u(x_i) - u(x_i - \Delta_x)}{\Delta x} & \equiv D_- u(x_i).    \end{aligned} \right .$$![Approximations to $u'(x)$](./fd_basic.png)If instead of allowing $\Delta x \rightarrow 0$ we come up with an approximation to the slope $u'(x_i)$ and hence our definitions of derivatives can directly be seen as approximations to derivatives when $\Delta x$ is perhaps small but non-zero.We will now delve into a more systematic way to derive these approximations as well as find higher order accurate approximations, higher order derivative approximations, and understand the error associated with the approximations. Taylor-Series MethodsOne way to derive finite difference approximations can be computed by using the Taylor series and the method of undetermined coefficients.$$u(x) = u(x_n) + (x - x_n) u'(x_n) + \frac{(x - x_n)^2}{2!} u''(x_n) + \frac{(x - x_n)^3}{3!} u'''(x_n) + \mathcal{O}((x - x_n)^4)$$Say we want to derive the second order accurate, first derivative approximation, this requires the values $(x_{n+1}, u(x_{n+1}))$ and $(x_{n-1}, u(x_{n-1}))$.  We can express these values via our Taylor series approximation above as$$\begin{aligned}    u(x_{n+1}) &= u(x_n) + (x_{n+1} - x_n) u'(x_n) + \frac{(x_{n+1} - x_n)^2}{2!} u''(x_n) + \frac{(x_{n+1} - x_n)^3}{3!} u'''(x_n) + \mathcal{O}((x_{n+1} - x_n)^4) \\    &= u(x_n) + \Delta x u'(x_n) + \frac{\Delta x^2}{2!} u''(x_n) + \frac{\Delta x^3}{3!} u'''(x_n) + \mathcal{O}(\Delta x^4) \end{aligned}$$and $$\begin{aligned}    u(x_{n-1}) &= u(x_n) + (x_{n-1} - x_n) u'(x_n) + \frac{(x_{n-1} - x_n)^2}{2!} u''(x_n) + \frac{(x_{n-1} - x_n)^3}{3!} u'''(x_n) + \mathcal{O}((x_{n-1} - x_n)^4) \\&= u(x_n) - \Delta x u'(x_n) + \frac{\Delta x^2}{2!} u''(x_n) - \frac{\Delta x^3}{3!} f'''(x_n) + \mathcal{O}(\Delta x^4) \end{aligned}$$Now to find out how to combine these into an expression for the derivative we assume our approximation looks like$$u'(x_n) + R(x_n) = A u(x_{n+1}) + B u(x_n) + C u(x_{n-1})$$where $R(x_n)$ is our error.  Plugging in the Taylor series approximations we find$$u'(x_n) + R(x_n) = A \left ( u(x_n) + \Delta x u'(x_n) + \frac{\Delta x^2}{2!} u''(x_n) + \frac{\Delta x^3}{3!} u'''(x_n) + \mathcal{O}(\Delta x^4)\right ) + B u(x_n) + C \left ( u(x_n) - \Delta x u'(x_n) + \frac{\Delta x^2}{2!} u''(x_n) - \frac{\Delta x^3}{3!} u'''(x_n) + \mathcal{O}(\Delta x^4) \right )$$Since we want $R(x_n) = \mathcal{O}(\Delta x^2)$ we want all terms lower than this to disappear except for those multiplying $u'(x_n)$ as those should sum to 1 to give us our approximation.  Collecting the terms with common derivatives $u^{(k)}(x_n)$ together we get a series of expressions for the coefficients $A$, $B$, and $C$ based on the fact we want an approximation to $u'(x_n)$.  The $n=0$ terms collected are $A + B + C$ and are set to 0 as we want the $u(x_n)$ term to disappear$$\begin{aligned}    u(x_n): & \quad A + B + C = 0 \\    u'(x_n): & \quad A \Delta x - C \Delta x = 1  \\    u''(x_n): & \quad A \frac{\Delta x^2}{2} + C \frac{\Delta x^2}{2} = 0 \end{aligned}$$This last equation $\Rightarrow A = -C$, using this in the second equation gives $A = \frac{1}{2 \Delta x}$ and $C = -\frac{1}{2 \Delta x}$.  The first equation then leads to $B = 0$.  Putting this altogether then gives us our previous expression including an estimate for the error:$$u'(x_n) + R(x_n) = \frac{u(x_{n+1}) - u(x_{n-1})}{2 \Delta x} + \frac{1}{2 \Delta x} \frac{\Delta x^3}{3!} u'''(x_n) + \mathcal{O}(\Delta x^4) + \frac{1}{2 \Delta x} \frac{\Delta x^3}{3!} u'''(x_n) + \mathcal{O}(\Delta x^4) $$$$R(x_n) = \frac{\Delta x^2}{3!} u'''(x_n) + \mathcal{O}(\Delta x^3) = \mathcal{O}(\Delta x^2)$$<jupyter_code>import numpy
import matplotlib.pyplot as plt

f = lambda x: numpy.sin(x)
f_prime = lambda x: numpy.cos(x)

# Use uniform discretization
#finely gridded x for plotting the exact function
x = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, 100)

N = 10
#coarsly grained x_hat to demonstrate forward difference estimates
x_hat = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, N)
delta_x = x_hat[1] - x_hat[0]

# Compute forward difference using a loop
f_prime_hat = numpy.empty(x_hat.shape)
for i in range(N - 1):
    f_prime_hat[i] = (f(x_hat[i+1]) - f(x_hat[i])) / delta_x
f_prime_hat[-1] = (f(x_hat[i]) - f(x_hat[i-1])) / delta_x

# Vector based calculation
# f_prime_hat[:-1] = (f(x_hat[1:]) - f(x_hat[:-1])) / (delta_x)

# Use first-order differences for points at edge of domain
f_prime_hat[-1] = (f(x_hat[-1]) - f(x_hat[-2])) / delta_x  # Backward Difference at x_N

fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)

axes.plot(x, f_prime(x), 'k')
axes.plot(x_hat + 0.5 * delta_x, f_prime_hat, 'ro')
axes.set_xlim((x[0], x[-1]))
axes.set_ylim((-1.1, 1.1))
axes.set_xlabel("x")
axes.set_ylabel(r"$f'(x)$")

plt.show()<jupyter_output><empty_output><jupyter_text>Example: Second Order DerivativeUsing our Taylor series approach lets derive the second order accurate second derivative formula.  Again we will use the same points and the Taylor series centered at $x = x_n$ so we end up with the same expression as before:$$\begin{aligned}    u''(x_n) + R(x_n) &= \quad A \left ( u(x_n) + \Delta x u'(x_n) + \frac{\Delta x^2}{2!} u''(x_n) + \frac{\Delta x^3}{3!} u'''(x_n) + \frac{\Delta x^4}{4!} u^{(4)}(x_n) + \mathcal{O}(\Delta x^5)\right ) \\    &\quad+ B u(x_n) \\    &\quad+ C \left ( u(x_n) - \Delta x u'(x_n) + \frac{\Delta x^2}{2!} u''(x_n) - \frac{\Delta x^3}{3!} u'''(x_n) + \frac{\Delta x^4}{4!} u^{(4)}(x_n) + \mathcal{O}(\Delta x^5) \right )\end{aligned}$$except this time we want to leave $u''(x_n)$ on the right hand side.  Doing the same trick as before we have the following expressions:$$\begin{aligned}    u(x_n): & \quad A + B + C = 0 \\    u'(x_n): & \quad A \Delta x - C \Delta x = 0 \\    u''(x_n): & \quad A \frac{\Delta x^2}{2} + C \frac{\Delta x^2}{2} = 1\end{aligned}$$The second equation implies $A = C$ which combined with the third implies$$A = C = \frac{1}{\Delta x^2}$$Finally the first equation gives$$B = -\frac{2}{\Delta x^2}$$leading to the final expression$$\begin{aligned}    u''(x_n) + R(x_n) &= \frac{u(x_{n+1}) - 2 u(x_n) + u(x_{n-1})}{\Delta x^2} \\&\quad+ \frac{1}{\Delta x^2} \left(\frac{\Delta x^3}{3!} u'''(x_n) + \frac{\Delta x^4}{4!} u^{(4)}(x_n)  - \frac{\Delta x^3}{3!} u'''(x_n) + \frac{\Delta x^4}{4!} u^{(4)}(x_n)  \right) + \mathcal{O}(\Delta x^5)\end{aligned}$$with$$R(x_n) = \frac{\Delta x^2}{12} u^{(4)}(x_n) + \mathcal{O}(\Delta x^3)$$<jupyter_code>f = lambda x: numpy.sin(x)
f_dubl_prime = lambda x: -numpy.sin(x)

# Compute derivative
f_dubl_prime_hat = numpy.empty(x_hat.shape)
f_dubl_prime_hat[1:-1] = (f(x_hat[2:]) -2.0 * f(x_hat[1:-1]) + f(x_hat[:-2])) / (delta_x**2)

# Use first-order differences for points at edge of domain
f_dubl_prime_hat[0] = (2.0 * f(x_hat[0]) - 5.0 * f(x_hat[1]) + 4.0 * f(x_hat[2]) - f(x_hat[3])) / delta_x**2
f_dubl_prime_hat[-1] = (2.0 * f(x_hat[-1]) - 5.0 * f(x_hat[-2]) + 4.0 * f(x_hat[-3]) - f(x_hat[-4])) / delta_x**2

fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)

axes.plot(x, f_dubl_prime(x), 'k')
axes.plot(x_hat, f_dubl_prime_hat, 'ro')
axes.set_xlim((x[0], x[-1]))
axes.set_ylim((-1.1, 1.1))
axes.set_xlabel("x")
axes.set_ylabel(r"$f''(x)$")

plt.show()<jupyter_output><empty_output><jupyter_text>General DerivationFor a general finite difference approximation located at $\bar{x}$ to the $k$th derivative with the arbitrary stencil $N \geq k + 1$ points $x_1, \ldots, x_N$ we can use some generalizations of the above method.  Note that although it is common that $\bar{x}$ is one of the stencil points this is not necessary.  We also assume that $u(x)$ is sufficiently smooth so that our Taylor series are valid.At each stencil point we have the approximation$$    u(x_i) = u(\bar{x}) + (x_i - \bar{x})u'(\bar{x}) + \cdots + \frac{1}{k!}(x_i - \bar{x})^k u^{(k)}(\bar{x}) + \cdots.$$Following our methodology above we want to find the linear combination of these Taylor series expansions such that$$    u^{(k)}(\bar{x}) + \mathcal{O}(\Delta x^p) = a_1 u(x_1) + a_2 u(x_2) + a_3 u(x_3) + \cdots + a_n u(x_n).$$Note that $\Delta x$ can vary in general and the asymptotic behavior of the method will be characterized by some sort of average distance or sometimes the maximum distance between the stencil points.Generalizing the approach above with the method of undetermined coefficients we want to eliminate the pieces of the above approximation that are in front of the derivatives less than order $k$.  The condition for this is$$    \frac{1}{(i - 1)!} \sum^N_{j=1} a_j (x_j - \bar{x})^{(i-1)} = \left \{ \begin{aligned}        1 & & \text{if} \quad i - 1 = k, \\        0 & & \text{otherwise}    \end{aligned} \right .$$for $i=1, \ldots, N$.  Assuming the $x_j$ are distinct we can write the system of equations in a Vandermonde system which will have a unique solution.<jupyter_code>import scipy.special

def finite_difference(k, x_bar, x):
    """Compute the finite difference stencil for the kth derivative"""
    
    N = x.shape[0]
    A = numpy.ones((N, N))
    x_row = x - x_bar
    for i in range(1, N):
        A[i, :] = x_row ** i / scipy.special.factorial(i)
    b = numpy.zeros(N)
    b[k] = 1.0
    c = numpy.linalg.solve(A, b)
    return c

print(finite_difference(2, 0.0, numpy.asarray([-1.0, 0.0, 1.0])))
print(finite_difference(1, 0.0, numpy.asarray([-1.0, 0.0, 1.0])))
print(finite_difference(1, -2.0, numpy.asarray([-2.0, -1.0, 0.0, 1.0, 2.0])))
print(finite_difference(2, 0.0, numpy.asarray([-1.0, 0.0, 1.0, 2.0, 3.0, 4.0])) * 12)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17245288920165766, 'fraction_numerical': 0.024786926264758778, 'mean_word_length': 3.633876811594203, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 17, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '274175', 'license': 'MIT', 'n_tokens_mistral': 9423, 'n_tokens_neox': 8961, 'n_words': 3310, 'path': 'Lectures/Lecture 23/.ipynb_checkpoints/Lecture23_FiniteDifference-checkpoint.ipynb', 'repo_name': 'astroarshn2000/PHYS305S20', 'size': '25578'}