docs: add pde section in mdbook (#110)

* docs: add spm battery model example

* add heat equation example

Cargo.toml

@@ -29,7 +29,7 @@ nalgebra = "0.33.2"
 nalgebra-sparse = { version = "0.10", features = ["io"] }
 num-traits = "0.2.17"
 serde = { version = "1.0.215", features = ["derive"] }
-diffsl = { package = "diffsl", version = "0.2.2", optional = true }
+diffsl = { package = "diffsl", version = "0.2.4", optional = true }
 petgraph = "0.6.4"
 faer = "0.19.4"
 suitesparse_sys = { version = "0.1.3", optional = true }

book/src/SUMMARY.md

@@ -9,6 +9,9 @@
     - [Example: Bouncing Ball](./primer/bouncing_ball.md)
   - [DAEs via the Mass Matrix](./primer/the_mass_matrix.md)
     - [Example: Electrical Circuits](./primer/electrical_circuits.md)
+  - [PDEs](./primer/pdes.md)
+    - [Example: Heat Equation](./primer/heat_equation.md)
+    - [Example: Physics-based Battery Simulation](./primer/physics_based_battery_simulation.md)
 - [DiffSol APIs for specifying problems](./specify/specifying_the_problem.md)
   - [ODE equations](./specify/ode_equations.md)
   - [Mass matrix](./specify/mass_matrix.md)

book/src/primer/first_order_odes.md

@@ -6,9 +6,9 @@ Ordinary Differential Equations (ODEs) are sometimes called rate equations becau
 \frac{dc}{dt} = k c
 \\]
 
-where \\(k\\) is a constant that describes the growth rate of the cells. This is a first order ODE, because it involves only the first derivative of the state variable \\(c\\) with respect to time. We can write down a general form for a first order ODE as:
+where \\(k\\) is a constant that describes the growth rate of the cells. This is a first order ODE, because it involves only the first derivative of the state variable \\(c\\) with respect to time.
 
