README

SOLVER Toolset - Free Pascal Console Mode Port

A port of the STAMS SOLVER Toolset v5.04 (W.S.C. Gurney, University of Strathclyde, 1995) from Sun Pascal / X11 to Free Pascal 3.2+ with console-mode I/O.

SOLVER is a template-based system for generating applications that solve sets of coupled first-order ordinary or delay differential equations. It was originally developed to accompany the book Ecological Dynamics by W.S.C. Gurney and R.M. Nisbet (Oxford University Press, 1998).

Prerequisites

• Free Pascal Compiler (fpc) 3.2.2 or later

Building

Compile units and utilities

make all

This compiles the four Pascal units (turbop, mxfiler, matrixunit, consoleui) and the utility programs (mxtx, txmx, joinmx, plotter, mkpg, mx2csv).

Compile and run a solver problem

cd examples/solver
../../solv.sh exp

See Using solv.sh below for details.

Overview

SOLVER solves initial value problems of the form:

dx_i/dt = f_i({x_j}, {c_k}, t)     i = 1..N

where x_i are the state variables, c_k are constant coefficients, and t is time. It also supports delay differential equations where the rates of change depend on historical values of the state variables.

The user defines a problem by writing two small Pascal files:

• <probname>.con - a constants file specifying array dimensions
• <probname>.def - a definitions file containing four procedures

The build script solv.sh combines these with the solver template, compiles, and (optionally) launches the resulting program.

Problem Definition

The Constants File (<probname>.con)

The constants file declares the dimensions of the problem:

CONST
    NoofStaVar = 1;     { Number of state variables }
    NoofCoeffs = 2;     { Number of coefficients }
    NoofPltVar = 1;     { Number of variables plotted on screen }
    NoofOptVar = 1;     { Number of variables output to file }
    NoofHisVar = 1;     { For delay problems (set to 1 if not used) }
    RingBufLen = 0;     { Ring buffer length (0 = no delay) }

Based on these values, SOLVER implicitly defines the types:

SVEC = ARRAY[1..NoofStaVar] OF DOUBLE;     { State vector }
CVEC = ARRAY[1..NoofCoeffs] OF DOUBLE;     { Coefficient vector }
HVEC = ARRAY[1..NoofHisVar] OF DOUBLE;     { History vector }
MPAD = RECORD                              { Memo pad }
         opt : ARRAY[1..NoofOptVar] OF DOUBLE;
         plt : ARRAY[1..NoofPltVar] OF DOUBLE;
         his : HVEC;
       END;

The Definitions File (<probname>.def)

The definitions file must define four procedures:

Get_Gradient

The main problem specification. Calculates the rates of change for the current state:

PROCEDURE Get_Gradient(VAR CurGrad: SVEC; CurState: SVEC;
                       Coeffs: CVEC; time: DOUBLE);
BEGIN
    CurGrad[1] := Coeffs[1] * CurState[1] * (1.0 - CurState[1] / Coeffs[2]);
END;

Get_Memo

Defines which variables are plotted on screen (plt[]), output to file (opt[]), and stored in the history buffer (his[]):

PROCEDURE Get_Memo(VAR memo: MPAD; CurState: SVEC;
                   Coeffs: CVEC; time: DOUBLE);
BEGIN
    WITH memo DO BEGIN
        opt[1] := CurState[1];
        plt[1] := CurState[1];
        his[1] := 0.0;
    END;
END;

Get_Initial_State

Optional automatic initialisation. Called before prompting the user for initial values. Any state variable set here will not be prompted for:

PROCEDURE Get_Initial_State(VAR Initial_State: SVEC; Coeffs: CVEC);
VAR i: INTEGER;
BEGIN
    FOR i := 1 TO NoofStaVar DO Initial_State[i] := 0.0;
END;

Leave as an empty stub if you want to enter all initial values interactively.

Get_History

For delay problems only. Returns the value of the history vector at times before t=0:

PROCEDURE Get_History(VAR History: HVEC; Coeffs: CVEC; time: DOUBLE);
BEGIN
    History[1] := 0.0;   { e.g., no recruitment before t=0 }
END;

Leave as an empty stub for non-delay problems.

Coding Techniques

For readability, you can declare named constants to index the state and coefficient vectors:

CONST
    N = 1;  P = 2;          { state variable indices }
    Lambda = 1;  Q = 2;     { coefficient indices }

PROCEDURE Get_Gradient(VAR CurGrad: SVEC; CurState: SVEC;
                       Coeffs: CVEC; time: DOUBLE);
