diff --git a/Cargo.toml b/Cargo.toml index 5b6f18ab..5282c1f9 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -11,43 +11,38 @@ repository = "https://github.com/martinjrobins/diffsol" # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html [features] -default = ["nalgebra", "faer"] +default = ["nalgebra", "faer", "diffsl"] faer = [] nalgebra = [] sundials = ["suitesparse_sys", "bindgen", "cc"] suitesparse = ["suitesparse_sys"] -diffsl-cranelift = ["diffsl-no-llvm", "diffsl"] -diffsl = [ ] +diffsl = ["dep:diffsl"] diffsl-llvm = [] -diffsl-llvm13 = ["diffsl13-0", "diffsl-llvm", "diffsl"] -diffsl-llvm14 = ["diffsl14-0", "diffsl-llvm", "diffsl"] -diffsl-llvm15 = ["diffsl15-0", "diffsl-llvm", "diffsl"] -diffsl-llvm16 = ["diffsl16-0", "diffsl-llvm", "diffsl"] -diffsl-llvm17 = ["diffsl17-0", "diffsl-llvm", "diffsl"] +diffsl-llvm13 = ["diffsl/llvm13-0", "diffsl", "diffsl-llvm"] +diffsl-llvm14 = ["diffsl/llvm14-0", "diffsl", "diffsl-llvm"] +diffsl-llvm15 = ["diffsl/llvm15-0", "diffsl", "diffsl-llvm"] +diffsl-llvm16 = ["diffsl/llvm16-0", "diffsl", "diffsl-llvm"] +diffsl-llvm17 = ["diffsl/llvm17-0", "diffsl", "diffsl-llvm"] [dependencies] -nalgebra = "0.33" +nalgebra = "0.33.2" nalgebra-sparse = { version = "0.10", features = ["io"] } num-traits = "0.2.17" -serde = { version = "1.0.196", features = ["derive"] } -diffsl-no-llvm = { package = "diffsl", version = "=0.2.0", optional = true } -diffsl13-0 = { package = "diffsl", version = "=0.2.0", features = ["llvm13-0"], optional = true } -diffsl14-0 = { package = "diffsl", version = "=0.2.0", features = ["llvm14-0"], optional = true } -diffsl15-0 = { package = "diffsl", version = "=0.2.0", features = ["llvm15-0"], optional = true } -diffsl16-0 = { package = "diffsl", version = "=0.2.0", features = ["llvm16-0"], optional = true } -diffsl17-0 = { package = "diffsl", version = "=0.2.0", features = ["llvm17-0"], optional = true } +serde = { version = "1.0.215", features = ["derive"] } +diffsl = { package = "diffsl", version = "0.2.2", optional = true } petgraph = "0.6.4" faer = "0.19.4" suitesparse_sys = { version = "0.1.3", optional = true } -thiserror = "1.0.63" +thiserror = "2.0.3" [dev-dependencies] -insta = { version = "1.34.0", features = ["yaml"] } +insta = { version = "1.41.1", features = ["yaml"] } criterion = { version = "0.5.1" } +skeptic = "0.13.7" [build-dependencies] bindgen = { version = "0.70.1", optional = true } -cc = { version = "1.0.99", optional = true } +cc = { version = "1.2.0", optional = true } [[bench]] name = "ode_solvers" diff --git a/book/Cargo.toml b/book/Cargo.toml new file mode 100644 index 00000000..0252ac9d --- /dev/null +++ b/book/Cargo.toml @@ -0,0 +1,10 @@ +[package] +name = "book" +version = "0.1.0" +edition = "2021" + +[dependencies] +diffsol = { path = "..", features = ["diffsl"] } +nalgebra = "0.33.2" +faer = "0.19.4" +plotly = "0.10.0" diff --git a/book/book.toml b/book/book.toml index 3ee2c540..5755e02d 100644 --- a/book/book.toml +++ b/book/book.toml @@ -10,8 +10,8 @@ mathjax-support = true [preprocessor.keeper] command = "mdbook-keeper" -manifest_dir = ".." -externs = ["diffsol", "nalgebra", "faer"] +manifest_dir = "." +externs = ["diffsol", "nalgebra", "faer", "plotly"] # is_workspace = true -build_features = ["diffsl-cranelift"] +# build_features = ["plotly"] diff --git a/book/src/SUMMARY.md b/book/src/SUMMARY.md index b554b9c4..e5f6cf53 100644 --- a/book/src/SUMMARY.md +++ b/book/src/SUMMARY.md @@ -1,17 +1,26 @@ # Summary - -- [Specifying the problem](./specifying_the_problem.md) - - [ODE equations](./ode_equations.md) - - [Mass matrix](./mass_matrix.md) - - [Root finding](./root_finding.md) - - [Forward Sensitivity](./forward_sensitivity.md) - - [Custom Problem Structs](./custom_problem_structs.md) - - [Non-linear functions](./non_linear_functions.md) - - [Constant functions](./constant_functions.md) - - [Linear functions](./linear_functions.md) - - [Putting it all together](./putting_it_all_together.md) - - [DiffSL](./diffsl.md) - - [Sparse problems](./sparse_problems.md) +- [Modelling with ODEs](./primer/modelling_with_odes.md) + - [Explicit First Order ODEs](./primer/first_order_odes.md) + - [Example: Population Dynamics](./primer/population_dynamics.md) + - [Higher Order ODEs](./primer/higher_order_odes.md) + - [Example: Spring-mass systems](./primer/spring_mass_systems.md) + - [Discrete Events](./primer/discrete_events.md) + - [Example: Compartmental models of Drug Delivery](./primer/compartmental_models_of_drug_delivery.md) + - [Example: Bouncing Ball](./primer/bouncing_ball.md) + - [DAEs via the Mass Matrix](./primer/the_mass_matrix.md) + - [Example: Electrical Circuits](./primer/electrical_circuits.md) +- [DiffSol APIs for specifying problems](./specify/specifying_the_problem.md) + - [ODE equations](./specify/ode_equations.md) + - [Mass matrix](./specify/mass_matrix.md) + - [Root finding](./specify/root_finding.md) + - [Forward Sensitivity](./specify/forward_sensitivity.md) + - [Custom Problem Structs](./specify/custom/custom_problem_structs.md) + - [Non-linear functions](./specify/custom/non_linear_functions.md) + - [Constant functions](./specify/custom/constant_functions.md) + - [Linear functions](./specify/custom/linear_functions.md) + - [Putting it all together](./specify/custom/putting_it_all_together.md) + - [DiffSL](./specify/diffsl.md) + - [Sparse problems](./specify/sparse_problems.md) - [Choosing a solver](./choosing_a_solver.md) - [Initialisation](./initialisation.md) - [Solving the problem](./solving_the_problem.md) diff --git a/book/src/diffsl.md b/book/src/diffsl.md index d37a0216..5e7097fb 100644 --- a/book/src/diffsl.md +++ b/book/src/diffsl.md @@ -1,59 +1 @@ # DiffSL - -Thus far we have used Rust code to specify the problem we want to solve. This is fine if you are using DiffSol from Rust, but what if you want to use DiffSol from a higher-level language like Python or R? -For this usecase we have designed a Domain Specific Language (DSL) called DiffSL that can be used to specify the problem. DiffSL is not a general purpose language but is tightly constrained to -the specification of a system of ordinary differential equations. It features a relativly simple syntax that consists of writing a series of tensors (dense or sparse) that represent the equations of the system. -For more detail on the syntax of DiffSL see the [DiffSL book](https://martinjrobins.github.io/diffsl/). This section will focus on how to use DiffSL to specify a problem in DiffSol. - - -## DiffSL Context - -The main struct that is used to specify a problem in DiffSL is the [`DiffSlContext`](https://docs.rs/diffsol/latest/diffsol/ode_solver/diffsl/struct.DiffSlContext.html) struct. Creating this struct -Just-In-Time (JIT) compiles your DiffSL code into a form that can be executed efficiently by DiffSol. - -```rust -# fn main() { -use diffsol::{DiffSl, CraneliftModule}; -type M = nalgebra::DMatrix; -type CG = CraneliftModule; - -let eqn = DiffSl::::compile(" - in = [r, k] - r { 1.0 } - k { 1.0 } - u { 0.1 } - F { r * u * (1.0 - u / k) } - out { u } -").unwrap(); -# } -``` - -Once you have created the `DiffSlContext` struct you can use it to create a problem using the `build_diffsl` method on the [`OdeBuilder`](https://docs.rs/diffsol/latest/diffsol/ode_solver/builder/struct.OdeBuilder.html) struct. - - -```rust -# fn main() { -# use diffsol::{DiffSl, CraneliftModule}; -use diffsol::{OdeBuilder, Bdf, OdeSolverMethod, OdeSolverState}; -# type M = nalgebra::DMatrix; -# type CG = CraneliftModule; - - -# let eqn = DiffSl::::compile(" -# in = [r, k] -# r { 1.0 } -# k { 1.0 } -# u { 0.1 } -# F { r * u * (1.0 - u / k) } -# out { u } -# ").unwrap(); -let problem = OdeBuilder::new() -.rtol(1e-6) -.p([1.0, 10.0]) -.build_from_eqn(eqn).unwrap(); -let mut solver = Bdf::default(); -let t = 0.4; -let state = OdeSolverState::new(&problem, &solver).unwrap(); -let _soln = solver.solve(&problem, state, t).unwrap(); -# } -``` diff --git a/book/src/lib.rs b/book/src/lib.rs new file mode 100644 index 00000000..e69de29b diff --git a/book/src/primer/bouncing_ball.md b/book/src/primer/bouncing_ball.md new file mode 100644 index 00000000..1c1b242a --- /dev/null +++ b/book/src/primer/bouncing_ball.md @@ -0,0 +1,145 @@ +# Example: Bouncing Ball + +Modelling a bouncing ball is a simple and intuitive example of a system with discrete events. The ball is dropped from a height \\(h\\) and bounces off the ground with a coefficient of restitution \\(e\\). When the ball hits the ground, its velocity is reversed and scaled by the coefficient of restitution, and the ball rises and then continues to fall until it hits the ground again. This process repeats until halted. + +The second order ODE that describes the motion of the ball is given by: + +\\[ +\frac{d^2x}{dt^2} = -g +\\] + +where \\(x\\) is the position of the ball, \\(t\\) is time, and \\(g\\) is the acceleration due to gravity. We can rewrite this as a system of two first order ODEs by introducing a new variable for the velocity of the ball: + +\\[ +\begin{align*} +\frac{dx}{dt} &= v \\\\ +\frac{dv}{dt} &= -g +\end{align*} +\\] + +where \\(v = \frac{dx}{dt}\\) is the velocity of the ball. This is a system of two first order ODEs, which can be written in vector form as: + +\\[ +\frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y}, t) +\\] + +where + +\\[ +\mathbf{y} = \begin{bmatrix} x \\\\ v \end{bmatrix} +\\] + +and + +\\[ +\mathbf{f}(\mathbf{y}, t) = \begin{bmatrix} v \\\\ -g \end{bmatrix} +\\] + +The initial conditions for the ball, including the height from which it is dropped and its initial velocity, are given by: + +\\[ +\mathbf{y}(0) = \begin{bmatrix} h \\\\ 0 \end{bmatrix} +\\] + +When the ball hits the ground, we need to update the velocity of the ball according to the coefficient of restitution, which is the ratio of the velocity after the bounce to the velocity before the bounce. The velocity after the bounce \\(v'\\) is given by: + +\\[ +v' = -e v +\\] + +where \\(e\\) is the coefficient of restitution. However, to implement this in our ODE solver, we need to detect when the ball hits the ground. We can do this by using DiffSol's event handling feature, which allows us to specify a function that is equal to zero when the event occurs, i.e. when the ball hits the ground. This function \\(g(\mathbf{y}, t)\\) is called an event or root function, and for our bouncing ball problem, it is given by: + +\\[ +g(\mathbf{y}, t) = x +\\] + +where \\(x\\) is the position of the ball. When the ball hits the ground, the event function will be zero and DiffSol will stop the integration, and we can update the velocity of the ball accordingly. + +In code, the bouncing ball problem can be solved using DiffSol as follows: + +```rust +# fn main() { +# use std::fs; +use diffsol::{ + DiffSl, CraneliftModule, OdeBuilder, Bdf, OdeSolverState, OdeSolverMethod, + OdeSolverStopReason, +}; +use plotly::{ + Plot, Scatter, common::Mode, layout::Layout, layout::Axis +}; +type M = nalgebra::DMatrix; +type CG = CraneliftModule; + +let eqn = DiffSl::::compile(" + g { 9.81 } h { 10.0 } + u_i { + x = h, + v = 0, + } + F_i { + v, + -g, + } + stop { + x, + } +").unwrap(); + +let e = 0.8; +let problem = OdeBuilder::new().build_from_eqn(eqn).unwrap(); +let mut solver = Bdf::default(); +let state = OdeSolverState::new(&problem, &solver).unwrap(); +solver.set_problem(state, &problem).unwrap(); + +let mut x = Vec::new(); +let mut v = Vec::new(); +let mut t = Vec::new(); +let final_time = 10.0; + +// save the initial state +x.push(solver.state().unwrap().y[0]); +v.push(solver.state().unwrap().y[1]); +t.push(0.0); + +// solve until the final time is reached +solver.set_stop_time(final_time).unwrap(); +loop { + match solver.step() { + Ok(OdeSolverStopReason::InternalTimestep) => (), + Ok(OdeSolverStopReason::RootFound(t)) => { + // get the state when the event occurred + let mut y = solver.interpolate(t).unwrap(); + + // update the velocity of the ball + y[1] *= -e; + + // make sure the ball is above the ground + y[0] = y[0].max(f64::EPSILON); + + // set the state to the updated state + solver.state_mut().unwrap().y.copy_from(&y); + solver.state_mut().unwrap().dy[0] = y[1]; + *solver.state_mut().unwrap().t = t; + }, + Ok(OdeSolverStopReason::TstopReached) => break, + Err(_) => panic!("unexpected solver error"), + } + x.push(solver.state().unwrap().y[0]); + v.push(solver.state().unwrap().y[1]); + t.push(solver.state().unwrap().t); +} +let mut plot = Plot::new(); +let x = Scatter::new(t.clone(), x).mode(Mode::Lines).name("x"); +let v = Scatter::new(t, v).mode(Mode::Lines).name("v"); +plot.add_trace(x); +plot.add_trace(v); + +let layout = Layout::new() + .x_axis(Axis::new().title("t")) + .y_axis(Axis::new()); +plot.set_layout(layout); +let plot_html = plot.to_inline_html(Some("bouncing-ball")); +# fs::write("../src/primer/images/bouncing-ball.html", plot_html).expect("Unable to write file"); +# } +``` +{{#include images/bouncing-ball.html}} \ No newline at end of file diff --git a/book/src/primer/compartmental_models_of_drug_delivery.md b/book/src/primer/compartmental_models_of_drug_delivery.md new file mode 100644 index 00000000..09f992d2 --- /dev/null +++ b/book/src/primer/compartmental_models_of_drug_delivery.md @@ -0,0 +1,153 @@ +# Example: Compartmental models of Drug Delivery + +The field of Pharmacokinetics (PK) provides a quantitative basis for describing the delivery of a drug to a patient, the diffusion of that drug through the plasma/body tissue, and the subsequent clearance of the drug from the patient's system. PK is used to ensure that there is sufficient concentration of the drug to maintain the required efficacy of the drug, while ensuring that the concentration levels remain below the toxic threshold. Pharmacokinetic (PK) models are often combined with Pharmacodynamic (PD) models, which model the positive effects of the drug, such as the binding of a drug to the biological target, and/or undesirable side effects, to form a full PKPD model of the drug-body interaction. This example will only focus on PK, neglecting the interaction with a PD model. + +![Fig 1](https://sabs-r3.github.io/software-engineering-projects/fig/pk1.jpg) + +PK enables the following processes to be quantified: + +- Absorption +- Distribution +- Metabolism +- Excretion + +These are often referred to as ADME, and taken together describe the drug concentration in the body when medicine is prescribed. These ADME processes are typically described by zeroth-order or first-order *rate* reactions modelling the dynamics of the quantity of drug $q$, with a given rate parameter $k$, for example: + +\\[ +\frac{dq}{dt} = -k^*, +\\] + +\\[ +\frac{dq}{dt} = -k q. +\\] + +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) + +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. + +\\[ +\frac{dq_c}{dt} = \text{Dose}(t) - \frac{q_c}{V_c} CL - Q_{p1} \left ( \frac{q_c}{V_c} - \frac{q_{p1}}{V_{p1}} \right ), +\\] + +\\[ +\frac{dq_{p1}}{dt} = Q_{p1} \left ( \frac{q_c}{V_c} - \frac{q_{p1}}{V_{p1}} \right ). +\\] + +This model describes an *intravenous bolus* dosing protocol, with a linear clearance from the central compartment (non-linear clearance processes are also possible, but not considered here). The dose function $\text{Dose}(t)$ will consist of instantaneous doses of $X$ ng of the drug at one or more time points. The other input parameters to the model are: +- \\(V_c\\) [mL], the volume of the central compartment +- \\(V_{p1}\\) [mL], the volume of the first peripheral compartment +- \\(CL\\) [mL/h], the clearance/elimination rate from the central compartment +- \\(Q_{p1}\\) [mL/h], the transition rate between central compartment and peripheral compartment 1 + +We will solve this system of ODEs using the DiffSol crate. Rather than trying to write down the dose function as a mathematical function, we will neglect the dose function from the equations and instead using DiffSol's API to specify the dose at specific time points. + +First lets write down the equations in the standard form of a first order ODE system: + +\\[ +\frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y}, t) +\\] + +where + +\\[ +\mathbf{y} = \begin{bmatrix} q_c \\\\ q_{p1} \end{bmatrix} +\\] + +and + +\\[ +\mathbf{f}(\mathbf{y}, t) = \begin{bmatrix} - \frac{q_c}{V_c} CL - Q_{p1} \left ( \frac{q_c}{V_c} - \frac{q_{p1}}{V_{p1}} \right ) \\\\ Q_{p1} \left ( \frac{q_c}{V_c} - \frac{q_{p1}}{V_{p1}} \right ) \end{bmatrix} +\\] + +We will also need to specify the initial conditions for the system: + +\\[ +\mathbf{y}(0) = \begin{bmatrix} 0 \\\\ 0 \end{bmatrix} +\\] + +For the dose function, we will specify a dose of 1000 ng at regular intervals of 6 hours. We will also specify the other parameters of the model: + +\\[ +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. + +```rust +# fn main() { +# use std::fs; +use diffsol::{ + DiffSl, CraneliftModule, OdeBuilder, Bdf, OdeSolverState, OdeSolverMethod, + OdeSolverStopReason, +}; +use plotly::{ + Plot, Scatter, common::Mode, layout::Layout, layout::Axis +}; +type M = nalgebra::DMatrix; +type CG = CraneliftModule; + +let eqn = DiffSl::::compile(" + Vc { 1000.0 } Vp1 { 1000.0 } CL { 100.0 } Qp1 { 50.0 } + u_i { + qc = 0, + qp1 = 0, + } + F_i { + - qc / Vc * CL - Qp1 * (qc / Vc - qp1 / Vp1), + Qp1 * (qc / Vc - qp1 / Vp1), + } +").unwrap(); + +let problem = OdeBuilder::new().build_from_eqn(eqn).unwrap(); +let mut solver = Bdf::default(); +let doses = vec![(0.0, 1000.0), (6.0, 1000.0), (12.0, 1000.0), (18.0, 1000.0)]; +let state = OdeSolverState::new(&problem, &solver).unwrap(); +solver.set_problem(state, &problem).unwrap(); + +let mut q_c = Vec::new(); +let mut q_p1 = Vec::new(); +let mut time = Vec::new(); + +// apply the first dose and save the initial state +solver.state_mut().unwrap().y[0] = doses[0].1; +q_c.push(solver.state().unwrap().y[0]); +q_p1.push(solver.state().unwrap().y[1]); +time.push(0.0); + +// solve and apply the remaining doses +for (t, dose) in doses.into_iter().skip(1) { + solver.set_stop_time(t).unwrap(); + loop { + let ret = solver.step(); + q_c.push(solver.state().unwrap().y[0]); + q_p1.push(solver.state().unwrap().y[1]); + time.push(solver.state().unwrap().t); + match ret { + Ok(OdeSolverStopReason::InternalTimestep) => continue, + Ok(OdeSolverStopReason::TstopReached) => break, + _ => panic!("unexpected solver error"), + } + } + solver.state_mut().unwrap().y[0] += dose; +} +let mut plot = Plot::new(); +let q_c = Scatter::new(time.clone(), q_c).mode(Mode::Lines).name("q_c"); +let q_p1 = Scatter::new(time, q_p1).mode(Mode::Lines).name("q_p1"); +plot.add_trace(q_c); +plot.add_trace(q_p1); + +let layout = Layout::new() + .x_axis(Axis::new().title("t [h]")) + .y_axis(Axis::new().title("amount [ng]")); +plot.set_layout(layout); +let plot_html = plot.to_inline_html(Some("drug-delivery")); +# fs::write("../src/primer/images/drug-delivery.html", plot_html).expect("Unable to write file"); +# } +``` +{{#include images/drug-delivery.html}} + + + diff --git a/book/src/primer/discrete_events.md b/book/src/primer/discrete_events.md new file mode 100644 index 00000000..f10d2e13 --- /dev/null +++ b/book/src/primer/discrete_events.md @@ -0,0 +1,9 @@ +# 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. + +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 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. + +The [Solving the Problem](../solving_the_problem.md) and [Root Finding](../specify/root_finding.md) sections provides an introduction to the API for solving ODEs and detecting events with DiffSol. In the next two sections, we will look at two examples of systems with discrete events: compartmental models of drug delivery and bouncing ball models, and solve them using DiffSol using this API. \ No newline at end of file diff --git a/book/src/primer/electrical_circuits.md b/book/src/primer/electrical_circuits.md new file mode 100644 index 00000000..ec3162d6 --- /dev/null +++ b/book/src/primer/electrical_circuits.md @@ -0,0 +1,150 @@ +# Example: Electrical Circuits + +Lets consider the following low-pass LRC filter circuit: + +```plaintext + +---L---+---C---+ + | | | +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: + +\\[ +V = R i_R +\\] + +where \\(i_R\\) is the current flowing through the resistor. The voltage across the inductor is given by: + +\\[ +\frac{di_L}{dt} = \frac{V_s - V}{L} +\\] + +where \\(di_L/dt\\) is the rate of change of current with respect to time. The voltage across the capacitor is the same as the voltage across the resistor and the equation for an ideal capacitor is: + +\\[ +\frac{dV}{dt} = \frac{i_C}{C} +\\] + +where \\(i_C\\) is the current flowing through the capacitor. The sum of the currents flowing into and out of the top-center node of the circuit must be zero according to Kirchhoff's current law: + + +\\[ +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: + +\\[ +M \frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y}, t) +\\] + +where + +\\[ +\mathbf{y} = \begin{bmatrix} i_R \\\\ i_L \\\\ i_C \\\\ V \end{bmatrix} +\\] + +and + +\\[ +\mathbf{f}(\mathbf{y}, t) = \begin{bmatrix} V - R i_R \\\\ \frac{V_s - V}{L} \\\\ i_L - i_R - i_C \\\\ \frac{i_C}{C} \end{bmatrix} +\\] + +The mass matrix \\(M\\) has one on the diagonal for the differential equations and zero for the algebraic equation. + +\\[ +M = \begin{bmatrix} 0 & 0 & 0 & 0 \\\\ 0 & 1 & 0 & 0 \\\\ 0 & 0 & 0 & 0 \\\\ 0 & 0 & 0 & 1 \end{bmatrix} +\\] + +Instead of providing the mass matrix explicitly, the DiffSL language specifies the multiplication of the mass matrix with the derivative term, \\(M \frac{d\mathbf{y}}{dt}\\), which is given by: + +\\[ +M \frac{d\mathbf{y}}{dt} = \begin{bmatrix} 0 \\\\ \frac{di_L}{dt} \\\\ 0 \\\\ \frac{dV}{dt} \end{bmatrix} +\\] + +The initial conditions for the system are: + +\\[ +\mathbf{y}(0) = \begin{bmatrix} 0 \\\\ 0 \\\\ 0 \\\\ 0 \end{bmatrix} +\\] + +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) +\\] + +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. + +```rust +# fn main() { +# use std::fs; +use diffsol::{ + DiffSl, CraneliftModule, OdeBuilder, Bdf, OdeSolverState, OdeSolverMethod +}; +use plotly::{ + Plot, Scatter, common::Mode, layout::Layout, layout::Axis +}; +type M = nalgebra::DMatrix; +type CG = CraneliftModule; + +let eqn = DiffSl::::compile(" + R { 100.0 } L { 1.0 } C { 0.001 } V0 { 10 } omega { 100.0 } + Vs { V0 * sin(omega * t) } + u_i { + iR = 0, + iL = 0, + iC = 0, + V = 0, + } + dudt_i { + diRdt = 0, + diLdt = 0, + diCdt = 0, + dVdt = 0, + } + M_i { + 0, + diLdt, + 0, + dVdt, + } + F_i { + V - R * iR, + (Vs - V) / L, + iL - iR - iC, + iC / C, + } + out_i { + iR, + } +").unwrap(); +let problem = OdeBuilder::new() + .build_from_eqn(eqn).unwrap(); +let mut solver = Bdf::default(); +let state = OdeSolverState::new(&problem, &solver).unwrap(); +let (ys, ts) = solver.solve(&problem, state, 1.0).unwrap(); + +let ir: Vec<_> = ys.row(0).into_iter().copied().collect(); +let t: Vec<_> = ts.into_iter().collect(); + +let ir = Scatter::new(t.clone(), ir).mode(Mode::Lines); + +let mut plot = Plot::new(); +plot.add_trace(ir); + +let layout = Layout::new() + .x_axis(Axis::new().title("t")) + .y_axis(Axis::new().title("current")); +plot.set_layout(layout); +let plot_html = plot.to_inline_html(Some("electrical-circuit")); +# fs::write("../src/primer/images/electrical-circuit.html", plot_html).expect("Unable to write file"); +# } +``` +{{#include images/electrical-circuit.html}} + diff --git a/book/src/primer/first_order_odes.md b/book/src/primer/first_order_odes.md new file mode 100644 index 00000000..c6b9c9e3 --- /dev/null +++ b/book/src/primer/first_order_odes.md @@ -0,0 +1,48 @@ +# 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: + +\\[ +\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: + +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: + +\\[ +\begin{align*} +\frac{dc_1}{dt} &= k_1 c_1 \\\\ +\frac{dc_2}{dt} &= k_2 c_2 +\end{align*} +\\] + +then combine them in a vector form as: + +\\[ +\begin{bmatrix} +\frac{dc_1}{dt} \\\\ +\frac{dc_2}{dt} +\end{bmatrix} = \begin{bmatrix} +k_1 c_1 \\\\ +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: + +\\[ +\frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y}, t) +\\] + +This is an explicit equation for the derivative of the state, \\(\frac{d\mathbf{y}}{dt}\\), as a function of only of the state variables \\(\mathbf{y}\\) and of time \\(t\\). + +We need one more piece of information to solve this system of ODEs: the initial conditions for the populations at time \\(t = 0\\). For example, if we started with a concentration of 10 for the first population and 5 for the second population, we would write: + +\\[ +\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 diff --git a/book/src/primer/higher_order_odes.md b/book/src/primer/higher_order_odes.md new file mode 100644 index 00000000..da5a8ea6 --- /dev/null +++ b/book/src/primer/higher_order_odes.md @@ -0,0 +1,38 @@ +# Higher Order ODEs + +The *order* of an ODE is the highest derivative that appears in the equation. So far, we have only looked at first order ODEs, which involve only the first derivative of the state variable with respect to time. However, many physical systems are described by higher order ODEs, which involve second or higher derivatives of the state variable. A simple example of a second order ODE is the motion of a mass under the influence of gravity. The equation of motion for the mass can be written as: + +\\[ +\frac{d^2x}{dt^2} = -g +\\] + +where \\(x\\) is the position of the mass, \\(t\\) is time, and \\(g\\) is the acceleration due to gravity. This is a second order ODE because it involves the second derivative of the position with respect to time. + +Higher order ODEs can always be rewritten as a system of first order ODEs by introducing new variables. For example, we can rewrite the second order ODE above as a system of two first order ODEs by introducing a new variable for the velocity of the mass: + +\\[ +\begin{align*} +\frac{dx}{dt} &= v \\\\ +\frac{dv}{dt} &= -g +\end{align*} +\\] + +where \\(v = \frac{dx}{dt}\\) is the velocity of the mass. This is a system of two first order ODEs, which can be written in vector form as: + +\\[ +\frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y}, t) +\\] + +where + +\\[ +\mathbf{y} = \begin{bmatrix} x \\\\ v \end{bmatrix} +\\] + +and + +\\[ +\mathbf{f}(\mathbf{y}, t) = \begin{bmatrix} v \\\\ -g \end{bmatrix} +\\] + +In the next section, we'll look at another example of a higher order ODE system: the spring-mass system, and solve this using DiffSol. \ No newline at end of file diff --git a/book/src/primer/images/bouncing-ball.html b/book/src/primer/images/bouncing-ball.html new file mode 100644 index 00000000..90c4d620 --- /dev/null +++ b/book/src/primer/images/bouncing-ball.html @@ -0,0 +1,4 @@ +
+ \ No newline at end of file diff --git a/book/src/primer/images/drug-delivery.html b/book/src/primer/images/drug-delivery.html new file mode 100644 index 00000000..1e9d94ff --- /dev/null +++ b/book/src/primer/images/drug-delivery.html @@ -0,0 +1,4 @@ +
+ \ No newline at end of file diff --git a/book/src/primer/images/electrical-circuit.html b/book/src/primer/images/electrical-circuit.html new file mode 100644 index 00000000..e96cc4f4 --- /dev/null +++ b/book/src/primer/images/electrical-circuit.html @@ -0,0 +1,4 @@ +
+ \ No newline at end of file diff --git a/book/src/primer/images/prey-predator.html b/book/src/primer/images/prey-predator.html new file mode 100644 index 00000000..02f8c303 --- /dev/null +++ b/book/src/primer/images/prey-predator.html @@ -0,0 +1,4 @@ +
+ \ No newline at end of file diff --git a/book/src/primer/images/prey-predator2.html b/book/src/primer/images/prey-predator2.html new file mode 100644 index 00000000..54385fa7 --- /dev/null +++ b/book/src/primer/images/prey-predator2.html @@ -0,0 +1,4 @@ +
+ \ No newline at end of file diff --git a/book/src/primer/images/spring-mass-system.html b/book/src/primer/images/spring-mass-system.html new file mode 100644 index 00000000..f64d79c8 --- /dev/null +++ b/book/src/primer/images/spring-mass-system.html @@ -0,0 +1,4 @@ +
+ \ No newline at end of file diff --git a/book/src/primer/modelling_with_odes.md b/book/src/primer/modelling_with_odes.md new file mode 100644 index 00000000..2207f23d --- /dev/null +++ b/book/src/primer/modelling_with_odes.md @@ -0,0 +1,15 @@ +# 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. + +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. + - [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. +- [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. + \ No newline at end of file diff --git a/book/src/primer/population_dynamics.md b/book/src/primer/population_dynamics.md new file mode 100644 index 00000000..3637e02a --- /dev/null +++ b/book/src/primer/population_dynamics.md @@ -0,0 +1,181 @@ +# Population Dynamics - Predator-Prey Model + +In this example, we will model the population dynamics of a predator-prey system using a set of first order ODEs. The [Lotka-Volterra equations](https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations) are a classic example of a predator-prey model, and describe the interactions between two species in an ecosystem. The first species is a prey species, which we will call \\(x\\), and the second species is a predator species, which we will call \\(y\\). + +The rate of change of the prey population is governed by two terms: growth and predation. The growth term represents the natural increase in the prey population in the absence of predators, and is proportional to the current population of prey. The predation term represents the rate at which the predators consume the prey, and is proportional to the product of the prey and predator populations. The rate of change of the prey population can be written as: + +\\[ +\frac{dx}{dt} = a x - b x y +\\] + +where \\(a\\) is the growth rate of the prey population, and \\(b\\) is the predation rate. + +The rate of change of the predator population is also governed by two terms: death and growth. The death term represents the natural decrease in the predator population in the absence of prey, and is proportional to the current population of predators. The growth term represents the rate at which the predators reproduce, and is proportional to the product of the prey and predator populations, since the predators need to consume the prey to reproduce. The rate of change of the predator population can be written as: + +\\[ +\frac{dy}{dt} = c x y - d y +\\] + +where \\(c\\) is the reproduction rate of the predators, and \\(d\\) is the death rate. + +The Lotka-Volterra equations are a simple model of predator-prey dynamics, and make several assumptions that may not hold in real ecosystems. For example, the model assumes that the prey population grows exponentially in the absence of predators, that the predator population decreases linearly in the absence of prey, and that the spatial distribution of the species has no effect. Despite these simplifications, the Lotka-Volterra equations capture some of the essential features of predator-prey interactions, such as oscillations in the populations and the dependence of each species on the other. When modelling with ODEs, it is important to consider the simplest model that captures the behaviour of interest, and to be aware of the assumptions that underlie the model. + +Putting the two equations together, we have a system of two first order ODEs: + +\\[ +\frac{x}{dt} = a x - b x y \\\\ +\frac{y}{dt} = c x y - d y +\\] + +which can be written in vector form as: + +\\[ +\begin{bmatrix} +\frac{dx}{dt} \\\\ +\frac{dy}{dt} +\end{bmatrix} = \begin{bmatrix} +a x - b x y \\\\ +c x y - d y +\end{bmatrix} +\\] + +or in the general form of a first order ODE system: + +\\[ +\frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y}, t) +\\] + +where + +\\[\mathbf{y} = \begin{bmatrix} x \\\\ y \end{bmatrix} \\] + +and + +\\[\mathbf{f}(\mathbf{y}, t) = \begin{bmatrix} a x - b x y \\\\ c x y - d y \end{bmatrix}\\] + +We also have initial conditions for the populations at time \\(t = 0\\). We can set both populations to 1 at the start like so: + +\\[ +\mathbf{y}(0) = \begin{bmatrix} 1 \\\\ 1 \end{bmatrix} +\\] + +Let's solve this system of ODEs using the DiffSol crate. We will use the [DiffSL](https://martinjrobins.github.io/diffsl/) domain-specific language to specify the problem. We could have also used [Rust closures](specify/ode_equations.md), but this allows us to illustrate the modelling process with a minimum of Rust syntax. + +```rust +# fn main() { +# use std::fs; +use diffsol::{ + DiffSl, CraneliftModule, OdeBuilder, Bdf, OdeSolverState, OdeSolverMethod +}; +use plotly::{ + Plot, Scatter, common::Mode, layout::Layout, layout::Axis +}; +type M = nalgebra::DMatrix; +type CG = CraneliftModule; + +let eqn = DiffSl::::compile(" + a { 2.0/3.0 } b { 4.0/3.0 } c { 1.0 } d { 1.0 } + u_i { + y1 = 1, + y2 = 1, + } + F_i { + a * y1 - b * y1 * y2, + c * y1 * y2 - d * y2, + } +").unwrap(); +let problem = OdeBuilder::new() + .build_from_eqn(eqn).unwrap(); +let mut solver = Bdf::default(); +let state = OdeSolverState::new(&problem, &solver).unwrap(); +let (ys, ts) = solver.solve(&problem, state, 40.0).unwrap(); + +let prey: Vec<_> = ys.row(0).into_iter().copied().collect(); +let predator: Vec<_> = ys.row(1).into_iter().copied().collect(); +let time: Vec<_> = ts.into_iter().collect(); + +let prey = Scatter::new(time.clone(), prey).mode(Mode::Lines).name("Prey"); +let predator = Scatter::new(time, predator).mode(Mode::Lines).name("Predator"); + +let mut plot = Plot::new(); +plot.add_trace(prey); +plot.add_trace(predator); + +let layout = Layout::new() + .x_axis(Axis::new().title("t")) + .y_axis(Axis::new().title("population")); +plot.set_layout(layout); +let plot_html = plot.to_inline_html(Some("prey-predator")); +# fs::write("../src/primer/images/prey-predator.html", plot_html).expect("Unable to write file"); +# } +``` +{{#include images/prey-predator.html}} + +A phase plane plot of the predator-prey system is a useful visualisation of the dynamics of the system. This plot shows the prey population on the x-axis and the predator population on the y-axis. Trajectories in the phase plane represent the evolution of the populations over time. Lets reframe the equations to introduce a new parameter \\(y_0\\) which is the initial predator and prey population. We can then plot the phase plane for different values of \\(y_0\\) to see how the system behaves for different initial conditions. + +Our initial conditions are now: + +\\[ +\mathbf{y}(0) = \begin{bmatrix} y_0 \\\\ y_0 \end{bmatrix} +\\] + +so we can solve this system for different values of \\(y_0\\) and plot the phase plane for each case. We will use similar code as above, but we will introduce our new parameter and loop over different values of \\(y_0\\) + +```rust +# fn main() { +# use std::fs; +use diffsol::{ + DiffSl, CraneliftModule, OdeBuilder, Bdf, OdeSolverState, OdeSolverMethod +}; +use plotly::{ + Plot, Scatter, common::Mode, layout::Layout, layout::Axis +}; +type M = nalgebra::DMatrix; +type CG = CraneliftModule; + +let eqn = DiffSl::::compile(" + in = [ y0 ] + y0 { 1.0 } + a { 2.0/3.0 } b { 4.0/3.0 } c { 1.0 } d { 1.0 } + u_i { + y1 = y0, + y2 = y0, + } + F_i { + a * y1 - b * y1 * y2, + c * y1 * y2 - d * y2, + } +").unwrap(); + +let mut problem = OdeBuilder::new().p([1.0]).build_from_eqn(eqn).unwrap(); +let mut solver = Bdf::default(); + +let mut plot = Plot::new(); +for y0 in (1..6).map(f64::from) { + problem.set_params(nalgebra::DVector::from_element(1, y0)).unwrap(); + + let state = OdeSolverState::new(&problem, &solver).unwrap(); + let (ys, _ts) = solver.solve(&problem, state, 40.0).unwrap(); + + let prey: Vec<_> = ys.row(0).into_iter().copied().collect(); + let predator: Vec<_> = ys.row(1).into_iter().copied().collect(); + + let phase = Scatter::new(prey, predator) + .mode(Mode::Lines).name(format!("y0 = {}", y0)); + plot.add_trace(phase); + + // release problem and state to set new parameters in the next iteration + solver.take_state().unwrap(); +} + +let layout = Layout::new() + .x_axis(Axis::new().title("x")) + .y_axis(Axis::new().title("y")); +plot.set_layout(layout); +let plot_html = plot.to_inline_html(Some("prey-predator2")); +# fs::write("../src/primer/images/prey-predator2.html", plot_html).expect("Unable to write file"); +# } +``` +{{#include images/prey-predator2.html}} + + diff --git a/book/src/primer/spatial_population_dynamics.md b/book/src/primer/spatial_population_dynamics.md new file mode 100644 index 00000000..068e6672 --- /dev/null +++ b/book/src/primer/spatial_population_dynamics.md @@ -0,0 +1 @@ +# Example: Spatial Population Dynamics diff --git a/book/src/primer/spring_mass_systems.md b/book/src/primer/spring_mass_systems.md new file mode 100644 index 00000000..1eee2fe3 --- /dev/null +++ b/book/src/primer/spring_mass_systems.md @@ -0,0 +1,73 @@ +# Example: Spring-mass systems + +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: + +\\[ +m \frac{d^2x}{dt^2} = -k x - c \frac{dx}{dt} +\\] + +where \\(x\\) is the position of the mass, \\(t\\) is time, and the negative sign on the right hand side indicates that the spring force and damping force act in the opposite direction to the displacement of the mass. + +We can convert this to a system of two first order ODEs by introducing a new variable for the velocity of the mass: + +\\[ +\begin{align*} +\frac{dx}{dt} &= v \\\\ +\frac{dv}{dt} &= -\frac{k}{m} x - \frac{c}{m} v +\end{align*} +\\] + +where \\(v = \frac{dx}{dt}\\) is the velocity of the mass. + +We can solve this system of ODEs using DiffSol with the following code: + +```rust +# fn main() { +# use std::fs; +use diffsol::{ + DiffSl, CraneliftModule, OdeBuilder, Bdf, OdeSolverState, OdeSolverMethod +}; +use plotly::{ + Plot, Scatter, common::Mode, layout::Layout, layout::Axis +}; +type M = nalgebra::DMatrix; +type CG = CraneliftModule; + +let eqn = DiffSl::::compile(" + k { 1.0 } m { 1.0 } c { 0.1 } + u_i { + x = 1, + v = 0, + } + F_i { + v, + -k/m * x - c/m * v, + } +").unwrap(); +let problem = OdeBuilder::new() + .build_from_eqn(eqn).unwrap(); +let mut solver = Bdf::default(); +let state = OdeSolverState::new(&problem, &solver).unwrap(); +let (ys, ts) = solver.solve(&problem, state, 40.0).unwrap(); + +let x: Vec<_> = ys.row(0).into_iter().copied().collect(); +let time: Vec<_> = ts.into_iter().collect(); + +let x_line = Scatter::new(time.clone(), x).mode(Mode::Lines); + +let mut plot = Plot::new(); +plot.add_trace(x_line); + +let layout = Layout::new() + .x_axis(Axis::new().title("t")) + .y_axis(Axis::new().title("x")); +plot.set_layout(layout); +let plot_html = plot.to_inline_html(Some("sping-mass-system")); +# fs::write("../src/primer/images/spring-mass-system.html", plot_html).expect("Unable to write file"); +# } +``` +{{#include images/spring-mass-system.html}} diff --git a/book/src/primer/the_mass_matrix.md b/book/src/primer/the_mass_matrix.md new file mode 100644 index 00000000..bc77a383 --- /dev/null +++ b/book/src/primer/the_mass_matrix.md @@ -0,0 +1,37 @@ +# DAEs via the Mass Matrix + +Differential-algebraic equations (DAEs) are a generalisation of ordinary differential equations (ODEs) that include algebraic equations, or equations that do not involve derivatives. Algebraic equations can arise in many physical systems and often are used to model constraints on the system, such as conservation laws or other relationships between the state variables. For example, in an electrical circuit, the current flowing into a node must equal the current flowing out of the node, which can be written as an algebraic equation. + +DAEs can be written in the general implicit form: + +\\[ +\mathbf{F}(\mathbf{y}, \mathbf{y}', t) = 0 +\\] + +where \\(\mathbf{y}\\) is the vector of state variables, \\(\mathbf{y}'\\) is the vector of derivatives of the state variables, and \\(\mathbf{F}\\) is a vector-valued function that describes the system of equations. However, for the purposes of this primer and the capabilities of DiffSol, we will focus on a specific form of DAEs called index-1 or semi-explicit DAEs, which can be written as a combination of differential and algebraic equations: + +\\[ +\begin{align*} +\frac{d\mathbf{y}}{dt} &= \mathbf{f}(\mathbf{y}, \mathbf{y}', t) \\\\ +0 = \mathbf{g}(\mathbf{y}, t) +\end{align*} +\\] + +where \\(\mathbf{f}\\) is the vector-valued function that describes the differential equations and \\(\mathbf{g}\\) is the vector-valued function that describes the algebraic equations. The key difference between DAEs and ODEs is that DAEs include algebraic equations that must be satisfied at each time step, in addition to the differential equations that describe the rate of change of the state variables. + +How does this relate to the standard form of an explicit ODE that we have seen before? Recall that an explicit ODE can be written as: + +\\[ +\frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y}, t) +\\] + +Lets update this equation to include a matrix \\(\mathbf{M}\\) that multiplies the derivative term: + +\\[ +M \frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y}, t) +\\] + +When \\(M\\) is the identity matrix (i.e. a matrix with ones along the diagonal), this reduces to the standard form of an explicit ODE. However, when \\(M\\) has diagonal entries that are zero, this introduces algebraic equations into the system and it reduces to the semi-explicit DAE equations show above. The matrix \\(M\\) is called the *mass matrix*. + +Thus, we now have a general form of a set of differential equations, that includes both ODEs and semi-explicit DAEs. This general form is used by DiffSol to allow users to specify a wide range of problems, from simple ODEs to more complex DAEs. In the next section, we will look at a few examples of DAEs and how to solve them using DiffSol and a mass matrix. + diff --git a/book/src/sparse_problems.md b/book/src/sparse_problems.md index df653b4a..3aff5e4b 100644 --- a/book/src/sparse_problems.md +++ b/book/src/sparse_problems.md @@ -1,185 +1 @@ # Sparse problems - -Lets consider a large system of equations that have a jacobian matrix that is sparse. For simplicity we will start with the logistic equation from the ["Specifying the Problem"](./specifying_the_problem.md) section, -but we will duplicate this equation 10 times to create a system of 10 equations. This system will have a jacobian matrix that is a diagonal matrix with 10 diagonal elements, and all other elements are zero. - -Since this system is sparse, we choose a sparse matrix type to represent the jacobian matrix. We will use the `diffsol::SparseColMat` type, which is a thin wrapper around `faer::sparse::SparseColMat`, a sparse compressed sparse column matrix type. - -```rust -# fn main() { -use diffsol::OdeBuilder; -type M = diffsol::SparseColMat; -type V = faer::Col; - -let problem = OdeBuilder::new() - .t0(0.0) - .rtol(1e-6) - .atol([1e-6]) - .p(vec![1.0, 10.0]) - .build_ode::( - |x, p, _t, y| { - for i in 0..10 { - y[i] = p[0] * x[i] * (1.0 - x[i] / p[1]); - } - }, - |x, p, _t, v , y| { - for i in 0..10 { - y[i] = p[0] * v[i] * (1.0 - 2.0 * x[i] / p[1]); - } - }, - |_p, _t| V::from_fn(10, |_| 0.1), - ).unwrap(); -# } -``` - -Note that we have not specified the jacobian itself, but instead we have specified the jacobian multiplied by a vector function \\(f'(y, p, t, v)\\). -DiffSol will use this function to generate a jacobian matrix, and since we have specified a sparse matrix type, DiffSol will attempt to -guess the sparsity pattern of the jacobian matrix and use this to efficiently generate the jacobian matrix. - -To illustrate this, we can calculate the jacobian matrix from the `rhs` function contained in the `problem` object: - -```rust -# use diffsol::OdeBuilder; -use diffsol::{OdeEquations, NonLinearOp, NonLinearOpJacobian, Matrix, ConstantOp}; - -# type M = diffsol::SparseColMat; -# type V = faer::Col; -# -# fn main() { -# let problem = OdeBuilder::new() -# .t0(0.0) -# .rtol(1e-6) -# .atol([1e-6]) -# .p(vec![1.0, 10.0]) -# .build_ode::( -# |x, p, _t, y| { -# for i in 0..10 { -# y[i] = p[0] * x[i] * (1.0 - x[i] / p[1]); -# } -# }, -# |x, p, _t, v , y| { -# for i in 0..10 { -# y[i] = p[0] * v[i] * (1.0 - 2.0 * x[i] / p[1]); -# } -# }, -# |_p, _t| V::from_fn(10, |_| 0.1), -# ).unwrap(); -let t0 = problem.t0; -let y0 = problem.eqn.init().call(t0); -let jacobian = problem.eqn.rhs().jacobian(&y0, t0); -for (i, j, v) in jacobian.triplet_iter() { - println!("({}, {}) = {}", i, j, v); -} -# } -``` - -which will print the jacobian matrix in triplet format: - -``` -(0, 0) = 0.98 -(1, 1) = 0.98 -(2, 2) = 0.98 -(3, 3) = 0.98 -(4, 4) = 0.98 -(5, 5) = 0.98 -(6, 6) = 0.98 -(7, 7) = 0.98 -(8, 8) = 0.98 -(9, 9) = 0.98 -``` - -DiffSol attempts to guess the sparsity pattern of your jacobian matrix by calling the \\(f'(y, p, t, v)\\) function repeatedly with different one-hot vectors \\(v\\) -with a `NaN` value at each index. The output of this function (i.e. which elements are `0` and which are `NaN`) is then used to determine the sparsity pattern of the jacobian matrix. -Due to the fact that for IEEE 754 floating point numbers, `NaN` is propagated through most operations, this method is able to detect which output elements are dependent on which input elements. - -However, this method is not foolproof, and it may fail to detect the correct sparsity pattern in some cases, particularly if values of `v` are used in control-flow statements. -If DiffSol does not detect the correct sparsity pattern, you can manually specify the jacobian. To do this, you need -to implement the [`diffsol::NonLinearOp`](https://docs.rs/diffsol/latest/diffsol/op/trait.NonLinearOp.html) trait for the rhs function. -This is described in more detail in the ["Custom Problem Structs"](./custom_problem_structs.md) section, but is illustrated below. - -```rust -# fn main() { -use std::rc::Rc; -use faer::sparse::{SparseColMat, SymbolicSparseColMat}; -use diffsol::{NonLinearOp, NonLinearOpJacobian, OdeSolverEquations, Op, UnitCallable, ConstantClosure, OdeBuilder}; - -type T = f64; -type V = faer::Col; -type M = diffsol::SparseColMat; - -struct MyProblem { - jacobian: SparseColMat, - p: Rc, -} - -impl MyProblem { - fn new(p: Rc) -> Self { - let mut triplets = Vec::new(); - for i in 0..10 { - triplets.push((i, i, 1.0)); - } - let jacobian = SparseColMat::try_new_from_triplets(10, 10, triplets.as_slice()).unwrap(); - MyProblem { p, jacobian } - } -} - -impl Op for MyProblem { - type T = T; - type V = V; - type M = M; - fn nstates(&self) -> usize { - 10 - } - fn nout(&self) -> usize { - 10 - } - fn nparams(&self) -> usize { - 2 - } -} - -impl NonLinearOp for MyProblem { - fn call_inplace(&self, x: &V, _t: T, y: &mut V) { - for i in 0..10 { - y[i] = self.p[0] * x[i] * (1.0 - x[i] / self.p[1]); - } - } - - -} -impl NonLinearOpJacobian for MyProblem { - fn jac_mul_inplace(&self, x: &V, _t: T, v: &V, y: &mut V) { - for i in 0..10 { - y[i] = self.p[0] * v[i] * (1.0 - 2.0 * x[i] / self.p[1]); - } - } - fn jacobian_inplace(&self, x: &Self::V, _t: Self::T, y: &mut Self::M) { - for i in 0..10 { - let row = y.faer().row_indices()[i]; - y.faer_mut().values_mut()[i] = self.p[0] * (1.0 - 2.0 * x[row] / self.p[1]); - } - } - fn jacobian_sparsity(&self) -> Option> { - Some(self.jacobian.symbolic().to_owned().unwrap()) - } -} - -let p_slice = [1.0, 10.0]; -let p = Rc::new(V::from_fn(p_slice.len(), |i| p_slice[i])); -let rhs = MyProblem::new(p.clone()); - -// use the provided constant closure to define the initial condition -let init_fn = |_p: &V, _t: T| V::from_fn(10, |_| 0.1); -let init = ConstantClosure::new(init_fn, p.clone()); - -// we don't have a mass matrix, root or output functions, so we can set to None -// we still need to give a placeholder type for these, so we use the diffsol::UnitCallable type -let mass: Option> = None; -let root: Option> = None; -let out: Option> = None; - -let eqn = OdeSolverEquations::new(rhs, mass, root, init, out, p); -let _problem = OdeBuilder::new().p(p_slice).build_from_eqn(eqn).unwrap(); -# } -``` - diff --git a/book/src/constant_functions.md b/book/src/specify/custom/constant_functions.md similarity index 100% rename from book/src/constant_functions.md rename to book/src/specify/custom/constant_functions.md diff --git a/book/src/custom_problem_structs.md b/book/src/specify/custom/custom_problem_structs.md similarity index 100% rename from book/src/custom_problem_structs.md rename to book/src/specify/custom/custom_problem_structs.md diff --git a/book/src/linear_functions.md b/book/src/specify/custom/linear_functions.md similarity index 100% rename from book/src/linear_functions.md rename to book/src/specify/custom/linear_functions.md diff --git a/book/src/non_linear_functions.md b/book/src/specify/custom/non_linear_functions.md similarity index 100% rename from book/src/non_linear_functions.md rename to book/src/specify/custom/non_linear_functions.md diff --git a/book/src/putting_it_all_together.md b/book/src/specify/custom/putting_it_all_together.md similarity index 100% rename from book/src/putting_it_all_together.md rename to book/src/specify/custom/putting_it_all_together.md diff --git a/book/src/specify/diffsl.md b/book/src/specify/diffsl.md new file mode 100644 index 00000000..56cad6ad --- /dev/null +++ b/book/src/specify/diffsl.md @@ -0,0 +1,59 @@ +# DiffSL + +Thus far we have used Rust code to specify the problem we want to solve. This is fine if you are using DiffSol from Rust, but what if you want to use DiffSol from a higher-level language like Python or R? +For this usecase we have designed a Domain Specific Language (DSL) called DiffSL that can be used to specify the problem. DiffSL is not a general purpose language but is tightly constrained to +the specification of a system of ordinary differential equations. It features a relatively simple syntax that consists of writing a series of tensors (dense or sparse) that represent the equations of the system. +For more detail on the syntax of DiffSL see the [DiffSL book](https://martinjrobins.github.io/diffsl/). This section will focus on how to use DiffSL to specify a problem in DiffSol. + + +## DiffSL Context + +The main struct that is used to specify a problem in DiffSL is the [`DiffSl`](https://docs.rs/diffsol/latest/diffsol/ode_solver/diffsl/struct.DiffSl.html) struct. Creating this struct +Just-In-Time (JIT) compiles your DiffSL code into a form that can be executed efficiently by DiffSol. + +```rust +# fn main() { +use diffsol::{DiffSl, CraneliftModule}; +type M = nalgebra::DMatrix; +type CG = CraneliftModule; + +let eqn = DiffSl::::compile(" + in = [r, k] + r { 1.0 } + k { 1.0 } + u { 0.1 } + F { r * u * (1.0 - u / k) } +").unwrap(); +# } +``` + +The `CG` parameter specifies the backend that you want to use to compile the DiffSL code. The `CraneliftModule` backend is the default backend and is behind the `diffsl` feature flag. If you want to use the faster LLVM backend you can use the `LlvmModule` backend, which is behind one of the `diffsl-llvm*` feature flags, depending on the version of LLVM you have installed. + +Once you have created the `DiffSl` struct you can use it to create a problem using the `build_from_eqn` method on the [`OdeBuilder`](https://docs.rs/diffsol/latest/diffsol/ode_solver/builder/struct.OdeBuilder.html) struct. + + +```rust +# fn main() { +# use diffsol::{DiffSl, CraneliftModule}; +use diffsol::{OdeBuilder, Bdf, OdeSolverMethod, OdeSolverState}; +# type M = nalgebra::DMatrix; +# type CG = CraneliftModule; + + +# let eqn = DiffSl::::compile(" +# in = [r, k] +# r { 1.0 } +# k { 1.0 } +# u { 0.1 } +# F { r * u * (1.0 - u / k) } +# ").unwrap(); +let problem = OdeBuilder::new() +.rtol(1e-6) +.p([1.0, 10.0]) +.build_from_eqn(eqn).unwrap(); +let mut solver = Bdf::default(); +let t = 0.4; +let state = OdeSolverState::new(&problem, &solver).unwrap(); +let _soln = solver.solve(&problem, state, t).unwrap(); +# } +``` diff --git a/book/src/forward_sensitivity.md b/book/src/specify/forward_sensitivity.md similarity index 100% rename from book/src/forward_sensitivity.md rename to book/src/specify/forward_sensitivity.md diff --git a/book/src/mass_matrix.md b/book/src/specify/mass_matrix.md similarity index 100% rename from book/src/mass_matrix.md rename to book/src/specify/mass_matrix.md diff --git a/book/src/ode_equations.md b/book/src/specify/ode_equations.md similarity index 100% rename from book/src/ode_equations.md rename to book/src/specify/ode_equations.md diff --git a/book/src/root_finding.md b/book/src/specify/root_finding.md similarity index 100% rename from book/src/root_finding.md rename to book/src/specify/root_finding.md diff --git a/book/src/specify/sparse_problems.md b/book/src/specify/sparse_problems.md new file mode 100644 index 00000000..df653b4a --- /dev/null +++ b/book/src/specify/sparse_problems.md @@ -0,0 +1,185 @@ +# Sparse problems + +Lets consider a large system of equations that have a jacobian matrix that is sparse. For simplicity we will start with the logistic equation from the ["Specifying the Problem"](./specifying_the_problem.md) section, +but we will duplicate this equation 10 times to create a system of 10 equations. This system will have a jacobian matrix that is a diagonal matrix with 10 diagonal elements, and all other elements are zero. + +Since this system is sparse, we choose a sparse matrix type to represent the jacobian matrix. We will use the `diffsol::SparseColMat` type, which is a thin wrapper around `faer::sparse::SparseColMat`, a sparse compressed sparse column matrix type. + +```rust +# fn main() { +use diffsol::OdeBuilder; +type M = diffsol::SparseColMat; +type V = faer::Col; + +let problem = OdeBuilder::new() + .t0(0.0) + .rtol(1e-6) + .atol([1e-6]) + .p(vec![1.0, 10.0]) + .build_ode::( + |x, p, _t, y| { + for i in 0..10 { + y[i] = p[0] * x[i] * (1.0 - x[i] / p[1]); + } + }, + |x, p, _t, v , y| { + for i in 0..10 { + y[i] = p[0] * v[i] * (1.0 - 2.0 * x[i] / p[1]); + } + }, + |_p, _t| V::from_fn(10, |_| 0.1), + ).unwrap(); +# } +``` + +Note that we have not specified the jacobian itself, but instead we have specified the jacobian multiplied by a vector function \\(f'(y, p, t, v)\\). +DiffSol will use this function to generate a jacobian matrix, and since we have specified a sparse matrix type, DiffSol will attempt to +guess the sparsity pattern of the jacobian matrix and use this to efficiently generate the jacobian matrix. + +To illustrate this, we can calculate the jacobian matrix from the `rhs` function contained in the `problem` object: + +```rust +# use diffsol::OdeBuilder; +use diffsol::{OdeEquations, NonLinearOp, NonLinearOpJacobian, Matrix, ConstantOp}; + +# type M = diffsol::SparseColMat; +# type V = faer::Col; +# +# fn main() { +# let problem = OdeBuilder::new() +# .t0(0.0) +# .rtol(1e-6) +# .atol([1e-6]) +# .p(vec![1.0, 10.0]) +# .build_ode::( +# |x, p, _t, y| { +# for i in 0..10 { +# y[i] = p[0] * x[i] * (1.0 - x[i] / p[1]); +# } +# }, +# |x, p, _t, v , y| { +# for i in 0..10 { +# y[i] = p[0] * v[i] * (1.0 - 2.0 * x[i] / p[1]); +# } +# }, +# |_p, _t| V::from_fn(10, |_| 0.1), +# ).unwrap(); +let t0 = problem.t0; +let y0 = problem.eqn.init().call(t0); +let jacobian = problem.eqn.rhs().jacobian(&y0, t0); +for (i, j, v) in jacobian.triplet_iter() { + println!("({}, {}) = {}", i, j, v); +} +# } +``` + +which will print the jacobian matrix in triplet format: + +``` +(0, 0) = 0.98 +(1, 1) = 0.98 +(2, 2) = 0.98 +(3, 3) = 0.98 +(4, 4) = 0.98 +(5, 5) = 0.98 +(6, 6) = 0.98 +(7, 7) = 0.98 +(8, 8) = 0.98 +(9, 9) = 0.98 +``` + +DiffSol attempts to guess the sparsity pattern of your jacobian matrix by calling the \\(f'(y, p, t, v)\\) function repeatedly with different one-hot vectors \\(v\\) +with a `NaN` value at each index. The output of this function (i.e. which elements are `0` and which are `NaN`) is then used to determine the sparsity pattern of the jacobian matrix. +Due to the fact that for IEEE 754 floating point numbers, `NaN` is propagated through most operations, this method is able to detect which output elements are dependent on which input elements. + +However, this method is not foolproof, and it may fail to detect the correct sparsity pattern in some cases, particularly if values of `v` are used in control-flow statements. +If DiffSol does not detect the correct sparsity pattern, you can manually specify the jacobian. To do this, you need +to implement the [`diffsol::NonLinearOp`](https://docs.rs/diffsol/latest/diffsol/op/trait.NonLinearOp.html) trait for the rhs function. +This is described in more detail in the ["Custom Problem Structs"](./custom_problem_structs.md) section, but is illustrated below. + +```rust +# fn main() { +use std::rc::Rc; +use faer::sparse::{SparseColMat, SymbolicSparseColMat}; +use diffsol::{NonLinearOp, NonLinearOpJacobian, OdeSolverEquations, Op, UnitCallable, ConstantClosure, OdeBuilder}; + +type T = f64; +type V = faer::Col; +type M = diffsol::SparseColMat; + +struct MyProblem { + jacobian: SparseColMat, + p: Rc, +} + +impl MyProblem { + fn new(p: Rc) -> Self { + let mut triplets = Vec::new(); + for i in 0..10 { + triplets.push((i, i, 1.0)); + } + let jacobian = SparseColMat::try_new_from_triplets(10, 10, triplets.as_slice()).unwrap(); + MyProblem { p, jacobian } + } +} + +impl Op for MyProblem { + type T = T; + type V = V; + type M = M; + fn nstates(&self) -> usize { + 10 + } + fn nout(&self) -> usize { + 10 + } + fn nparams(&self) -> usize { + 2 + } +} + +impl NonLinearOp for MyProblem { + fn call_inplace(&self, x: &V, _t: T, y: &mut V) { + for i in 0..10 { + y[i] = self.p[0] * x[i] * (1.0 - x[i] / self.p[1]); + } + } + + +} +impl NonLinearOpJacobian for MyProblem { + fn jac_mul_inplace(&self, x: &V, _t: T, v: &V, y: &mut V) { + for i in 0..10 { + y[i] = self.p[0] * v[i] * (1.0 - 2.0 * x[i] / self.p[1]); + } + } + fn jacobian_inplace(&self, x: &Self::V, _t: Self::T, y: &mut Self::M) { + for i in 0..10 { + let row = y.faer().row_indices()[i]; + y.faer_mut().values_mut()[i] = self.p[0] * (1.0 - 2.0 * x[row] / self.p[1]); + } + } + fn jacobian_sparsity(&self) -> Option> { + Some(self.jacobian.symbolic().to_owned().unwrap()) + } +} + +let p_slice = [1.0, 10.0]; +let p = Rc::new(V::from_fn(p_slice.len(), |i| p_slice[i])); +let rhs = MyProblem::new(p.clone()); + +// use the provided constant closure to define the initial condition +let init_fn = |_p: &V, _t: T| V::from_fn(10, |_| 0.1); +let init = ConstantClosure::new(init_fn, p.clone()); + +// we don't have a mass matrix, root or output functions, so we can set to None +// we still need to give a placeholder type for these, so we use the diffsol::UnitCallable type +let mass: Option> = None; +let root: Option> = None; +let out: Option> = None; + +let eqn = OdeSolverEquations::new(rhs, mass, root, init, out, p); +let _problem = OdeBuilder::new().p(p_slice).build_from_eqn(eqn).unwrap(); +# } +``` + diff --git a/book/src/specifying_the_problem.md b/book/src/specify/specifying_the_problem.md similarity index 68% rename from book/src/specifying_the_problem.md rename to book/src/specify/specifying_the_problem.md index 5efaecb6..f6f9e04a 100644 --- a/book/src/specifying_the_problem.md +++ b/book/src/specify/specifying_the_problem.md @@ -1,6 +1,6 @@ -# Specifying the problem +# DiffSol APIs for specifying problems -Most of the DiffSol user-facing API revolves around specifying the problem you want to solve, thus a large part of this book will be dedicated to explaining how to specify a problem. +Most of the DiffSol user-facing API revolves around specifying the problem you want to solve, thus a large part of this book will be dedicated to explaining how to specify a problem. All the examples presented in [the primer](../primer/modelling_with_odes.md) used the DiffSL DSL to specify the problem, but DiffSol also provides a pure Rust API for specifying problems. This API is sometimes more verbose than the DSL, but is more flexible, more performant, and easier to use if you have a model already written in Rust or another language that you can easily port or call from Rust. ## ODE equations @@ -30,9 +30,8 @@ DiffSol has three main APIs for specifying problems: 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 [`DiffSlContext`](https://docs.rs/diffsol/latest/diffsol/ode_solver/diffsl/struct.DiffSlContext.html) struct, where the user can specify the functions above using the [DiffSl](https://martinjrobins.github.io/diffsl/) - Domain Specific Language (DSL). This API requires a local LLVM installation, and is behind a feature flag, but has the best API if you want to use DiffSL from a higher-level language like Python or R - while still having the performance of Rust. +- The [`DiffSlContext`](https://docs.rs/diffsol/latest/diffsol/ode_solver/diffsl/struct.DiffSlContext.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. diff --git a/book/theme/head.hbs b/book/theme/head.hbs new file mode 100644 index 00000000..a6d1d16c --- /dev/null +++ b/book/theme/head.hbs @@ -0,0 +1 @@ + diff --git a/build.rs b/build.rs index 9f8d4961..efd6aa95 100644 --- a/build.rs +++ b/build.rs @@ -188,54 +188,52 @@ mod if_sundials { } } -#[cfg(feature = "sundials")] fn main() -> Result<(), String> { - use if_sundials::*; - - // return if we are not using the sundials features - if !cfg!(feature = "sundials") { - return Ok(()); - } - let sundials = Library::new_sundials(); - let suitesparse = Library::new_suitesparse(); - - // generate sundials bindings - let bindings_rs = generate_bindings(&suitesparse, &sundials)?; - let major = get_sundials_version_major(bindings_rs).unwrap_or(5); - println!("cargo:rustc-cfg=sundials_version_major=\"{}\"", major); - println!("cargo::rustc-check-cfg=cfg(sundials_version_major, values(\"5\", \"6\", \"7\"))"); - - // compile sundials benches - let mut files = compile_benches(&sundials, &suitesparse); - files.push("benches/sundials_benches.rs".to_string()); - files.push("build.rs".to_string()); - files.push("benches/idaHeat2d_klu_v5.c".to_string()); - files.push("benches/idaHeat2d_klu_v6.c".to_string()); - files.push("benches/idaFoodWeb_bnd_v5.c".to_string()); - files.push("benches/idaFoodWeb_bnd_v6.c".to_string()); - for name in files.into_iter() { - println!("cargo:rerun-if-changed={}", name); - } + #[cfg(feature = "sundials")] + { + use if_sundials::*; + // return if we are not using the sundials features + if !cfg!(feature = "sundials") { + return Ok(()); + } + let sundials = Library::new_sundials(); + let suitesparse = Library::new_suitesparse(); + + // generate sundials bindings + let bindings_rs = generate_bindings(&suitesparse, &sundials)?; + let major = get_sundials_version_major(bindings_rs).unwrap_or(5); + println!("cargo:rustc-cfg=sundials_version_major=\"{}\"", major); + println!("cargo::rustc-check-cfg=cfg(sundials_version_major, values(\"5\", \"6\", \"7\"))"); + + // compile sundials benches + let mut files = compile_benches(&sundials, &suitesparse); + files.push("benches/sundials_benches.rs".to_string()); + files.push("build.rs".to_string()); + files.push("benches/idaHeat2d_klu_v5.c".to_string()); + files.push("benches/idaHeat2d_klu_v6.c".to_string()); + files.push("benches/idaFoodWeb_bnd_v5.c".to_string()); + files.push("benches/idaFoodWeb_bnd_v6.c".to_string()); + for name in files.into_iter() { + println!("cargo:rerun-if-changed={}", name); + } - // link to sundials - if let Some(dir) = sundials.lib.as_ref() { - println!("cargo:rustc-link-search=native={}", dir); - } - let library_type = "dylib"; - for lib_name in LINK_SUNDIALS_LIBRARIES { - println!( - "cargo:rustc-link-lib={}=sundials_{}", - library_type, lib_name - ); - } - // link to klu - if let Some(dir) = suitesparse.lib.as_ref() { - println!("cargo:rustc-link-search=native={}", dir); - println!("cargo:rustc-env=LD_LIBRARY_PATH={}", dir); + // link to sundials + if let Some(dir) = sundials.lib.as_ref() { + println!("cargo:rustc-link-search=native={}", dir); + } + let library_type = "dylib"; + for lib_name in LINK_SUNDIALS_LIBRARIES { + println!( + "cargo:rustc-link-lib={}=sundials_{}", + library_type, lib_name + ); + } + // link to klu + if let Some(dir) = suitesparse.lib.as_ref() { + println!("cargo:rustc-link-search=native={}", dir); + println!("cargo:rustc-env=LD_LIBRARY_PATH={}", dir); + } + println!("cargo:rustc-link-lib={}=klu", library_type); } - println!("cargo:rustc-link-lib={}=klu", library_type); Ok(()) } - -#[cfg(not(feature = "sundials"))] -fn main() {} diff --git a/src/error.rs b/src/error.rs index ac80ea30..8940de76 100644 --- a/src/error.rs +++ b/src/error.rs @@ -69,7 +69,7 @@ pub enum OdeSolverError { StepSizeTooSmall { time: f64 }, #[error("Sensitivity requested but equations do not support it")] SensitivityNotSupported, - #[error("Failed to get mutable reference to equations, is there a solver created with this problem?")] + #[error("Failed to get mutable reference to equations. If there is a solver created with this problem, call solver.take_state() to release the problem")] FailedToGetMutableReference, #[error("Builder error: {0}")] BuilderError(String), diff --git a/src/jacobian/mod.rs b/src/jacobian/mod.rs index 36dd49ea..ddd690cf 100644 --- a/src/jacobian/mod.rs +++ b/src/jacobian/mod.rs @@ -86,6 +86,7 @@ gen_find_non_zeros_linear!( LinearOpTranspose ); +#[derive(Clone)] pub struct JacobianColoring { dst_indices_per_color: Vec<::Index>, src_indices_per_color: Vec<::Index>, diff --git a/src/lib.rs b/src/lib.rs index 6d3fb6a6..9bd6642b 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -126,19 +126,7 @@ //! - For vectors: [Vector], [VectorIndex], [VectorView], [VectorViewMut], and [VectorCommon]. //! -#[cfg(feature = "diffsl-llvm13")] -pub extern crate diffsl13_0 as diffsl; -#[cfg(feature = "diffsl-llvm14")] -pub extern crate diffsl14_0 as diffsl; -#[cfg(feature = "diffsl-llvm15")] -pub extern crate diffsl15_0 as diffsl; -#[cfg(feature = "diffsl-llvm16")] -pub extern crate diffsl16_0 as diffsl; -#[cfg(feature = "diffsl-llvm17")] -pub extern crate diffsl17_0 as diffsl; -#[cfg(feature = "diffsl-cranelift")] -pub extern crate diffsl_no_llvm as diffsl; -#[cfg(feature = "diffsl-cranelift")] +#[cfg(feature = "diffsl")] pub use diffsl::CraneliftModule; #[cfg(feature = "diffsl-llvm")] pub use diffsl::LlvmModule; diff --git a/src/ode_solver/bdf.rs b/src/ode_solver/bdf.rs index 53fe8f31..28c5205f 100644 --- a/src/ode_solver/bdf.rs +++ b/src/ode_solver/bdf.rs @@ -794,6 +794,9 @@ where self.state.as_ref().map(|state| state.as_ref()) } fn take_state(&mut self) -> Option> { + self.ode_problem = None; + self.op = None; + self.s_op = None; Option::take(&mut self.state) } @@ -889,7 +892,21 @@ where let mut convergence_fail = false; if self.is_state_modified { + // reinitalise root finder if needed + if let Some(root_fn) = problem.eqn.root() { + let state = self.state.as_ref().unwrap(); + self.root_finder + .as_ref() + .unwrap() + .init(&root_fn, &state.y, state.t); + } + // reinitialise diff matrix self.initialise_to_first_order(); + + // reinitialise tstop if needed + if let Some(t_stop) = self.tstop { + self.set_stop_time(t_stop)?; + } } self._predict_forward(); @@ -1206,7 +1223,7 @@ mod test { }, tests::{ test_checkpointing, test_interpolate, test_no_set_problem, test_ode_solver, - test_ode_solver_adjoint, test_ode_solver_no_sens, test_state_mut, + test_ode_solver_adjoint, test_ode_solver_no_sens, test_param_sweep, test_state_mut, test_state_mut_on_problem, }, }, @@ -1678,4 +1695,39 @@ mod test { let y = test_ode_solver_no_sens(&mut s, &problem, soln, None, false); assert!(abs(y[0] - 0.6) < 1e-6, "y[0] = {}", y[0]); } + + #[test] + fn test_param_sweep_bdf() { + let s = Bdf::default(); + let (problem, _soln) = exponential_decay_problem::(false); + let mut ps = Vec::new(); + for y0 in (1..10).map(f64::from) { + ps.push(nalgebra::DVector::::from_vec(vec![0.1, y0])); + } + test_param_sweep(s, problem, ps); + } + + #[cfg(feature = "diffsl")] + #[test] + fn test_ball_bounce_bdf() { + type M = nalgebra::DMatrix; + type LS = crate::NalgebraLU; + type Nls = crate::NewtonNonlinearSolver; + type Eqn = crate::DiffSl; + let s = Bdf::::default(); + let (x, v, t) = crate::ode_solver::tests::test_ball_bounce(s); + + let expected_x = [ + 0.003751514915514589, + 0.00750117409999241, + 0.015370589755655079, + ]; + let expected_v = [11.202428570923361, 11.19914432101355, 11.192247396202946]; + let expected_t = [1.4281779078441663, 1.4285126937676944, 1.4292157442071036]; + for (i, ((x, v), t)) in x.iter().zip(v.iter()).zip(t.iter()).enumerate() { + assert!((x - expected_x[i]).abs() < 1e-4); + assert!((v - expected_v[i]).abs() < 1e-4); + assert!((t - expected_t[i]).abs() < 1e-4); + } + } } diff --git a/src/ode_solver/diffsl.rs b/src/ode_solver/diffsl.rs index d15c814e..0d42b0e7 100644 --- a/src/ode_solver/diffsl.rs +++ b/src/ode_solver/diffsl.rs @@ -14,37 +14,6 @@ pub type T = f64; /// Context for the ODE equations specified using the [DiffSL language](https://martinjrobins.github.io/diffsl/). /// /// This contains the compiled code and the data structures needed to evaluate the ODE equations. -/// Note that the example below uses the LLVM backend (requires one of the `diffsl-llvm` features), -/// but the Cranelift backend can also be used using -/// `diffsl::CraneliftModule` instead of `diffsl::LlvmModule` (requires `diffsl-cranelift` feature). -/// -/// # Example -/// -/// ```rust -/// use diffsol::{OdeBuilder, Bdf, OdeSolverState, OdeSolverMethod, DiffSl, diffsl::LlvmModule}; -/// -/// // dy/dt = -ay -/// // y(0) = 1 -/// let eqn = DiffSl::, LlvmModule>::compile(" -/// in = [a] -/// a { 1 } -/// u { 1.0 } -/// F { -a*u } -/// out { u } -/// ").unwrap(); -/// let problem = OdeBuilder::new() -/// .rtol(1e-6) -/// .p([0.1]) -/// .build_from_eqn(eqn).unwrap(); -/// let mut solver = Bdf::default(); -/// let t = 0.4; -/// let state = OdeSolverState::new(&problem, &solver).unwrap(); -/// solver.set_problem(state, &problem); -/// while solver.state().unwrap().t <= t { -/// solver.step().unwrap(); -/// } -/// let y = solver.interpolate(t); -/// ``` pub struct DiffSlContext, CG: CodegenModule> { compiler: Compiler, data: RefCell>, @@ -53,6 +22,7 @@ pub struct DiffSlContext, CG: CodegenModule> { nstates: usize, nroots: usize, nparams: usize, + has_mass: bool, nout: usize, } @@ -62,7 +32,7 @@ impl, CG: CodegenModule> DiffSlContext { pub fn new(text: &str) -> Result { let compiler = Compiler::from_discrete_str(text).map_err(|e| DiffsolError::Other(e.to_string()))?; - let (nstates, nparams, nout, _ndata, nroots) = compiler.get_dims(); + let (nstates, nparams, nout, _ndata, nroots, has_mass) = compiler.get_dims(); let data = RefCell::new(compiler.get_new_data()); let ddata = RefCell::new(compiler.get_new_data()); let tmp = RefCell::new(M::V::zeros(nstates)); @@ -76,13 +46,14 @@ impl, CG: CodegenModule> DiffSlContext { tmp, nroots, nout, + has_mass, }) } pub fn recompile(&mut self, text: &str) -> Result<(), DiffsolError> { self.compiler = Compiler::from_discrete_str(text).map_err(|e| DiffsolError::Other(e.to_string()))?; - let (nstates, nparams, nout, _ndata, nroots) = self.compiler.get_dims(); + let (nstates, nparams, nout, _ndata, nroots, has_mass) = self.compiler.get_dims(); self.data = RefCell::new(self.compiler.get_new_data()); self.ddata = RefCell::new(self.compiler.get_new_data()); self.tmp = RefCell::new(M::V::zeros(nstates)); @@ -90,6 +61,7 @@ impl, CG: CodegenModule> DiffSlContext { self.nstates = nstates; self.nout = nout; self.nroots = nroots; + self.has_mass = has_mass; Ok(()) } } @@ -400,7 +372,7 @@ impl, CG: CodegenModule> OdeEquations for DiffSl { } fn mass(&self) -> Option> { - self.context.compiler.has_mass().then_some(DiffSlMass(self)) + self.context.has_mass.then_some(DiffSlMass(self)) } fn root(&self) -> Option> { diff --git a/src/ode_solver/method.rs b/src/ode_solver/method.rs index bfab1b4e..51176845 100644 --- a/src/ode_solver/method.rs +++ b/src/ode_solver/method.rs @@ -308,11 +308,10 @@ where while adjoint_solver.step()? != OdeSolverStopReason::TstopReached {} // correct the adjoint solution for the initial conditions + let adjoint_problem = adjoint_solver.problem().unwrap().clone(); let mut state = adjoint_solver.take_state().unwrap(); let state_mut = state.as_mut(); - adjoint_solver - .problem() - .unwrap() + adjoint_problem .eqn .correct_sg_for_init(t0, state_mut.s, state_mut.sg); diff --git a/src/ode_solver/mod.rs b/src/ode_solver/mod.rs index 504a43b6..9265b165 100644 --- a/src/ode_solver/mod.rs +++ b/src/ode_solver/mod.rs @@ -33,9 +33,10 @@ mod tests { use crate::matrix::Matrix; use crate::op::unit::UnitCallable; use crate::{ - op::OpStatistics, NonLinearOpJacobian, OdeEquations, OdeEquationsAdjoint, - OdeEquationsImplicit, OdeEquationsRef, OdeEquationsSens, OdeSolverMethod, OdeSolverProblem, - OdeSolverState, OdeSolverStopReason, + op::OpStatistics, CraneliftModule, DenseMatrix, DiffSl, MatrixCommon, NonLinearOpJacobian, + OdeBuilder, OdeEquations, OdeEquationsAdjoint, OdeEquationsImplicit, OdeEquationsRef, + OdeEquationsSens, OdeSolverMethod, OdeSolverProblem, OdeSolverState, OdeSolverStopReason, + VectorView, }; use crate::{ConstantOp, DefaultDenseMatrix, DefaultSolver, NonLinearOp, Op, Vector}; use num_traits::One; @@ -297,11 +298,10 @@ mod tests { while adjoint_solver.state().unwrap().t.abs() > t0 { adjoint_solver.step().unwrap(); } + let adjoint_problem = adjoint_solver.problem().unwrap().clone(); let mut state = adjoint_solver.take_state().unwrap(); let state_mut = state.as_mut(); - adjoint_solver - .problem() - .unwrap() + adjoint_problem .eqn .correct_sg_for_init(t0, state_mut.s, state_mut.sg); @@ -516,6 +516,86 @@ mod tests { ); } + #[cfg(feature = "diffsl")] + pub fn test_ball_bounce(mut solver: Method) -> (Vec, Vec, Vec) + where + M: Matrix, + M: DefaultSolver, + M::V: DefaultDenseMatrix, + Method: OdeSolverMethod>, + { + let eqn = DiffSl::compile( + " + g { 9.81 } h { 10.0 } + u_i { + x = h, + v = 0, + } + F_i { + v, + -g, + } + stop { + x, + } + ", + ) + .unwrap(); + + let e = 0.8; + let problem = OdeBuilder::new().build_from_eqn(eqn).unwrap(); + let state = OdeSolverState::new(&problem, &solver).unwrap(); + solver.set_problem(state, &problem).unwrap(); + + let final_time = 2.5; + + // solve and apply the remaining doses + solver.set_stop_time(final_time).unwrap(); + loop { + match solver.step() { + Ok(OdeSolverStopReason::InternalTimestep) => (), + Ok(OdeSolverStopReason::RootFound(t)) => { + // get the state when the event occurred + let mut y = solver.interpolate(t).unwrap(); + + // update the velocity of the ball + y[1] *= -e; + + // make sure the ball is above the ground + y[0] = y[0].max(f64::EPSILON); + + // set the state to the updated state + solver.state_mut().unwrap().y.copy_from(&y); + solver.state_mut().unwrap().dy[0] = y[1]; + *solver.state_mut().unwrap().t = t; + + break; + } + Ok(OdeSolverStopReason::TstopReached) => break, + Err(_) => panic!("unexpected solver error"), + } + } + // do three more steps after the 1st bound and many sure they are correct + let mut x = vec![]; + let mut v = vec![]; + let mut t = vec![]; + for _ in 0..3 { + let ret = solver.step(); + x.push(solver.state().unwrap().y[0]); + v.push(solver.state().unwrap().y[1]); + t.push(solver.state().unwrap().t); + match ret { + Ok(OdeSolverStopReason::InternalTimestep) => (), + Ok(OdeSolverStopReason::RootFound(_)) => { + panic!("should be an internal timestep but found a root") + } + Ok(OdeSolverStopReason::TstopReached) => break, + _ => panic!("should be an internal timestep"), + } + } + (x, v, t) + } + pub fn test_checkpointing( mut solver1: Method, mut solver2: Method, @@ -563,6 +643,35 @@ mod tests { } } + pub fn test_param_sweep( + mut s: Method, + mut problem: OdeSolverProblem, + ps: Vec, + ) where + Method: OdeSolverMethod, + Eqn: OdeEquationsImplicit, + Eqn::M: DefaultSolver, + Eqn::V: DefaultDenseMatrix, + { + let mut old_soln = None; + for p in ps { + problem.set_params(p).unwrap(); + let state = OdeSolverState::new(&problem, &s).unwrap(); + let (ys, _ts) = s.solve(&problem, state, Eqn::T::from(10.0)).unwrap(); + // check that the new solution is different from the old one + if let Some(old_soln) = &mut old_soln { + let new_soln = ys.column(ys.ncols() - 1).into_owned(); + let diff = (new_soln - &*old_soln) + .squared_norm(old_soln, &problem.atol, problem.rtol) + .sqrt(); + assert!(diff > Eqn::T::from(1.0e-6), "diff: {}", diff); + } + old_soln = Some(ys.column(ys.ncols() - 1).into_owned()); + s.take_state().unwrap(); + assert!(s.problem().is_none()); + } + } + pub fn test_state_mut_on_problem( mut s: Method, problem: OdeSolverProblem, diff --git a/src/ode_solver/sdirk.rs b/src/ode_solver/sdirk.rs index 38d56e9b..072ddac7 100644 --- a/src/ode_solver/sdirk.rs +++ b/src/ode_solver/sdirk.rs @@ -466,6 +466,9 @@ where } fn take_state(&mut self) -> Option> { + self.problem = None; + self.op = None; + self.s_op = None; Option::take(&mut self.state) } @@ -544,6 +547,23 @@ where } let n = self.state.as_ref().unwrap().y.len(); + if self.is_state_mutated { + // reinitalise root finder if needed + if let Some(root_fn) = self.problem.as_ref().unwrap().eqn.root() { + let state = self.state.as_ref().unwrap(); + self.root_finder + .as_ref() + .unwrap() + .init(&root_fn, &state.y, state.t); + } + // reinitialise tstop if needed + if let Some(t_stop) = self.tstop { + self.set_stop_time(t_stop)?; + } + + self.is_state_mutated = false; + } + // optionally do the first step let start = if self.is_sdirk { 0 } else { 1 }; let mut updated_jacobian = false; @@ -816,8 +836,6 @@ where let new_h = self._update_step_size(factor)?; self._jacobian_updates(new_h, SolverState::StepSuccess); - self.is_state_mutated = false; - // update statistics self.statistics.number_of_linear_solver_setups = self.op.as_ref().unwrap().number_of_jac_evals(); @@ -1076,7 +1094,7 @@ mod test { }, tests::{ test_checkpointing, test_interpolate, test_no_set_problem, test_ode_solver, - test_ode_solver_adjoint, test_ode_solver_no_sens, test_state_mut, + test_ode_solver_adjoint, test_ode_solver_no_sens, test_param_sweep, test_state_mut, test_state_mut_on_problem, }, }, @@ -1362,4 +1380,33 @@ mod test { let y = test_ode_solver_no_sens(&mut s, &problem, soln, None, false); assert!(abs(y[0] - 0.6) < 1e-6, "y[0] = {}", y[0]); } + + #[test] + fn test_param_sweep_tr_bdf2() { + let s = Sdirk::tr_bdf2(); + let (problem, _soln) = exponential_decay_problem::(false); + let mut ps = Vec::new(); + for y0 in (1..10).map(f64::from) { + ps.push(nalgebra::DVector::::from_vec(vec![0.1, y0])); + } + test_param_sweep(s, problem, ps); + } + + #[cfg(feature = "diffsl")] + #[test] + fn test_ball_bounce_tr_bdf2() { + type M = nalgebra::DMatrix; + type LS = crate::NalgebraLU; + type Eqn = crate::DiffSl; + let s = Sdirk::::tr_bdf2(); + let (x, v, t) = crate::ode_solver::tests::test_ball_bounce(s); + let expected_x = [6.375884661615263]; + let expected_v = [0.6878538646461059]; + let expected_t = [2.5]; + for (i, ((x, v), t)) in x.iter().zip(v.iter()).zip(t.iter()).enumerate() { + assert!((x - expected_x[i]).abs() < 1e-4); + assert!((v - expected_v[i]).abs() < 1e-4); + assert!((t - expected_t[i]).abs() < 1e-4); + } + } } diff --git a/src/vector/mod.rs b/src/vector/mod.rs index 2354bc17..3b3f25e1 100644 --- a/src/vector/mod.rs +++ b/src/vector/mod.rs @@ -13,7 +13,7 @@ mod nalgebra_serial; #[cfg(feature = "sundials")] pub mod sundials; -pub trait VectorIndex: Sized + Index + Debug { +pub trait VectorIndex: Sized + Index + Debug + Clone { fn zeros(len: IndexType) -> Self; fn len(&self) -> IndexType; fn clone_as_vec(&self) -> Vec; diff --git a/src/vector/sundials.rs b/src/vector/sundials.rs index 2da18a69..cbe9a723 100644 --- a/src/vector/sundials.rs +++ b/src/vector/sundials.rs @@ -104,7 +104,7 @@ impl SundialsVectorView<'_> { } } -#[derive(Debug)] +#[derive(Debug, Clone)] pub struct SundialsIndexVector(Vec); impl Display for SundialsIndexVector {