Working principles

This page provides an introduction to FONSim’s internal workings. Its aim is to provide the reader with an overview such that they can decide whether FONSim is the correct tool for their use case as well as reason about how they can use it the most effectively. Readers are also welcome to read the FONSim paper.

This document deliberately does not go into specifics or practicalities. For practical how-to guides on how to use FONSim, please have a look at the tutorials. There is also the formal code documentation which discusses all important functions and classes in FONSim. This page assumes that you have read the first paragraph of the introduction.

Variables and terminals

A system in FONSim is defined as a bidirectional graph of components connected by their terminals and each terminal has one or more variables. Two terminals are considered to be connected when they are connected to a common node.

A variable corresponds with a scalar time dependent value, for example a mass flow or a pressure. Following classical systems modeling, a variable in a terminal is considered to be over or through. Variables are also used to denote a component state (for example, the amount of fluid inside).

A terminal is a set of variables that together sufficiently describe the interaction trough that terminal. For example, for isothermal fluidics a possible choice is a pressure (over variable) and a mass flow (through variable). As another example, mechanical rotational interactions (such as a drive shaft) can be described by torque (over variable) and rotational velocity (through variable).

Components and their equations

A component has one or more terminals as well as two functions. The terminals define the component’s possibilities to interact with other components, while the functions describe its intrinsic behavior. These two functions are the update_state function (unless the component is stateless) and the evaluate function. They are expressed in plain Python and can therefore contain any desired Python expression.

Both functions take as arguments the variable values. The update_state function takes as additional arguments the current state and the timestep size and it returns a new state. (Any continuous differential equations must therefore be discretized to a difference equation first.) The evaluate function returns the equation residuals that should equal zero at each simulation step.

Derivatives can be defined explicitly, or can be left out in which case FONSim will estimate them using finite differences. A more practical introduction to components that goes deeper into the usage specifics is provided by this tutorial.

Simulation and solving

Once a system has been constructed, FONSim’s simulation module allows to derive numerical results from it. The default FONSim solver is a fairly straightforward implementation of Newton-Raphson with backward Euler discretization for improved stability. (Given that this backward discretization is not apparent in the component definition, it is relatively easy to employ forward Euler instead.) In addition, the default solver is equipped with automatic timestep adjustment. In usage, simulating a defined system takes minimal effort (see, for example, this tutorial) because FONSim makes abstraction of the simulation procedure.


When asked to simulate a system, FONSim initializes the candidate solution vector \(\mathbf{x_a}\), which contains the candidate values of all variables in terminals, and the square system matrix \(\mathbf{J}\). Hashtables (Python dictionaries) map variables in the object oriented system structure to matrix indices in \(O(1)\) time. FONSim also builds using Kirchoff laws the network equations that describe the relations between variables of different components. These equations are linear and constant for a given system.

The total residual vector \(\mathbf{f}\) takes as argument the candidate solution vector \(\mathbf{x}\) and returns the residual. If the solution were to be exactly correct, the residual would be a vector of zeros. The system matrix \(\mathbf{J}\) is the derivative of this residuals vector \(\mathbf{f}\) to the solution vector \(\mathbf{x}\) (i.e. Jacobian matrix).

The upper part \(\mathbf{J_n}\) of the system matrix expresses the network equations and is therefore constant over a simulation. On the other hand, the lower part \(\mathbf{J_a}\) expresses the derivatives of the component equations. As those equations tend to be nonlinear, \(\mathbf{J_a}\) has to be updated in every iteration of the Newton outer loop.


After these preparations, FONSim starts iterating over the timesteps, solving timestep by timestep. At each timestep FONSim solves the nonlinear residual function \(\mathbf{f}\) using the Newton-Raphson method. Each Newton outer loop iteration starts with evaluating this residual vector with the candidate solution \(\mathbf{x}\) For the component residuals, this is done by iterating over all components and evaluating their update_state and evaluate functions. (FONSim does not do any symbolic math on the given expressions.) The network residuals follow from multiplying \(\mathbf{J_n}\) with \(\mathbf{x}\). Subsequently, this residual is backprojected onto \(\mathbf{x}\) using a highly optimized Numpy linear system solver. This outer loop is repeated until all residuals are lower than set thresholds. (These thresholds depend on the variables’ ranges to avoid scale effects.)

While evaluating the component residuals of each component, FONSim also gathers the derivatives of the evaluate and the update_state functions. The product of these two is then written to the appropriate location in \(\mathbf{J_a}\), Had forward Euler been used, only the derivatives of the evaluate function would be used.

The timestep adjustment is handled using an interpolation criterion. Given \(\mathbf{x}\) solved at three subsequent simulation steps \(k-2\), \(k-1\) and \(k\), the vector \(\mathbf{x}^{k-1}\) must be within the desired precision from the \((\mathbf{x}^{k-2} + \mathbf{x}^{k})/2\). If this is not the case, the timestep is reduced, while if it is far closer than the desired precision, the timestep is increased. This desired precision can be specified as relative to the variable’s range. Each time the timestep is reduced, the simulator returns two steps to ensure proper sampling.

Once the simulation has finished, the simulation results are easily accessed for visualization or other purposes. One can access them via the components because their internal variables are mapped to the simulation solution matrix.


This page explains the working principles of FONSim. The behavior of the interconnected components is defined using variables, terminals and equations. Together the connected components form a system, which is modeled as a bidirectional graph. Components have terminals that define how they can interact with other components and functions that describe their intrinsic behavior, which can be expressed in plain Python. FONSim’s default solver is a straightforward implementation of the Newton-Raphson method with backward Euler discretization, and it has automatic timestep adjustment. The system matrix expresses the network equations and the derivatives of component equations. Simulating a defined system takes minimal effort because FONSim makes abstraction of the simulation procedure.

For further reading, it is suggested to have a look at the tutorials and the formal code documentation.