BEGIN
    CurGrad[N] := Coeffs[Lambda] * CurState[N] * (1.0 - CurState[N]);
    CurGrad[P] := CurState[N] * (1.0 - exp(-Coeffs[Q] * CurState[P]));
END;

You can also use local variables for clarity:

PROCEDURE Get_Gradient(VAR CurGrad: SVEC; CurState: SVEC;
                       Coeffs: CVEC; time: DOUBLE);
VAR R, K, pop, dNdt: DOUBLE;
BEGIN
    R := Coeffs[1];  K := Coeffs[2];  pop := CurState[1];
    dNdt := R * pop * (1.0 - pop / K);
    CurGrad[1] := dNdt;
END;

Using solv.sh

The solv.sh script compiles a solver problem from .con and .def files, producing an executable. Run it from the directory containing your problem files:

../../solv.sh [switches] <probname>

Command-Line Switches

What solv.sh Does

1. Copies <probname>.con to solvcons.inc
2. Copies <probname>.def to solvdefs.inc
3. For SWV (-s): copies <probname>.swt to solvswch.inc
4. Copies the selected integrator to solvint.inc
5. Compiles the solver template with fpc
6. Produces executable <probname>.solver
7. Launches the executable

Compile Only (No Launch)

To compile without launching, run fpc directly:

SCRIPTDIR="../../"
cp myproblem.con "$SCRIPTDIR/toolsrc/solvcons.inc"
cp myproblem.def "$SCRIPTDIR/toolsrc/solvdefs.inc"
cp "$SCRIPTDIR/toolsrc/solvint_rk4.inc" "$SCRIPTDIR/toolsrc/solvint.inc"
fpc -Fu"$SCRIPTDIR/units" -Fi"$SCRIPTDIR/toolsrc" -FE. -o"myproblem.solver" "$SCRIPTDIR/toolsrc/solver_fsv.pas"
rm -f "$SCRIPTDIR/toolsrc/solvcons.inc" "$SCRIPTDIR/toolsrc/solvdefs.inc" "$SCRIPTDIR/toolsrc/solvint.inc"

Running SOLVER

When a solver program is launched, it follows an interactive console session:

1. File Setup

You are prompted for:

• Results file - name of an MxFile (.mx) for output data. Press Enter for no output file.
• Parameters file - name of a text file for parameter logging. Press Enter to skip.

2. Calculation Parameters

For the fixed timestep solver (-f, default):

• Time step (dt) - the integration time step
• Start time - when integration begins (default 0)
• Output increment - interval between file outputs (0 = every step)
• Halt time - when integration stops

For the adaptive timestep solver (-a):

• Start time - when integration begins (default 0)
• Stop time - when integration stops
• Minimum timestep - the smallest allowed timestep
• Maximum timestep - the largest allowed timestep (expressed as a multiple of the minimum)
• Output interval - time interval between file outputs (0 = output every step)
• Tolerable error - the maximum acceptable relative error per step (controls accuracy)

3. Coefficients

You are prompted for each coefficient C[1], C[2], ..., C[NoofCoeffs].

4. Initial State

Unless Get_Initial_State has set them, you are prompted for each state variable S[1], S[2], ..., S[NoofStaVar].

5. Ring Buffer (Delay Problems)

If RingBufLen > 0, you are asked to confirm the ring buffer length.

6. Display Mode

Choose the display mode:

• T (TimePlot) - plots plt[1]..plt[NoofPltVar] against time (default)
• P (PhasePlot) - plots even plt[] components against odd components
• V (VectorPlot) - plots state vector components against their index

7. Ordinate Scale

Set the maximum and minimum values for the Y-axis display.

8. Integration

The integration runs until the halt time is reached or a key is pressed. After halting, you can:

• Continue to a new endpoint
• Start a new run (parameters default to previous values)
• Quit

YAML Batch Mode

Instead of entering parameters interactively, you can supply them in a YAML file:

./exp.solver -y exp.yaml

The YAML file specifies all parameters that would otherwise be entered interactively:

# exp.yaml - Exponential decay example
results: exp_out.mx
params: ""                # optional, "" for no param file
dt: 0.1                   # integration time step
start_time: 0.0
output_interval: 0.1
halt_time: 10.0
coefficients:
  - 0.5                   # C[1] = decay rate
initial_state:
  - 1.0                   # S[1] = initial value
display_mode: T            # T = timeplot, P = phaseplot, V = vectorplot
ordinate_min: 0.0
ordinate_max: 1.2

YAML Parameters

In YAML mode, the solver runs a single calculation and exits. The initial_state list provides values for state variables not already set by Get_Initial_State; template-initialised variables are not overridden.

File Output Format

Each row of the output MxFile contains:

