The Freeode project provides a simple programming language, to numerically solve differential equations. The compiler produces source code in the interpreted language Python, that performs the numerical computations. Both dynamic and steady state simulations can be done. The SciPy and NumPy libraries are used for the numerical computations. Currently only ODE can be solved.
The generated Python code can be used as a stand-alone program, or in an interactive Python session. Especially can the generated Python objects be used as building blocks of more complex programs. The programmer is freed from the relatively mechanical task of implementing the simulation function, and can concentrate on the higher level aspects of the problem.
Python promises to be highly platform independent. Therefore everything should work on Windows exactly as it works on Linux.
The free SciPy and NumPy libraries turn the Python language into a prototyping environment for numerical algorithms, quite similar to the commercial language Matlab. There are several other free numerical prototyping languages; examples are Octave and Scilab. Python however is a general-purpose programming language with a rich standard library. Most notably here, is Python's powerfull string processing and file I/O.
The compiler is licensed under the GPL. The runtime libraries are licensed under the LGPL so they can be linked to commercial applications. The generated program can off course be licensed under any license you wish.
This rough example should give a first impression of the usage of the compiler and of the generated program. Also it should show what code the compiler emits.
A simple biological reactor should be simulated. The simulation has two state variables: X the concentration of biomass, and S the concentration of sugar. The growth speed µ is an algebraic variable.
dX/dt = µ*X | Change of biomass concentration |
dS/dt = -1/Yxs*µ*X | Change of sugar concentration |
with: | |
µ = µmax * S/(S+Ks) | Growth speed of biomass |
Initial values and the values of the parameters have been omitted for brevity.
#Biological reactor with no inflow or outfow class Batch: #Define values that stay constant during the simulation. data mu_max, Ks, Yxs: Float param #Define values that change during the simulation. data mu, X, S: Float #Initialize the simulation. func initialize(this): #Specify options for the simulation algorithm. solution_parameters(duration=20, reporting_interval=0.1) #Give values to the parameters mu_max = 0.32; #max growth speed Ks = 0.01; #at this sugar concentration growth speed is 0.5*mu_max Yxs = 0.5; #one g sugar gives this much biomass #Give initial values to the state variables. X = 0.1; #initial biomass concentration S = 20; #initial sugar concentration #compute dynamic behaviour - the system's 'equations' func dynamic(this): mu = mu_max * S/(S+Ks); #growth speed (of biomass) $X = mu*X; #change of biomass concentration $S = -1/Yxs*mu*X; #change of sugar concentration #show results func final(this): graph(mu, X, S); compile Batch
This is a complete SIML program to to solve the system of differential equations. The differential equations are in the dynamic function. The init function is invoked once at the beginning of the simulation, the final function is invoked at the end.
The source distribution contains a syntax coloring file, for editors based on the Katepart. Therfore the editors Kate and Kwrite (and Kdevelop) can show SIML code colored like the example above. (The HTML of both code examples was generated by the Kate editor's HTML export feature.)
def dynamic(self, time, state, returnAlgVars=False): ''' Compute time derivative of state variables. This function will be called by the solver repeatedly. ''' #take the state variables out of the state vector v_S = state[0] v_X = state[1] #do computations v_mu = self.p_mu_max*v_S/(v_S + self.p_Ks) v_X_dt = v_mu*v_X v_S_dt = -1.0/self.p_Yxs*v_mu*v_X if returnAlgVars: #put algebraic variables into array algVars = array([v_mu, ], 'float64') return algVars else: #assemble the time derivatives into the return vector stateDt = array([v_S_dt, v_X_dt, ], 'float64') return stateDt
This is the generated simulation function. The compiler emits code with comments, so that humans have a chance to understand it. Much of the generated Python code has been omitted. The omitted parts initialize the simulation object, and show the results.
The compiler does what a human programmer would also (have to) do. The problem dependent part is only three lines (below the comment: do computations). The rest of the code takes numbers out of arrays, or puts them into arrays. When the length of the state vector changes, a human programmer would have to update the array indices in the whole program. A tedious and error prone work.
The simulation function can return either derivatives of the state variables, or algebraic variables. The algebraic variables have to be computed separately, because the solver algorithm does only see the derivatives.
These are the (bash) commands to edit the program, compile it, and run it.
$> kwrite bioreactor_simple.siml & #Edit Siml file $> simlc bioreactor_simple.siml #Run compiler $> ./bioreactor_simple.py #Run generated program
The compiler can also run the generated Program. This is useful for the development of simulation programs.
$> kwrite bioreactor_simple.siml & #Edit Siml file $> simlc bioreactor_simple.siml -r all #Run compiler
After the commands have been executed, a window opens, that shows the simulation results: