Tip

An interactive online version of this notebook is available, which can be accessed via

---

Alternatively, you may download this notebook and run it offline.

Attention

You are viewing this notebook on the latest version of the documentation, where these notebooks may not be compatible with the stable release of PyBaMM since they can contain features that are not yet released. We recommend viewing these notebooks from the stable version of the documentation. To install the latest version of PyBaMM that is compatible with the latest notebooks, build PyBaMM from source.

Creating a simple ODE model

In this series of notebooks, we will run through the steps involved in creating a new model within pybamm. Before using pybamm we recommend following the Getting Started guides.

In this notebook we create and solve the following simple ODE model:

\[\frac{\textrm{d}x}{\textrm{d} t} = 4x - 2y, \quad x(0) = 1,\]

\[\frac{\textrm{d}y}{\textrm{d} t} = 3x - y, \quad y(0) = 2.\]

We begin by importing the PyBaMM library into this notebook, along with NumPy and Matplotlib, which we use for plotting:

%pip install "pybamm[plot,cite]" -q    # install PyBaMM if it is not installed
import pybamm
import numpy as np
import matplotlib.pyplot as plt

Note: you may need to restart the kernel to use updated packages.

Setting up the model

We first initialise the model using the BaseModel class. This sets up the required structure for our model.

model = pybamm.BaseModel()

Next, we define the variables in the model using the Variable class. In more complicated models we can give the variables more informative string names, but here we simply name the variables “x” and “y”

x = pybamm.Variable("x")
y = pybamm.Variable("y")

We can now use the symbols we have created for our variables to write out our governing equations. Note that the governing equations must be provied in the explicit form d/dt = rhs since pybamm only stores the right hand side (rhs) and assumes that the left hand side is the time derivative.

dxdt = 4 * x - 2 * y
dydt = 3 * x - y

The governing equations must then be added to the dictionary model.rhs. The dictionary stores key and item pairs, where the key is the variable which is governed by the equation stored in the corresponding item. Note that the keys are the symbols that represent the variables and are not the variable names (e.g. the key is x, not the string “x”).

model.rhs = {x: dxdt, y: dydt}

The initial conditions are also stored in a dictionary, model.initial_conditions, which again uses the variable as the key

model.initial_conditions = {x: pybamm.Scalar(1), y: pybamm.Scalar(2)}

Finally, we can add any variables of interest to our model. Note that these can be things other than the variables that are solved for. For example, we may want to store the variable defined by \(z=x+4y\) as a model output. Variables are added to the model using the model.variables dictionary as follows:

model.variables = {"x": x, "y": y, "z": x + 4 * y}

Note that the keys of this dictionary are strings (i.e. the names of the variables). The string names can be different from the variable symbol, and should in general be something informative. The model is now completely defined and is ready to be discretised and solved!

Using the model

We first discretise the model using the pybamm.Discretisation class. Calling the method process_model turns the model variables into a pybamm.StateVector object that can be passed to a solver. Since the model is a system of ODEs we do not need to provide a mesh or worry about any spatial dependence, so we can use the default discretisation. Details on how to provide a mesh will be covered in the following notebook.

disc = pybamm.Discretisation()  # use the default discretisation
disc.process_model(model);

Now that the model has been discretised it is ready to be solved. Here we choose the ODE solver pybamm.ScipySolver and solve, returning the solution at 20 time points in the interval \(t \in [0, 1]\)

solver = pybamm.ScipySolver()
t = np.linspace(0, 1, 20)
solution = solver.solve(model, t)

After solving, we can extract the variables from the solution. These are automatically post-processed so that the solutions can be called at any time \(t\) using interpolation. The times at which the model was solved at are stored in solution.t and the statevectors at those times are stored in solution.y

t_sol, y_sol = solution.t, solution.y  # get solution times and states
x = solution["x"]  # extract and process x from the solution
y = solution["y"]  # extract and process y from the solution

We then plot the numerical solution against the exact solution.

t_fine = np.linspace(0, t[-1], 1000)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4))
ax1.plot(t_fine, 2 * np.exp(t_fine) - np.exp(2 * t_fine), t_sol, x(t_sol), "o")
ax1.set_xlabel("t")
ax1.legend(["2*exp(t) - exp(2*t)", "x"], loc="best")

ax2.plot(t_fine, 3 * np.exp(t_fine) - np.exp(2 * t_fine), t_sol, y(t_sol), "o")
ax2.set_xlabel("t")
ax2.legend(["3*exp(t) - exp(2*t)", "y"], loc="best")

plt.tight_layout()
plt.show()

In the next notebook we show how to create, discretise and solve a PDE model in pybamm.

References

The relevant papers for this notebook are:

pybamm.print_citations()

[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.
[2] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.
[3] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). ECSarXiv. February, 2020. doi:10.1149/osf.io/67ckj.
[4] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, and others. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17(3):261–272, 2020. doi:10.1038/s41592-019-0686-2.