Integration Methods

Runge-Kutta 4th Order (-rk4, default)

A conservative, generally reliable algorithm that allows reasonably large time steps. Recommended for most problems.

Runge-Kutta 2nd Order (-rk2)

Sometimes preferable for problems with numerous slope discontinuities, where smaller time steps compensate for the lack of high-order continuity without the computational overhead of RK4.

Euler's Method (-eul)

Should not be used for integrating ODEs as it is notoriously unstable and inaccurate. It is useful for:

• Integrating stochastic differential equations (the integral is rigorously interpretable as an Ito integral)
• Producing time series from iterated maps (use fixed timesteps with dt=1)

Delay Differential Equations

SOLVER supports problems where the rate of change depends on historical values of the state variables. To use delay features:

1. Set the Ring Buffer

In the .con file, set RingBufLen to a non-zero value (typically 1000). This determines how many past time steps are stored. The buffer covers a time span of RingBufLen * dt.

2. Define History Variables

Set NoofHisVar to the number of quantities you want to track historically. In Get_Memo, assign values to memo.his[1]..memo.his[NoofHisVar].

3. Initialise the History Buffer

In Get_History, return the value of the history vector for times before t=0. This procedure is called with values of time in the range (-dt * RingBufLen)..0.

4. Access Past Values

Use the built-in function:

FUNCTION PastValue(N: HisVarNo; time: DOUBLE): DOUBLE;

This returns the value of history variable N at the specified time. If the requested time does not exactly match a calculation point, quadratic interpolation is used. The time must lie within the range covered by the ring buffer.

Example: Nicholson's Blowflies

{ blowfly.con }
CONST
    NoofStaVar = 1;  NoofCoeffs = 4;
    NoofOptVar = 1;  NoofPltVar = 1;
    NoofHisVar = 1;  RingBufLen = 1000;

{ blowfly.def }
PROCEDURE Get_History(VAR History: HVEC; Coeffs: CVEC; time: DOUBLE);
BEGIN  History[1] := 0.0;  END;

PROCEDURE Get_Initial_State(VAR Initial_State: SVEC; Coeffs: CVEC);
BEGIN  { do nothing - prompt user }  END;

PROCEDURE Get_Gradient(VAR CurGrad: SVEC; CurState: SVEC;
                       Coeffs: CVEC; time: DOUBLE);
VAR N, AdultRecruitment, Delta, tau: DOUBLE;
BEGIN
    Delta := Coeffs[3];  tau := Coeffs[4];  N := CurState[1];
    AdultRecruitment := PastValue(1, time - tau);
    CurGrad[1] := AdultRecruitment - Delta * N;
END;

PROCEDURE Get_Memo(VAR memo: MPAD; CurState: SVEC;
                   Coeffs: CVEC; time: DOUBLE);
VAR EggProd, P, N0, N: DOUBLE;
BEGIN
    P := Coeffs[1];  N0 := Coeffs[2];  N := CurState[1];
    EggProd := P * N * exp(-N / N0);
    WITH memo DO BEGIN
        plt[1] := N;   opt[1] := N;   his[1] := EggProd;
    END;
END;

Here, coefficients are: C[1]=P (max egg production rate), C[2]=N0 (population scale), C[3]=delta (death rate), C[4]=tau (development delay).

Switching Function Solver (SWV)

The switching function variant (-s) extends the fixed timestep solver to handle problems with state-dependent discontinuities. A "switch" is a condition that turns on or off when a user-defined function of the state crosses zero. When a switch transition is detected during a time step, the solver uses quadratic interpolation to locate the exact transition point, steps to it, calls user-defined reset code, and then continues.

The Switch File (<probname>.swt)

In addition to .con and .def, SWV problems require a .swt file defining two procedures:

Get_Switch_Functions

Evaluates the switch functions for the current state. Each switch function's sign determines the switch state: non-negative means ON (True), negative means OFF (False).

PROCEDURE Get_Switch_Functions(VAR SwFuncVal: SWFVEC; State: SVEC;
                                                      t: DOUBLE; C: CVEC);
BEGIN
    SwFuncVal[1] := threshold - PastValue(1, t - delay);
END;

Reset_State

Called at each transition point. Allows the state vector to be modified when a switch changes. HiToLo is True when the switch is going from ON to OFF.

PROCEDURE Reset_State(SwNo: INTEGER; TransTime: DOUBLE; HiToLo: BOOLEAN;
                                                    VAR S: SVEC; C: CVEC);
BEGIN
    { modify S as needed at the transition }
END;

Constants File for SWV

The .con file must declare NoofSwitches in addition to the standard constants:

