I encourage you to play around with the juptyer notebook for this lecture - you can copy the code with the this notebook button on the side of this page.
Warning
These notes are mainly a record of what we discussed and are not a substitute for attending the lectures and reading books! If anything is unclear/wrong, let me know and I will update the notes.
We will consider the following abstract formulation of a boundary value problem. Suppose we have a differential operator \mathcal L, a “nice” domain \Omega \subset \mathbb R^d, and dataf (source or sink) and v (boundary values). Seek u : \Omega \to \mathbb R such that
\begin{align}
\mathcal L[u] &= f &&\text{in } \Omega \nonumber\\
u &= v &&\text{on } \partial \Omega .
\tag{BVP}
\end{align}
This is a problem with Dirichlet boundary conditions (the values of the solution are prescribed on the boundary).
Example 12.1 (Linear Advection-Diffusion-Reaction Equations) Consider
\begin{align}
\mathcal L[u] := -\Delta u + p \cdot \nabla u + q u
\end{align}
In the context of semiconductor physics, u(x) represents the concentration of charge carriers at a position x \in \Omega. Each of these terms relates to the physics in question:
-\Delta : diffusion; electrons move from areas of high concentration to low,
p \cdot \nabla u : advection (drift); the solution is driven with velocity p (due to e.g. an external electric field),
q u : reaction; electrons and electron holes are being created or destroyed.
The same equation also models e.g pollutant transport in a river, drug delivery in a blood vessel, asset price in finance, heat transfer in fluids, ……
Here, we solve this problem numerically in one spatial dimension:
functionadr_1d(n, p, q, f) h =1/n x = h:h:(1-h) # Finite Difference Coefficients# -u'' ≈ (-u_{i-1} + 2u_i - u_{i+1}) / h^2# pu' ≈ p * (u_{i+1} - u_{i-1}) / (2h) [Central Difference] d = (2/h^2+ q) *ones(n-1) d₋ = (-1/h^2- p/(2h)) *ones(n-2) d₊ = (-1/h^2+ p/(2h)) *ones(n-2) A =Tridiagonal(d₋, d, d₊) b =f.(x) u = A \ b# Dirichlet boundary conditions u = [0.0, u..., 0.0] x = [0.0, x..., 1.0]return x, uend# u'' = diffusion# p = advection# q = reaction# f = source f(x) =1.0# diffusionp, q =0., 0.x, u =adr_1d(100, p, q, f)P1 =plot(x, u, title=L"Diffusion: $-u′′ + %$p u′ + %$q u = 1$", xlabel=L"Position $x$", ylabel=L"Concentration $u$", label=L"u(x)", lw=2, m=1)display(P1)# Advection (drift)p, q =50., 0.x, u =adr_1d(100, p, q, f)P2 =plot(x, u, title=L"Advection: $-u′′ + %$p u′ + %$q u = 1$", xlabel=L"Position $x$", ylabel=L"Concentration $u$", label=L"u(x)", lw=2, m=1)display(P2)# Reaction (decay)p, q =0., 100.x, u =adr_1d(100, p, q, f)P3 =plot(x, u, title=L"Reaction: $-u′′ + %$p u′ + %$q u = 1$", xlabel=L"Position $x$", ylabel=L"Concentration $u$", label=L"u(x)", lw=2, m=1)display(P3)# Advection + Reactionp, q =100., 100.x, u =adr_1d(100, p, q, f)P4 =plot(x, u, title=L"$-u′′ + %$p u′ + %$q u = 1$", xlabel=L"Position $x$", ylabel=L"Concentration $u$", label=L"u(x)", lw=2, m=1)display(P4)
by solving a linear system of equations (more generally, p and q can also be functions of x).
In order to do this, we needed to approximate the linear differential operatorsu \mapsto u' and u \mapsto u'' with matrices. Suppose u : [0,1] \to \mathbb R solves the equation above with the boundary conditions u(0) = u(1) = 0. We try to approximate the solution of this equation on the grid \{x_j\}_{j=0}^n where x_j = \frac{j}{n}. We will let h:= \frac{1}{n} be the mesh size and let u_j \in \mathbb R^{n+1} be an approximation to u(x_j). Since u(0) = u(1) = 0, we can simply set u_0 = u_n = 0 and we want to write an equation for the remaining variables u_1,\dots,u_{n-1}.
(these are backwards and central differences, resp.)
Exercise 12.2 Show that \partial^0 = \frac{1}{2}( \partial^- + \partial^+).
Definition 12.1 We say a finite difference is consistent if the error E(h) between the exact derivative and the finite difference goes to zero as h \to 0. Moreover, we say the finite difference has order of accuracy p if E(h)=O(h^p) as h \to 0.
Exercise 12.3 What is the order of accuracy of \partial^+, \partial^-, \partial^0?
Exercise 12.4 Given u_1,\dots,u_{n-1} (recall that we are considering u_0 = u_n = 0), one can approximate the derivative of u at x_j by computing \partial^- u_j, \partial^+ u_j, or \partial^0 u_j (there are of course other choices!). Write these operations as matrices so that u'(x_j) \approx (D^- u)_j, (D^+ u)_j, and (D^0 u)_j.
Hint, these differentiation matrices are (n-1)\times(n-1).