-We can also use the same form to solve a coupled set of equations, say we had *two* populations of cells in the same dish that grow with different rates. We could define the state of the system as the concentrations of the two cell populations, and assign these states to variables \\(c_1\\) and \\(c_2\\). We could then write down both equations as:
+We can also solve a coupled set of equations, say we had *two* populations of cells in the same dish that grow with different rates. We could define the state of the system as the concentrations of the two cell populations, and assign these states to variables \\(c_1\\) and \\(c_2\\). We could then write down both equations as:
 
 \\[
 \begin{align*}

book/src/primer/heat_equation.md

@@ -0,0 +1,228 @@
+# Example: Heat equation
+
+Lets consider a simple example, the heat equation. The heat equation is a PDE that describes how the temperature of a material changes over time. In one dimension, the heat equation is
+
+\\[
+\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2}
+\\]
+
+where \\(u(x, t)\\) is the temperature of the material at position \\(x\\) and time \\(t\\), and \\(D\\) is the thermal diffusivity of the material. To solve this equation, we need to discretize it in space and time. We can use a finite difference method to do this.
+
+## Finite difference method
+
+The finite difference method is a numerical method for discretising a spatial derivative like \\(\frac{\partial^2 u}{\partial x^2}\\). It approximates this *continuous* term by a *discrete* term, in this case the multiplication of a matrix by a vector. We can use this discretisation method to convert the heat equation into a system of ODEs suitable for DiffSol.
+
+We will not go into the details of the finite difference method here but mearly derive a single finite difference approximation for the term \\(\frac{\partial^2 u}{\partial x^2}\\), or \\(u_{xx}\\) using more compact notation. 
+
+The central FD approximation of \\(u_{xx}\\) is:
+
+\\[
+u_{xx} \approx \frac{u(x + h) - 2u(x) + u(x-h)}{h^2}
+\\]
+
+where \\(h\\) is the spacing between points along the x-axis.
+
+We will discretise \\(u_{xx} = 0\\) at \\(N\\) regular points along \\(x\\) from 0 to 1, given by \\(x_1\\), \\(x_2\\), ...
+
+              +----+----+----------+----+> x
+              0   x_1  x_2    ... x_N   1
+
+Using this set of point and the discretised equation, this gives a set of \\(N\\) equations at each interior point on the domain:
+
+\\[
+\frac{v_{i+1} - 2v_i + v_{i-1}}{h^2} \text{ for } i = 1...N
+\\]
+
+where \\(v_i \approx u(x_i)\\)
+
+We will need additional equations at \\(x=0\\) and \\(x=1\\), known as the *boundary conditions*. For this example we will use \\(u(x) = g(x)\\) at \\(x=0\\) and \\(x=1\\) (also known as a non-homogenous Dirichlet bc), so that \\(v_0 = g(0)\\), and \\(v\_{N+1} = g(1)\\), and the equation at \\(x_1\\) becomes:
+
+\\[
+\frac{v_{i+1} - 2v_i + g(0)}{h^2}
+\\]
+
+and the equation at \\(x_N\\) becomes:
+
+\\[
+\frac{g(1) - 2v_i + v_{i-1}}{h^2}
+\\]
+
+We can therefore represent the final \\(N\\) equations in matrix form like so:
+
+\\[
+\frac{1}{h^2}
+\begin{bmatrix} -2      & 1      &         &   &     \\\\
+ 1      & -2     & 1       &       & \\\\
+&\ddots & \ddots  &  \ddots &\\\\
+&        & 1      &  -2     &  1     \\\\
+&        &        &   1     & -2     \end{bmatrix}
+\begin{bmatrix} v_1    \\\\
+v_2    \\\\
+\vdots \\\\
+v_{N-1}\\\\
+v_{N}
+\end{bmatrix} + \frac{1}{h^2} \begin{bmatrix} g(0)    \\\\
+0    \\\\
+\vdots \\\\
+0    \\\\
+g(1)
+\end{bmatrix}
+\\]
+
+The relevant sparse matrix here is \\(A\\), given by
+
+\\[
+A = \begin{bmatrix} -2      & 1      &         &   &     \\\\
+ 1      & -2     & 1       &       & \\\\
+&\ddots & \ddots  &  \ddots &\\\\
+&        & 1      &  -2     &  1     \\\\
+&        &        &   1     & -2     \end{bmatrix}
+\\]
+
+As you can see, the number of non-zero elements grows linearly with the size \\(N\\), so a sparse matrix format is much preferred over a dense matrix holding all \\(N^2\\) elements!
+The additional vector that encodes the boundary conditions is \\(b\\), given by
+
+\\[
+b = \begin{bmatrix} g(0)    \\\\
+0    \\\\
+\vdots \\\\
+0    \\\\
+g(1)
+\end{bmatrix}
+\\]
+
+
+## Method of Lines Approximation
+
+We can use our FD approximation of the spatial derivative to convert the heat equation into a system of ODEs. We can write the heat equation as:
+
+\\[
+\frac{du}{dt} = D \frac{d^2 u}{dx^2} \approx \frac{D}{h^2} (A u + b)
+\\]
+
+where \\(u\\) is a vector of temperatures at each point in space, \\(A\\) and \\(b\\) is the sparse matrix and vector we derived above. This is a system of ODEs that we can solve using DiffSol.
+
+## DiffSol Implementation
+
+```rust
+use diffsol::{
+    DiffSl, OdeBuilder, CraneliftModule, SparseColMat, FaerSparseLU, 
+    OdeSolverMethod
+};
+use plotly::{
+    layout::{Axis, Layout}, Plot, Surface
+};
+# use std::fs;
+# fn main() {
+type M = SparseColMat<f64>;
+type LS = FaerSparseLU<f64>;
+type CG = CraneliftModule;
+
+let eqn = DiffSl::<M, CG>::compile("
+    D { 1.0 }
+    h { 1.0 }
+    g { 0.0 }
+    m { 1.0 }
+    A_ij {
+        (0, 0): -2.0,
+        (0, 1): 1.0,
+        (1, 0): 1.0,
+        (1, 1): -2.0,
+        (1, 2): 1.0,
+        (2, 1): 1.0,
+        (2, 2): -2.0,
+        (2, 3): 1.0,
+        (3, 2): 1.0,
+        (3, 3): -2.0,
+        (3, 4): 1.0,
+        (4, 3): 1.0,
+        (4, 4): -2.0,
+        (4, 5): 1.0,
+        (5, 4): 1.0,
+        (5, 5): -2.0,
+        (5, 6): 1.0,
+        (6, 5): 1.0,
+        (6, 6): -2.0,
+        (6, 7): 1.0,
+        (7, 6): 1.0,
+        (7, 7): -2.0,
+        (7, 8): 1.0,
+        (8, 7): 1.0,
+        (8, 8): -2.0,
+        (8, 9): 1.0,
+        (9, 8): 1.0,
+        (9, 9): -2.0,
+        (9, 10): 1.0,
+        (10, 9): 1.0,
+        (10, 10): -2.0,
+        (10, 11): 1.0,
+        (11, 10): 1.0,
+        (11, 11): -2.0,
+        (11, 12): 1.0,
+        (12, 11): 1.0,
+        (12, 12): -2.0,
+        (12, 13): 1.0,
+        (13, 12): 1.0,
+        (13, 13): -2.0,
+        (13, 14): 1.0,
+        (14, 13): 1.0,
+        (14, 14): -2.0,
+        (14, 15): 1.0,
+        (15, 14): 1.0,
+        (15, 15): -2.0,
+        (15, 16): 1.0,
+        (16, 15): 1.0,
+        (16, 16): -2.0,
+        (16, 17): 1.0,
+        (17, 16): 1.0,
+        (17, 17): -2.0,
+        (17, 18): 1.0,
+        (18, 17): 1.0,
+        (18, 18): -2.0,
+        (18, 19): 1.0,
+        (19, 18): 1.0,
+        (19, 19): -2.0,
+        (19, 20): 1.0,
+        (20, 19): 1.0,
+        (20, 20): -2.0,
+    }
+    b_i { 
+        (0): g,
+        (1:20): 0.0,
+        (20): g,
+    }
+    u_i {
+        (0:5): g,
+        (5:15): g + m,
+        (15:21): g,
+    }
+    heat_i { A_ij * u_j }
+    F_i {
+        D * (heat_i + b_i) / (h * h)
+    }
+").unwrap();
+
+
+let problem = OdeBuilder::<M>::new()
+    .build_from_eqn(eqn)
+    .unwrap();
+let times = (0..50).map(|i| i as f64).collect::<Vec<f64>>();
+let mut solver = problem.bdf::<LS>().unwrap();
+let sol = solver.solve_dense(&times).unwrap();
+
+let x = (0..21).map(|i| i as f64).collect::<Vec<f64>>();
+let y = times;
+let z = sol.col_iter().map(|v| v.iter().copied().collect::<Vec<f64>>()).collect::<Vec<Vec<f64>>>();
+let trace = Surface::new(z).x(x).y(y);
+let mut plot = Plot::new();
+plot.add_trace(trace);
+let layout = Layout::new()
+    .x_axis(Axis::new().title("x"))
+    .y_axis(Axis::new().title("t"))
+    .z_axis(Axis::new().title("u"));
+plot.set_layout(layout);
+let plot_html = plot.to_inline_html(Some("heat-equation"));
+# fs::write("../src/primer/images/heat-equation.html", plot_html).expect("Unable to write file");
+# }
+```
+{{#include images/heat-equation.html}}

book/src/primer/pdes.md

@@ -0,0 +1,21 @@
+# Partial Differential Equations (PDEs)
+
+DiffSol is an ODE solver, but it can also solve PDEs. The idea is to discretize the PDE in space and time, and then solve the resulting system of ODEs. This is called the method of lines.
+
+Discretizing a PDE is a large topic, and there are many ways to do it. Common methods include finite difference, finite volume, finite element, and spectral methods. Finite difference methods are the simplest to understand and implement, so some of the examples in this book will demonstrate this method to give you a flavour of how to solve PDEs with DiffSol. However, in general we recommend that you use another package to discretise you PDE, and then import the resulting ODE system into DiffSol for solving.
+
+## Some useful packages
+
+There are many packages in the Python and Julia ecosystems that can help you discretise your PDE. Here are a few that we know of, but there are many more out there:
+
+Python
+- [FEniCS](https://fenicsproject.org/): A finite element package. Uses the Unified Form Language (UFL) to specify PDEs.
+- [FireDrake](https://firedrakeproject.org/): A finite element package, uses the same UFL as FEniCS.
+- [FiPy](https://www.ctcms.nist.gov/fipy/): A finite volume package.
+- [scikit-fdiff](https://scikit-fdiff.readthedocs.io/en/latest/): A finite difference package.
+
+Julia:
+- [MethodOfLines.jl](https://github.com/SciML/MethodOfLines.jl): A finite difference package.
+- [Gridap.jl](https://github.com/gridap/Gridap.jl): A finite element package.
+
+