CONST
    NoofStaVar   = 1;
    NoofHisVar   = 1;
    NoofOptVar   = 1;
    NoofPltVar   = 1;
    NoofCoeffs   = 1;
    RingBufLen   = 500;
    NoofSwitches = 1;

Using Switch_State in Definitions

In the .def file, use the built-in function Switch_State(n) to query whether switch n is currently ON (True) or OFF (False):

PROCEDURE Get_Gradient(VAR Gradient: SVEC; State: SVEC; Coeffs: CVEC; t: DOUBLE);
BEGIN
    IF Switch_State(1) THEN
        Gradient[1] := 1.0 - Coeffs[1] * State[1]
    ELSE
        Gradient[1] := -Coeffs[1] * State[1];
END;

Switch Report File

When running interactively, SWV prompts for a switch report file (.rep). This logs every switch transition with the switch number, new state, and exact time. In YAML batch mode, use the switch_report key:

switch_report: kelp.rep

Set to empty string or omit to skip the report file.

Partial Differential Equations

SOLVER can be used for PDEs by discretising the spatial dimension. Identify each grid point with a state variable, and compute spatial derivatives using finite differences in Get_Gradient. Use Get_Initial_State to set up the initial spatial distribution. VectorPlot display mode is useful for visualising the distribution directly.

Utility Programs

mxtx - MxFile to Text

Convert binary MxFiles to human-readable text format:

./utils/mxtx result

Reads result.mx and outputs result.tx. Supports multiple files and precision control:

./utils/mxtx -p10 file1 file2 file3

The -p switch sets the number of significant figures (default: 6).

mx2csv - MxFile to CSV/TSV

Convert binary MxFiles to standard CSV or TSV format, suitable for import into R, Python, Excel, and other tools:

./utils/mx2csv result.mx

Reads result.mx and outputs result.csv. Options:

Examples:

./utils/mx2csv result.mx                     # result.mx -> result.csv
./utils/mx2csv -d tab -h result.mx           # result.mx -> result.tsv with header
./utils/mx2csv -p 10 -o output.csv result.mx # 10 sig figs, custom output name
./utils/mx2csv file1.mx file2.mx file3.mx    # batch conversion

txmx - Text to MxFile

Convert text data files to MxFile binary format:

./utils/txmx mydata

Reads mydata (or mydata.tx) and outputs mydata.mx. Uppercase filenames are converted to lowercase for the output file.

joinmx - Join MxFiles

Join two MxFiles side by side into a single file:

./utils/joinmx file1 file2

Reads file1.mx and file2.mx, producing union.mx with columns from both files. Options:

./utils/joinmx file1 file2 -o output    # Output to output.mx instead of union.mx
./utils/joinmx -d5 file1 file2          # Align first row of file2 at row 5

If the files have different numbers of rows, empty slots are filled with NoEntry markers.

plotter - Data to PostScript/PDF Plot

Generate publication-quality PostScript or PDF plots from MxFile, CSV, TSV, or whitespace-delimited text data. The input format is auto-detected by file extension:

PDF output is generated natively using fcl-pdf -- no external tools required.

./utils/plotter -x 1 -y 2 -t "Exponential Decay" -xl "time" -yl "S(1)" -o exp.ps result.mx
./utils/plotter -x 1 -y 2,3,4 -l "S,I,R" -t "SIR Model" -pdf -o sir.pdf sir_out.csv

Options:

By default, plotter produces self-contained .ps files viewable in any PostScript viewer (e.g. open, gs, Preview.app). Use -f for fragment mode when assembling multiple plots with mkpg.

PDF output is available via the -pdf flag or by using a .pdf extension with -o. PDF files are generated natively using the fcl-pdf library -- no external tools are required.

Text input notes: For CSV and TSV files, if the first field of the first line is non-numeric (e.g. a column header), the entire first line is treated as a header and skipped. Empty fields are treated as missing data (NoEntry), so gaps in the data are handled the same way as in MxFile input.

Multiple curves example:

./utils/plotter -y 2,3 -t "Two Curves" -o multi.ps result.mx

Marks without lines:

./utils/plotter -m triangle -n -o marks.ps result.mx

mkpg - Page Assembly

Assemble multiple PostScript plot fragments into a single page:

./utils/mkpg ./utils/PsPlot.eps plot1.ps plot2.ps plot3.ps

This creates page.epsf with the plots automatically arranged. Options:

Plotting Workflow

The typical workflow for generating publication-quality plots:

solver --> result.mx --> plotter --> plot.ps/plot.pdf --> viewer
                  |              |
                  |              +--> plotter -f --> frag.ps --> mkpg --> page.epsf
                  |
                  +--> mx2csv --> result.csv --> R / Python / Excel
                                     |
                                     +--> plotter --> plot.pdf

