From 9f29e24cd68de4250a4bcd0cb88ceade68f5316a Mon Sep 17 00:00:00 2001 From: martinjrobins Date: Tue, 17 Dec 2024 09:46:32 +0000 Subject: [PATCH] docs: misc book fixes --- book/src/benchmarks/python.md | 26 +- book/src/choosing_a_solver.md | 2 +- .../compartmental_models_of_drug_delivery.md | 4 +- book/src/primer/discrete_events.md | 4 +- book/src/primer/electrical_circuits.md | 8 +- book/src/primer/first_order_odes.md | 12 +- book/src/primer/heat_equation.md | 17 +- book/src/primer/images/Mass_spring_damper.svg | 162 +++++++ book/src/primer/images/pk2.svg | 402 ++++++++++++++++++ book/src/primer/modelling_with_odes.md | 13 +- book/src/primer/pdes.md | 4 +- .../physics_based_battery_simulation.md | 6 +- book/src/primer/spring_mass_systems.md | 2 +- book/src/primer/the_mass_matrix.md | 4 +- book/src/solving_the_problem.md | 4 +- book/src/specify/ode_equations.md | 4 +- book/src/specify/specifying_the_problem.md | 16 +- 17 files changed, 630 insertions(+), 60 deletions(-) create mode 100644 book/src/primer/images/Mass_spring_damper.svg create mode 100644 book/src/primer/images/pk2.svg diff --git a/book/src/benchmarks/python.md b/book/src/benchmarks/python.md index 250ec76..b1c8722 100644 --- a/book/src/benchmarks/python.md +++ b/book/src/benchmarks/python.md @@ -1,38 +1,38 @@ # Python (Diffrax & Casadi) -[Diffrax](https://docs.kidger.site/diffrax/) is a Python library for solving ODEs and SDEs implemented using JAX. [Casadi](https://web.casadi.org/) is a C++ library with Python and MATLAB bindings, for solving ODEs and DAEs, nonlinear optimisation and algorithmic differentiation. In this benchmark we compare the performance of the DiffSol implementation with the Diffrax and Casadi libraries. +[Diffrax](https://docs.kidger.site/diffrax/) is a Python library for solving ODEs and SDEs implemented using JAX. [Casadi](https://web.casadi.org/) is a C++ library with Python and MATLAB bindings for solving ODE and nonlinear optimisation problems. In this benchmark we compare the performance of the DiffSol implementation with the Diffrax and Casadi libraries. -As well as demonstrating the performance of the DiffSol library, this benchmark also serves as an example of how to wrap and use DiffSol in other languages. The code for this benchmark can be found [here](https://github.com/martinjrobins/diffsol_python_benchmark). The [maturin](https://www.maturin.rs/) library was used to generate a template for the Python bindings and the CI/CD pipeline neccessary to build the bindings, run pytest tests and build the wheels ready for distribution on PyPI. The [pyo3](https://github.com/PyO3/pyo3) library was used to wrap the DiffSol library in Python. +As well as demonstrating the performance of the DiffSol library, this benchmark also serves as an example of how to wrap and use DiffSol in other languages. The code for this benchmark can be found [here](https://github.com/martinjrobins/diffsol_python_benchmark). The [maturin](https://www.maturin.rs/) library was used to generate a template for the Python bindings and the CI/CD pipeline neccessary to build the wheels ready for distribution on PyPI. The [pyo3](https://github.com/PyO3/pyo3) library was used to wrap the DiffSol library in Python. ## Problem setup -We will use the `robertson_ode` problem for this benchmark. This is a stiff ODE system with 3 equations and 3 unknowns, and is a common benchmark problem for ODE solvers. To illustrate the performance over a range of problem sizes we duplicated the equations by a factor of `ngroups`, so the number of equations is `3 * ngroups`. +We will use the `robertson_ode` problem for this benchmark. This is a stiff ODE system with 3 equations and 3 unknowns, and is a common benchmark problem. To illustrate the performance over a range of problem sizes, the Robertson equations were duplicated a factor of `ngroups`, so the total number of equations solved is `3 * ngroups`. -For the Diffrax implementation we based this on the [example](https://docs.kidger.site/diffrax/examples/stiff_ode/) from the Diffrax documentation, extending this to include the `ngroups` parameter. As with the example, we used the `Kvaerno5` method for the solver. You can see the final implementation of the model [here](https://github.com/martinjrobins/diffsol_python_benchmark/blob/main/diffsol_python_benchmark/diffrax_models.py). +The Diffrax implementation was based this on the [example](https://docs.kidger.site/diffrax/examples/stiff_ode/) in the Diffrax documentation, which was further extending to include the `ngroups` parameter. As is already used in the example, the `Kvaerno5` method was used for the solver. You can see the final implementation of the model [here](https://github.com/martinjrobins/diffsol_python_benchmark/blob/main/diffsol_python_benchmark/diffrax_models.py). -For the Casadi implementation we wrote this from scratch using the libraries Python API. You can see the final implementation of the model [here](https://github.com/martinjrobins/diffsol_python_benchmark/blob/main/diffsol_python_benchmark/casadi_models.py). +The Casadi implementation was written from scratch using Casadi's python API. You can see the final implementation of the model [here](https://github.com/martinjrobins/diffsol_python_benchmark/blob/main/diffsol_python_benchmark/casadi_models.py). -The DiffSol implementation of the model was done using the DiffSL language, and you can see the final implementation of the model [here](https://github.com/martinjrobins/diffsol_python_benchmark/blob/main/diffsol_python_benchmark/diffsol_models.py). +The DiffSol implementation of the model written using the DiffSL language, you can see the final implementation of the model [here](https://github.com/martinjrobins/diffsol_python_benchmark/blob/main/diffsol_python_benchmark/diffsol_models.py). -The final implementation of the benchmark using these models is done [here](https://github.com/martinjrobins/diffsol_python_benchmark/blob/main/bench/bench.py). The DiffSol benchmark is done using the `bdf` solver. For `ngroup` < 20 it uses the `nalgebra` dense matrix and LU solver, and for `ngroups` >= 20 the `faer` sparse matrix and LU solver. +The full implementation of the benchmark presented below can be seen [here](https://github.com/martinjrobins/diffsol_python_benchmark/blob/main/bench/bench.py). The DiffSol benchmark is performed using the `bdf` solver. For `ngroup` < 20 it uses the `nalgebra` dense matrix and LU solver, and for `ngroups` >= 20 the `faer` sparse matrix and LU solver are used. ## Differences between implementations -There are a number of differences between the Diffrax, Casadi and DiffSol implementations that may affect the performance of the solvers. The main differences are: +There are a few key differences between the Diffrax, Casadi and DiffSol implementations that may affect the performance of the solvers. The main differences are: - The Casadi implementation uses sparse matrices, whereas the DiffSol implementation uses dense matrices for `ngroups` < 20, and sparse matrices for `ngroups` >= 20. This will provide an advantage for DiffSol for smaller problems. -- I'm unsure if the Diffrax implementation uses sparse or dense matrices, but it is most likely dense, as JAX only has experimental support for sparse matrices. This will provide an advantage for DiffSol for larger problems. -- The Diffrax implementation uses the `Kvaerno5` method, which is a 5th order implicit Runge-Kutta method. This is different from the BDF method used by both the Casadi and DiffSol implementations. -- Each library was allowed to use multiple threads according to their default settings. The only part of the DiffSol implementation that takes advantage of multiple threads is the `faer` sparse LU solver and matrix. Both the `nalgebra` LU solver / matrix, and the DiffSL generated code are single-threaded only. Diffrax uses JAX, which takes advantage of multiple threads (CPU only, no GPUs were used in these benchmarks). The Casadi implementation also uses multiple threads, but I'm unsure of the details. +- I'm unsure if the Diffrax implementation uses sparse or dense matrices, but it is most likely dense as JAX only has experimental support for sparse matrices. Treating the Jacobian as dense will be a disadvantage for Diffrax for larger problems as the Jacobian is very sparse. +- The Diffrax implementation uses the `Kvaerno5` method (a 5th order implicit Runge-Kutta method). This is different from the BDF method used by both the Casadi and DiffSol implementations. +- Each library was allowed to use multiple threads according to their default settings. The only part of the DiffSol implementation that takes advantage of multiple threads is the `faer` sparse LU solver and matrix. Both the `nalgebra` LU solver, matrix, and the DiffSL generated code are all single-threaded. Diffrax uses JAX, which takes advantage of multiple threads (CPU only, no GPUs were used in these benchmarks). The Casadi implementation also uses multiple threads. ## Results -The benchmarks were run on a Dell PowerEdge R7525 2U rack server, with dual AMD EPYC 7343 3.2Ghz 16C CPU and 128GB Memory. Each benchmark was run using both a low (1e-8) and high (1e-4) tolerances for both `rtol` and `atol`, and with `ngroup` ranging between 1 - 10,000. The results are presented in the following graphs, where the x-axis is the size of the problem `ngroup` and the y-axis is the time taken to solve the problem relative to the time taken by the DiffSol implementation (so `10^0` is the same time as DiffSol, `10^1` is 10 times slower etc.) +The benchmarks were run on a Dell PowerEdge R7525 2U rack server, with dual AMD EPYC 7343 3.2Ghz 16C CPU and 128GB Memory. Each benchmark was run using both a low (1e-8) and high (1e-4) tolerances for both `rtol` and `atol`, and with `ngroup` ranging between 1 - 10,000. The results are presented in the following graph, where the x-axis is the size of the problem `ngroup` and the y-axis is the time taken to solve the problem relative to the time taken by the DiffSol implementation (so `10^0` is the same time as DiffSol, `10^1` is 10 times slower etc.). ![Python](./images/python_plot.svg) DiffSol is faster than both the Casadi and Diffrax implementations over the range of problem sizes and tolerances tested, although the Casadi and DiffSol implementations converge to be similar for larger problems (`ngroups` > 100). -The region that DiffSol really outperforms the other implementations is for smaller problems (`ngroups` < 5), at `ngroups` = 1, Casadi and Diffrax are between 3 - 40 times slower than DiffSol. This small size region are where the dense matrix and solver used is more appropriate for the problem, and the overhead of the other libraries is more significant. The Casadi library needs to traverse a graph of operations to calculate each rhs or jacobian evaluation, whereas the DiffSL JIT compiler will compile to native code using the LLVM backend, along with low-level optimisations that are not available to Casadi. Diffrax as well is significantly slower than DiffSol for smaller problems, and this might be due to (a) Diffrax being a ML library and not optimised for solving stiff ODEs, and (b) double precision is used, which again is not a common use case for ML libraries. +The DiffSol implementation outperforms the other implementations significantly for small problem sizes (`ngroups` < 5). E.g. at `ngroups` = 1, Casadi and Diffrax are between 3 - 40 times slower than DiffSol. At these small problem sizes, the dense matrix and solver used by DiffSol provide an advantage over the sparse solver used by Casadi. Casadi also has additional overhead to evaluate each function evaluation, as it needs to traverse a graph of operations to calculate each rhs or jacobian evaluation, whereas the DiffSL JIT compiler will compile to native code using the LLVM backend, along with low-level optimisations that are not available to Casadi. Diffrax is also significantly slower than DiffSol for smaller problems, this might be due to (a) Diffrax being a ML library and not optimised for solving stiff ODEs, or (b) double precision is used, which again is not a common use case for ML libraries, or (c) perhaps the different solver methods (Kvaerno5 vs BDF) are causing the difference. As the problem sizes get larger, the performance of Diffrax and Casadi improve rapidly relative to DiffSol, but after `ngroups` > 10 the performance of Diffrax drops off again, probably due to JAX not taking advantage of the sparse structure of the problem. The performance of Casadi continues to improve, and for `ngroups` > 100 it is comparable to DiffSol. By the time `ngroups` = 10,000, the performance of Casadi is identical to DiffSol. \ No newline at end of file diff --git a/book/src/choosing_a_solver.md b/book/src/choosing_a_solver.md index e040c83..c7e9815 100644 --- a/book/src/choosing_a_solver.md +++ b/book/src/choosing_a_solver.md @@ -63,7 +63,7 @@ Each solver's state struct implements the [`OdeSolverState`](https://docs.rs/dif For example, say that you wish to bypass the initialisation of the state as you already have the algebraic constraints and so don't need to solve for them. You can use the `new_without_initialise` method on the `OdeSolverState` trait to create a new state without initialising it. You can then use the `as_mut` method to get a mutable reference to the state and set the values manually. -Note that each state struct has a [`as_ref`](https://docs.rs/diffsol/latest/diffsol/ode_solver/state/trait.OdeSolverState.html#tymethod.as_ref) and [`as_mut`](https://docs.rs/diffsol/latest/diffsol/ode_solver/state/trait.OdeSolverState.html#tymethod.as_mut) methods that return a [`StateRef`](https://docs.rs/diffsol/latest/diffsol/ode_solver/state/struct.StateRef.html) or ['StateRefMut`](https://docs.rs/diffsol/latest/diffsol/ode_solver/state/struct.StateRefMut.html) struct respectively. These structs provide a solver-independent way to access the state values so you can use the same code with different solvers. +Note that each state struct has a [`as_ref`](https://docs.rs/diffsol/latest/diffsol/ode_solver/state/trait.OdeSolverState.html#tymethod.as_ref) and [`as_mut`](https://docs.rs/diffsol/latest/diffsol/ode_solver/state/trait.OdeSolverState.html#tymethod.as_mut) methods that return a [`StateRef`](https://docs.rs/diffsol/latest/diffsol/ode_solver/state/struct.StateRef.html) or [`StateRefMut`](https://docs.rs/diffsol/latest/diffsol/ode_solver/state/struct.StateRefMut.html) struct respectively. These structs provide a solver-independent way to access the state values so you can use the same code with different solvers. ```rust # use diffsol::OdeBuilder; diff --git a/book/src/primer/compartmental_models_of_drug_delivery.md b/book/src/primer/compartmental_models_of_drug_delivery.md index 0507a4c..b1687db 100644 --- a/book/src/primer/compartmental_models_of_drug_delivery.md +++ b/book/src/primer/compartmental_models_of_drug_delivery.md @@ -24,7 +24,7 @@ These are often referred to as ADME, and taken together describe the drug concen The body itself is modelled as one or more *compartments*, each of which is defined as a kinetically homogeneous unit (these compartments do not relate to specific organs in the body, unlike Physiologically based pharmacokinetic, PBPK, modeling). There is typically a main *central* compartment into which the drug is administered and from which the drug is excreted from the body, combined with zero or more *peripheral* compartments to which the drug can be distributed to/from the central compartment (See Fig 2). Each peripheral compartment is only connected to the central compartment. -![Fig 2](https://sabs-r3.github.io/software-engineering-projects/fig/pk2.svg) +![Fig 2](images/pk2.svg) The following example PK model describes the two-compartment model shown diagrammatically in the figure above. The time-dependent variables to be solved are the drug quantity in the central and peripheral compartments, $q_c$ and $q_{p1}$ (units: [ng]) respectively. @@ -74,7 +74,7 @@ For the dose function, we will specify a dose of 1000 ng at regular intervals of V_c = 1000 \text{ mL}, \quad V_{p1} = 1000 \text{ mL}, \quad CL = 100 \text{ mL/h}, \quad Q_{p1} = 50 \text{ mL/h} \\] -Let's now solve this system of ODEs using DiffSol. +Let's now solve this system of ODEs using DiffSol. To implement the discrete dose events, we set a stop time for the simulation at each dose event using the [OdeSolverMethod::set_stop_time](https://docs.rs/diffsol/latest/diffsol/ode_solver/method/trait.OdeSolverMethod.html#tymethod.set_stop_time) method. During timestepping we can check the return value of the [OdeSolverMethod::step](https://docs.rs/diffsol/latest/diffsol/ode_solver/method/trait.OdeSolverMethod.html#tymethod.step) method to see if the solver has reached the stop time. If it has, we can apply the dose and continue the simulation. ```rust # fn main() { diff --git a/book/src/primer/discrete_events.md b/book/src/primer/discrete_events.md index f10d2e1..78e2025 100644 --- a/book/src/primer/discrete_events.md +++ b/book/src/primer/discrete_events.md @@ -1,8 +1,8 @@ # Discrete Events -ODEs describe the continuous evolution of a system over time, but many systems also involve discrete events that occur at specific times. For example, in a compartmental model of drug delivery, the administration of a drug is a discrete event that occurs at a specific time. In a bouncing ball model, the collision of the ball with the ground is a discrete event that changes the state of the system. It is normally difficult to model these events using ODEs alone, as they require a different approach to handle the discontinuities in the system. While we can represent discrete events mathematically using delta functions, many ODE solvers are not designed to handle these discontinuities, and may produce inaccurate results or fail to converge during the integration. +ODEs describe the continuous evolution of a system over time, but many systems also involve discrete events that occur at specific times. For example, in a compartmental model of drug delivery, the administration of a drug is a discrete event that occurs at a specific time. In a bouncing ball model, the collision of the ball with the ground is a discrete event that changes the state of the system. It is normally difficult to model these events using ODEs alone, as they require a different approach to handle the discontinuities in the system. While we can represent discrete events mathematically using delta functions, many ODE solvers are not designed to handle discontinuities, and may produce inaccurate results or fail to converge during the integration. -DiffSol provides a way to model discrete events in a system of ODEs by maintaining an internal state of each solver that can be updated when a discrete event occurs. Each solver has an internal state that holds information such as the current time \\(t\\), the current state of the system \\(\mathbf{y}\\), and other solver-specific information. When a discrete event occurs, the user can update the internal state of the solver to reflect the change in the system, and then continue the integration of the ODEs as normal. +DiffSol provides a way to model discrete events in a system of ODEs by allowing users to manipulate the internal state of each solver during the time-stepping. Each solver has an internal state that holds information such as the current time \\(t\\), the current state of the system \\(\mathbf{y}\\), and other solver-specific information. When a discrete event occurs, the user can update the internal state of the solver to reflect the change in the system, and then continue the integration of the ODE as normal. DiffSol also provides a way to stop the integration of the ODEs, either at a specific time or when a specific condition is met. This can be useful for modelling systems with discrete events, as it allows the user to control the integration of the ODEs and to handle the events in a flexible way. diff --git a/book/src/primer/electrical_circuits.md b/book/src/primer/electrical_circuits.md index 7cf4e8d..7cf4fd9 100644 --- a/book/src/primer/electrical_circuits.md +++ b/book/src/primer/electrical_circuits.md @@ -10,7 +10,7 @@ V_s = R | +-------+-------+ ``` -The circuit consists of a resistor \\(R\\), an inductor \\(L\\), and a capacitor \\(C\\) connected in series to a voltage source \\(V_s\\). The voltage across the resistor \\(V\\) is given by Ohm's law: +The circuit consists of a resistor \\(R\\), an inductor \\(L\\), and a capacitor \\(C\\) connected to a voltage source \\(V_s\\). The voltage across the resistor \\(V\\) is given by Ohm's law: \\[ V = R i_R @@ -35,7 +35,7 @@ where \\(i_C\\) is the current flowing through the capacitor. The sum of the cur i_L = i_R + i_C \\] -Thus we have a system of two differential equations and two algebraic equation that describe the evolution of the currents through the resistor, inductor, and capacitor, and the voltage across the resistor. We can write these equations in the general form: +Thus we have a system of two differential equations and two algebraic equation that describe the evolution of the currents through the resistor, inductor, and capacitor; and the voltage across the resistor. We can write these equations in the general form: \\[ M \frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y}, t) @@ -74,10 +74,10 @@ The initial conditions for the system are: The voltage source \\(V_s\\) acts as a forcing function for the system, and we can specify this as sinusoidal function of time. \\[ - V_s(t) = V_0 \sin(omega t) + V_s(t) = V_0 \sin(\omega t) \\] -where \\(omega\\) is the angular frequency of the source. Since this is a low-pass filter, we will choose a high frequency for the source, say \\(omega = 200\\), to demonstrate the filtering effect of the circuit. +where \\(\omega\\) is the angular frequency of the source. Since this is a low-pass filter, we will choose a high frequency for the source, say \\(\omega = 200\\), to demonstrate the filtering effect of the circuit. We can solve this system of equations using DiffSol and plot the current and voltage across the resistor as a function of time. diff --git a/book/src/primer/first_order_odes.md b/book/src/primer/first_order_odes.md index 953be84..2881375 100644 --- a/book/src/primer/first_order_odes.md +++ b/book/src/primer/first_order_odes.md @@ -1,6 +1,6 @@ # Explicit First Order ODEs -Ordinary Differential Equations (ODEs) are sometimes called rate equations because they describe how the *rate of change* of a system depends on its current *state*. For example, lets assume we wish to model the growth of a population of cells within a petri dish. We could define the *state* of the system as the concentration of cells in the dish, and assign this state to a variable \\(c\\). The *rate of change* of the system would then be the rate at which the concentration of cells changes with time, which we could denote as \\(\frac{dc}{dt}\\). We know that our cells will grow at a rate proportional to the current concentration of cells, so we could write this as: +Ordinary Differential Equations (ODEs) are often called rate equations because they describe how the *rate of change* of a system depends on its current *state*. For example, lets assume we wish to model the growth of a population of cells within a petri dish. We could define the *state* of the system as the concentration of cells in the dish, and assign this state to a variable \\(c\\). The *rate of change* of the system would then be the rate at which the concentration of cells changes with time, which we could denote as \\(\frac{dc}{dt}\\). We know that our cells will grow at a rate proportional to the current concentration of cells, so this can be written as: \\[ \frac{dc}{dt} = k c @@ -8,7 +8,7 @@ Ordinary Differential Equations (ODEs) are sometimes called rate equations becau 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 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 extend this further to solve multiple equations simultaineously, in order to model the rate of change of more than one quantity. For example, 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\\). could then write down both equations as: \\[ \begin{align*} @@ -17,7 +17,7 @@ We can also solve a coupled set of equations, say we had *two* populations of ce \end{align*} \\] -then combine them in a vector form as: +and then combine them in a vector form as: \\[ \begin{bmatrix} @@ -29,8 +29,7 @@ k_2 c_2 \end{bmatrix} \\] - -By defining a new *vector* of state variables \\(\mathbf{y} = [c_1, c_2]\\) and a a function \\(\mathbf{f}(\mathbf{y}, t) = \begin{bmatrix} k_1 c_1 \\\\ k_2 c_2 \end{bmatrix}\\), we are left with the standard form of a *explicit* first order ODE system: +By defining a new *vector* of state variables \\(\mathbf{y} = [c_1, c_2]\\) and a vector valued function \\(\mathbf{f}(\mathbf{y}, t) = \begin{bmatrix} k_1 c_1 \\\\ k_2 c_2 \end{bmatrix}\\), we are left with the standard form of a *explicit* first order ODE system: \\[ \frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y}, t) @@ -44,5 +43,4 @@ We need one more piece of information to solve this system of ODEs: the initial \mathbf{y}(0) = \begin{bmatrix} 10 \\\\ 5 \end{bmatrix} \\] - -Many ODE solver libraries, like DiffSol, require users to provide their ODEs in the form of a set of explicit first order ODEs. Given both the system of ODEs and the initial conditions, the solver can then integrate the equations forward in time to find the solution \\(\mathbf{y}(t)\\). This is the general process for solving ODEs, so it is important to know how to translate your problem into a set of first order ODEs, and thus to the general form of a explicit first order ODE system shown above. In the next two sections, we will look at two examples of first order ODE systems: population dynamics and infectious disease, and solve them using DiffSol. \ No newline at end of file +Many ODE solver libraries, like DiffSol, require users to provide their ODEs in the form of a set of explicit first order ODEs. Given both the system of ODEs and the initial conditions, the solver can then integrate the equations forward in time to find the solution \\(\mathbf{y}(t)\\). This is the general process for solving ODEs, so it is important to know how to translate your problem into a set of first order ODEs, and thus to the general form of a explicit first order ODE system shown above. In the next two sections, we will look at an example of a first order ODE system in the area of population dynamics, and then solve it using DiffSol. \ No newline at end of file diff --git a/book/src/primer/heat_equation.md b/book/src/primer/heat_equation.md index 5845290..c2d66f6 100644 --- a/book/src/primer/heat_equation.md +++ b/book/src/primer/heat_equation.md @@ -6,7 +6,7 @@ Lets consider a simple example, the heat equation. The heat equation is a PDE th \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. +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 discretise the spatial derivative, and then solve the resulting system of ODEs using DiffSol. ## Finite difference method @@ -27,7 +27,7 @@ We will discretise \\(u_{xx} = 0\\) at \\(N\\) regular points along \\(x\\) from +----+----+----------+----+> 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: +Using this set of points and the discrete approximation, 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 @@ -94,13 +94,20 @@ g(1) ## 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: +We can use our FD approximation of the spatial derivative to convert the heat equation into a system of ODEs. Starting from our original definition of the heat equation: \\[ -\frac{du}{dt} = D \frac{d^2 u}{dx^2} \approx \frac{D}{h^2} (A u + b) +\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} +\\] + +and using our finite difference approximation and definition of the sparse matrix \\(A\\) and vector \\(b\\), this becomes: + + +\\[ +\frac{du}{dt} = \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. +where \\(u\\) is a vector of temperatures at each point in space. This is a system of ODEs that we can solve using DiffSol. ## DiffSol Implementation diff --git a/book/src/primer/images/Mass_spring_damper.svg b/book/src/primer/images/Mass_spring_damper.svg new file mode 100644 index 0000000..30db6c8 --- /dev/null +++ b/book/src/primer/images/Mass_spring_damper.svg @@ -0,0 +1,162 @@ + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + m + k + x + + + + + + c + diff --git a/book/src/primer/images/pk2.svg b/book/src/primer/images/pk2.svg new file mode 100644 index 0000000..bee5ba4 --- /dev/null +++ b/book/src/primer/images/pk2.svg @@ -0,0 +1,402 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + Central Compartment + + Peripheral Compartment + + Dose + + Clearance + + + + + + + + + + + + + + + + + + + + + diff --git a/book/src/primer/modelling_with_odes.md b/book/src/primer/modelling_with_odes.md index 2207f23..540120b 100644 --- a/book/src/primer/modelling_with_odes.md +++ b/book/src/primer/modelling_with_odes.md @@ -1,15 +1,18 @@ # Modelling with ODEs -Ordinary Differential Equations (ODEs) are a powerful tool for modelling a wide range of physical systems. Unlike purely data-driven models, ODEs are based on the underlying physics, biology, or chemistry of the system being modelled. This makes them particularly useful for predicting the behaviour of a system under conditions that have not been observed before. In this section, we will introduce the basics of ODE modelling, and illustrate their use with a series of examples written using the DiffSol crate. +Ordinary Differential Equations (ODEs) are a powerful tool for modelling a wide range of physical systems. Unlike purely data-driven models, ODEs are based on the underlying physics, biology, or chemistry of the system being modelled. This makes them particularly useful for predicting the behaviour of a system under conditions that have not been observed. In this section, we will introduce the basics of ODE modelling, and illustrate their use with a series of examples written using the DiffSol crate. The topics covered in this section are: - [First Order ODEs](./primer/first_order_odes.md): First order ODEs are the simplest type of ODE. Any ODE system can be written as a set of first order ODEs, so libraries like DiffSol are designed such that the user provides their equations in this form. - [Example: Population Dynamics](./primer/population_dynamics.md): A simple example of a first order ODE system, modelling the interaction of predator and prey populations. -- [Higher Order ODEs](./primer/higher_order_odes.md): Higher order ODEs are ODEs that involve derivatives of order higher than one. These can be converted to a system of first order ODEs, which is the form that DiffSol expects. +- [Higher Order ODEs](./primer/higher_order_odes.md): Higher order ODEs are equations that involve derivatives of order greater than one. These can be converted to a system of first order ODEs, which is the form that DiffSol expects. - [Example: Spring-mass systems](./primer/spring_mass_systems.md): A simple example of a higher order ODE system, modelling the motion of a damped spring-mass system. -- [Discrete Events](./primer/discrete_events.md): Discrete events are events that occur at specific times or when the system is in a particular state, rather than continuously. These can be modelled using ODEs by treating the events as changes in the system's state. DiffSol provides an API to detect and handle these events. - - [Example: Compartmental models of Drug Delivery](./primer/compartmental_models_of_drug_delivery.md): Pharmacokinetic models are a common example of systems with discrete events, where the drug is absorbed, distributed, metabolised, and excreted by the body. The drug is often administered at discrete times, and the model must account for this. - - [Example: Bouncing Ball](./primer/bouncing_ball.md): A simple example of a system where the discrete event occurs when the ball hits the ground, instead of a specific time. +- [Discrete Events](./primer/discrete_events.md): Discrete events are events that occur at specific times or when the system is in a particular state, rather than continuously. These can be modelled by treating the events as changes in the ODE system's state. DiffSol provides an API to detect and handle these events. + - [Example: Compartmental models of Drug Delivery](./primer/compartmental_models_of_drug_delivery.md): Pharmacokinetic models describe how a drug is absorbed, distributed, metabolised, and excreted by the body. They are a common example of systems with discrete events, as the drug is often administered at discrete times. + - [Example: Bouncing Ball](./primer/bouncing_ball.md): A simple example of a system where the discrete event occurs when the ball hits the ground, instead of at a specific time. - [DAEs via the Mass Matrix](./primer/the_mass_matrix.md): Differential Algebraic Equations (DAEs) are a generalisation of ODEs that include algebraic equations as well as differential equations. DiffSol can solve DAEs by treating them as ODEs with a mass matrix. This section explains how to use the mass matrix to solve DAEs. - [Example: Electrical Circuits](./primer/electrical_circuits.md): Electrical circuits are a common example of DAEs, here we will model a simple low-pass LRC filter circuit. +- [PDEs](./primer/pdes.md): Partial Differential Equations (PDEs) are a generalisation of ODEs that involve derivatives with respect to more than one variable (e.g. a spatial variable). DiffSol can be used to solver PDEs using the method of lines, where the spatial derivatives are discretised to form a system of ODEs. + - [Example: Heat Equation](./primer/heat_equation.md): The heat equation describes how heat diffuses in a domain over time. We will solve the heat equation in a 1D domain with Dirichlet boundary conditions. + - [Example: Physics-based Battery Simulation](./primer/physics_based_battery_simulation.md): A more complex example of a PDE system, modelling the charge and discharge of a lithium-ion battery. For this example we will use the PyBaMM library to form the ODE system, and DiffSol to solve it. \ No newline at end of file diff --git a/book/src/primer/pdes.md b/book/src/primer/pdes.md index f809535..b079586 100644 --- a/book/src/primer/pdes.md +++ b/book/src/primer/pdes.md @@ -2,11 +2,11 @@ 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. +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 your 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: +There are many packages in the Python and Julia ecosystems that can help you discretise your PDE. Here are a few, 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. diff --git a/book/src/primer/physics_based_battery_simulation.md b/book/src/primer/physics_based_battery_simulation.md index ed0a7cf..292e101 100644 --- a/book/src/primer/physics_based_battery_simulation.md +++ b/book/src/primer/physics_based_battery_simulation.md @@ -69,9 +69,11 @@ V - V_{\text{min}} = 0, ## Solving the Single Particle Model using DiffSol -Rather than write down the equations ourselves, we will rely on the [PyBaMM library](https://pybamm.org/) to generate the equations for us. PyBaMM is a Python library that can generate a wide variety of physics-based battery models, using different parameterisations, physics and operating conditions. Combined with [a tool](https://github.com/martinjrobins/pybamm2diffsl) that takes a PyBaMM model and writes it out in the DiffSL language, we can generate [a DiffSL file](src/spm.ds) that can be used to solve the equations of the SPM model described above. We can then use the DiffSol crate to solve the model and calculate the terminal voltage of the battery over a range of current rates. +The equations above describe the Single Particle Model of a lithium-ion battery, but they are relativly complex and difficult to discretise compared with the simple heat equation PDE that we saw in the [Heat Equation](./heat_equation.md) section. -The code below reads in the DiffSL file, compiles it, and then solves the equation for different current rates. We wish to stop the simulation when either the final time is reached, or when one of the stopping conditions is met. We will output the terminal voltage of the battery at regular intervals during the simulation, because the terminal voltage can change more rapidly than the state variables $c_n$ and $c_p$, particularly during the "knee" of the discharge curve. +Rather than derive and write down the discretised equations outselves, we will instead rely on the [PyBaMM library](https://pybamm.org/) to generate the equations for us. PyBaMM is a Python library that can generate a wide variety of physics-based battery models, using different parameterisations, physics and operating conditions. Combined with [a tool](https://github.com/martinjrobins/pybamm2diffsl) that takes a PyBaMM model and writes it out in the DiffSL language, we can generate [a DiffSL file](src/spm.ds) that can be used to solve the equations of the SPM model described above. We can then use the DiffSol crate to solve the model and calculate the terminal voltage of the battery over a range of current rates. + +The code below reads in the DiffSL file, compiles it, and then solves the equation for different current rates. We wish to stop the simulation when either the final time is reached, or when one of the stopping conditions is met. We will output the terminal voltage of the battery at regular intervals during the simulation, because the terminal voltage can change more rapidly than the state variables \\(c_n\\) and \\(c_p\\), particularly during the "knee" of the discharge curve. The discretised equations result in sparse matrices, so we use the sparse matrix and linear solver modules provided by the [faer](https://github.com/sarah-quinones/faer-rs) crate to solve the equations efficiently. diff --git a/book/src/primer/spring_mass_systems.md b/book/src/primer/spring_mass_systems.md index 3321acc..6df22b4 100644 --- a/book/src/primer/spring_mass_systems.md +++ b/book/src/primer/spring_mass_systems.md @@ -2,7 +2,7 @@ We will model a [damped spring-mass system](https://en.wikipedia.org/wiki/Mass-spring-damper_model) using a second order ODE. The system consists of a mass \\(m\\) attached to a spring with spring constant \\(k\\), and a damping force proportional to the velocity of the mass with damping coefficient \\(c\\). -![Spring-mass system](https://upload.wikimedia.org/wikipedia/commons/thumb/4/45/Mass_spring_damper.svg/330px-Mass_spring_damper.svg.png) + The equation of motion for the mass can be written as: diff --git a/book/src/primer/the_mass_matrix.md b/book/src/primer/the_mass_matrix.md index bc77a38..db4ab0b 100644 --- a/book/src/primer/the_mass_matrix.md +++ b/book/src/primer/the_mass_matrix.md @@ -12,8 +12,8 @@ where \\(\mathbf{y}\\) is the vector of state variables, \\(\mathbf{y}'\\) is th \\[ \begin{align*} -\frac{d\mathbf{y}}{dt} &= \mathbf{f}(\mathbf{y}, \mathbf{y}', t) \\\\ -0 = \mathbf{g}(\mathbf{y}, t) +\frac{d\mathbf{y}}{dt} &= \mathbf{f}(\mathbf{y}, t) \\\\ +0 &= \mathbf{g}(\mathbf{y}, t) \end{align*} \\] diff --git a/book/src/solving_the_problem.md b/book/src/solving_the_problem.md index 3e18bb4..720e22c 100644 --- a/book/src/solving_the_problem.md +++ b/book/src/solving_the_problem.md @@ -7,8 +7,8 @@ Each solver implements the [`OdeSolverMethod`](https://docs.rs/diffsol/latest/di DiffSol has a few high-level solution functions on the `OdeSolverMethod` trait that are the easiest way to solve your equations: - [`solve`](https://docs.rs/diffsol/latest/diffsol/ode_solver/method/trait.OdeSolverMethod.html#method.solve) - solve the problem from an initial state up to a specified time, returning the solution at all the internal timesteps used by the solver. - [`solve_dense`](https://docs.rs/diffsol/latest/diffsol/ode_solver/method/trait.OdeSolverMethod.html#method.solve_dense) - solve the problem from an initial state, returning the solution at a `Vec` of times provided by the user. -- ['solve_dense_sensitivities`](https://docs.rs/diffsol/latest/diffsol/ode_solver/method/trait.OdeSolverMethod.html#method.solve_dense_sensitivities) - solve the forward sensitivity problem from an initial state, returning the solution at a `Vec` of times provided by the user. -- ['solve_adjoint'](https://docs.rs/diffsol/latest/diffsol/ode_solver/method/trait.OdeSolverMethod.html#method.solve_adjoint) - solve the adjoint sensitivity problem from an initial state to a final time, returning the integration of the output function over time as well as its gradient with respect to the initial state. +- [`solve_dense_sensitivities`](https://docs.rs/diffsol/latest/diffsol/ode_solver/method/trait.OdeSolverMethod.html#method.solve_dense_sensitivities) - solve the forward sensitivity problem from an initial state, returning the solution at a `Vec` of times provided by the user. +- [`solve_adjoint`](https://docs.rs/diffsol/latest/diffsol/ode_solver/method/trait.OdeSolverMethod.html#method.solve_adjoint) - solve the adjoint sensitivity problem from an initial state to a final time, returning the integration of the output function over time as well as its gradient with respect to the initial state. The following example shows how to solve a simple ODE problem using the `solve` method on the `OdeSolverMethod` trait. diff --git a/book/src/specify/ode_equations.md b/book/src/specify/ode_equations.md index 8e23cee..75e14ed 100644 --- a/book/src/specify/ode_equations.md +++ b/book/src/specify/ode_equations.md @@ -1,8 +1,6 @@ # ODE equations -The simplest way to create a new ode problem is to use the [`OdeBuilder`](https://docs.rs/diffsol/latest/diffsol/ode_solver/builder/struct.OdeBuilder.html) struct. -You can use methods to set the equations to be solve, initial time, initial step size, relative & absolute tolerances, and parameters, or leave them at their default values. -Then, call the `build` method to create a new problem. +To create a new ode problem, use the [`OdeBuilder`](https://docs.rs/diffsol/latest/diffsol/ode_solver/builder/struct.OdeBuilder.html) struct. You can create a builder using the `OdeBuilder::new` method, and then chain methods to set the equations to be solve, initial time, initial step size, relative & absolute tolerances, and parameters, or leave them at their default values. Then, call the `build` method to create a new problem. Below is an example of how to create a new ODE problem using the `OdeBuilder` struct. The specific problem we will solve is the logistic equation diff --git a/book/src/specify/specifying_the_problem.md b/book/src/specify/specifying_the_problem.md index 053ddd4..7ddbfae 100644 --- a/book/src/specify/specifying_the_problem.md +++ b/book/src/specify/specifying_the_problem.md @@ -23,15 +23,13 @@ The user can also optionally specify a root function \\(r(t, y, p)\\) that can b ## DiffSol problem APIs -DiffSol has three main APIs for specifying problems: -- The [`OdeBuilder`](https://docs.rs/diffsol/latest/diffsol/ode_solver/builder/struct.OdeBuilder.html) struct, where the user can specify the functions above using Rust closures. - This is the easiest API to use from Rust, and is the recommended API for most users. -- The [`OdeEquations`](https://docs.rs/diffsol/latest/diffsol/ode_solver/equations/trait.OdeEquations.html) trait - where the user can implement the functions above on their own structs. - This API is more flexible than the `OdeBuilder` API, but is more complex to use. It is useful if you have custom data structures and code that you want to use to evaluate - your functions that does not fit within the `OdeBuilder` API. -- The [`DiffSl`](https://docs.rs/diffsol/latest/diffsol/ode_solver/diffsl/struct.DiffSl.html) struct, where the user can specify the functions above using the [DiffSL](https://martinjrobins.github.io/diffsl/) - Domain Specific Language (DSL). This API is behind a feature flag (`diffsl` if you want to use the slower cranelift backend, `diffsl-llvm*` if you want to use the faster LLVM backend), but has the best API if you want to use DiffSL from a higher-level language like Python or R while still having similar performance. +Specifying a problem in DiffSol is done via the [`OdeBuilder`](https://docs.rs/diffsol/latest/diffsol/ode_solver/builder/struct.OdeBuilder.html) struct, using the `OdeBuilder::new` method to create a new builder, and then chaining methods to set the equations to be solved, initial time, initial step size, relative & absolute tolerances, and parameters, or leaving them at their default values. Then, call the `build` method to create a new problem. + +Users can specify the equations to be solved via three main APIs, ranging from the simplest to the most complex (but also the most flexible): +- The [`DiffSl`](https://docs.rs/diffsol/latest/diffsol/ode_solver/diffsl/struct.DiffSl.html) struct allows users to specify the equations above using the [DiffSL](https://martinjrobins.github.io/diffsl/) + Domain Specific Language (DSL). This API is behind a feature flag (`diffsl` if you want to use the slower cranelift backend, `diffsl-llvm*` if you want to use the faster LLVM backend). This is the easiest API to use as it can use automatic differentiation to calculate the neccessary gradients, and is the best API if you want to use DiffSL from a higher-level language like Python or R while still having similar performance to Rust. +- The [`OdeBuilder`](https://docs.rs/diffsol/latest/diffsol/ode_solver/builder/struct.OdeBuilder.html) struct also has methods to set the equations using rust closures (see e.g. [OdeBuilder::rhs_implicit](https://docs.rs/diffsol/latest/diffsol/ode_solver/builder/struct.OdeBuilder.html#method.rhs_implicit)). This API is convenient if you want to stick to pure rust code without using DiffSL and the JIT compiler, but requires you to calculate the gradients of the equations yourself. +- Implementing the [`OdeEquations`](https://docs.rs/diffsol/latest/diffsol/ode_solver/equations/trait.OdeEquations.html) trait allows users to implement the equations on their own structs. This API is the most flexible as it allows users to use custom data structures and code that might not fit within the `OdeBuilder` API. However, it is the most verbose API and requires users to be more familiar with the various DiffSol traits.