Single plot (self-contained):

cd examples/solver && ../../solv.sh exp
# (run solver, save results to result.mx)
../../utils/plotter -x 1 -y 2 -t "Exponential Decay" -xl "time" -yl "S(1)" -o exp.ps result.mx
open exp.ps

Multi-plot page assembly:

../../utils/plotter -f -x 1 -y 2 -o plot1.ps result1.mx
../../utils/plotter -f -x 1 -y 2 -o plot2.ps result2.mx
../../utils/mkpg ../../utils/PsPlot.eps plot1.ps plot2.ps
open page.epsf

MxFiles

MxFiles are the binary data format used by SOLVER for all file output. Each MxFile stores a two-dimensional matrix of double-precision floating-point numbers (8 bytes per element). Key properties:

• Elements are either a meaningful number or a special NoEntry sentinel value
• The file carries metadata indicating the number of rows and columns
• Elements can be accessed directly by row and column index
• File names use the .mx extension
• All data output from SOLVER is automatically stored in MxFile format

To inspect the contents of an MxFile, use mxtx to convert it to text.

Available Examples

solver/

Basic examples from the original distribution:

models/

Models from the textbook:

advmods/

Advanced models:

ecodyn/

Examples organised by textbook chapter (chapters 2-9), from Ecological Dynamics by Gurney and Nisbet.

Directory Structure

solver-fpc/
  Makefile              Build units and utilities
  solv.sh               Compile-and-run script for solver problems
  units/
    turbop.pas          String/number conversion utilities
    mxfiler.pas         MxFile binary I/O
    matrixunit.pas      Linear algebra (eigenvalues, determinants)
    consoleui.pas       Console UI (replaces X11 WinLib/PlotLib)
    psplot.pas          PostScript plotting extension
    pdfplot.pas         Native PDF plotting (fcl-pdf)
    yamlinput.pas       YAML input file reader
    otyaml/             otYaml YAML parser library (bundled)
  toolsrc/
    solver_fsv.pas      Fixed time-step solver template
    solver_asv.pas      Adaptive time-step solver template
    solver_swv.pas      Switching function solver template
    solvcom.inc         Common solver procedures
    ringbuf.inc         Ring buffer for delay DEs
    solvint_rk4.inc     Runge-Kutta 4th order integrator
    solvint_rk2.inc     Runge-Kutta 2nd order integrator
    solvint_eul.inc     Euler integrator
  utils/
    mxtx.pas            MxFile to text converter
    mx2csv.pas          MxFile to CSV/TSV converter
    txmx.pas            Text to MxFile converter
    joinmx.pas          MxFile joiner
    plotter.pas         Plot generator (MxFile, CSV, TSV, text input)
    mkpg.pas            PostScript page assembly
    PsPlot.eps          PostScript prolog (used by mkpg)
    psprn.sh            Convenience print script
  examples/
    solver/             Basic examples (exp, nich)
    models/             Textbook models (12 problems)
    advmods/            Advanced models (12 problems)
    ecodyn/             Ecological Dynamics examples (chapters 2-9)

Differences from the Original

This port makes the following changes from the original Sun Pascal / X11 version:

• Compiler: Sun Pascal / DEC Pascal replaced by Free Pascal (fpc)
• GUI: X11 windowing (WinLib/PlotLib) replaced by console ReadLn/WriteLn
• Graphics: On-screen plotting is replaced by text messages; results are written to MxFiles which can be analysed with mxtx or external tools
• Build system: The solv shell script is replaced by solv.sh with the same command-line interface
• Module system: Sun Pascal module declarations replaced by Free Pascal unit declarations
• Preprocessor: C preprocessor directives (#include, #ifdef) replaced by Free Pascal equivalents ({$I}, {$IFDEF})
• Solver variants: The fixed timestep solver (-f), adaptive timestep solver (-a), and switching function solver (-s) are all ported
• Platform: Originally Solaris/SunOS/OSF1; this port targets macOS (aarch64) but should work on any platform supported by Free Pascal

References

• Gurney, W.S.C. and Nisbet, R.M. (1998). Ecological Dynamics. Oxford University Press.
• Nisbet, R.M. and Gurney, W.S.C. (1982). Modelling Fluctuating Populations. John Wiley & Sons.
• Original SOLVER toolset: STAMS, University of Strathclyde, Glasgow G1 1XH, Scotland.

Third-Party Code

• The YAML parser in units/otyaml/ is from otYaml, an event-based YAML parser for Free Pascal.