odesolve.c File Reference

Adaptive integration of ordinary differential equations. More...

#include <stdlib.h>
#include <math.h>
#include "odesolve.h"

Include dependency graph for odesolve.c:

Macros
#define	MAXSTP   10000
#define	TINY   1.0e-30
#define	SAFETY   0.9
#define	PGROW   -0.2
#define	PSHRNK   -0.25
#define	ERRCON   1.89e-4

Functions
int	rkqs (double *x, int n, double htry, double eps, double *hdid, double *hnext, void(*derivs)(double, double *, double *))
	Reject and shrink trial RK45 steps until scaled local error is acceptable.
void	rkck (double x, int n, double h, void(*derivs)(double, double *, double *))
	Compute one Cash-Karp trial state and embedded truncation error.
int	odesolve_open (int n)
	Allocate the legacy solver's shared workspace for n equations.
void	odesolve_close ()
	Free every legacy ODE workspace array and reset capacity to zero.
int	odesolve_integrate (double ystart[], int n, double x1, double x2, double eps, double h1, void(*derivs)(double, double *, double *))
	Integrate an ODE system with adaptive Cash-Karp RK45 steps.

Variables
int	nmax
double *	y
double *	yscal
double *	yerr
double *	ytemp
double *	dydx
double *	ak

Detailed Description

Adaptive integration of ordinary differential equations.

Author

Caleb Buahin caleb.nosp@m..bua.nosp@m.hin@g.nosp@m.mail.nosp@m..com (OpenSWMM maintenance and extensions)

Copyright

Copyright (c) 2026 Caleb Buahin. OpenSWMM contributions only; original EPA SWMM credits are retained.

License

Apache-2.0 for OpenSWMM contributions; original EPA SWMM provenance is retained.

Macro Definition Documentation

ERRCON

#define ERRCON   1.89e-4

MAXSTP

#define MAXSTP   10000

PGROW

#define PGROW   -0.2

PSHRNK

#define PSHRNK   -0.25

SAFETY

#define SAFETY   0.9

TINY

#define TINY   1.0e-30

Function Documentation

odesolve_close()

void odesolve_close

(

void

)

Free every legacy ODE workspace array and reset capacity to zero.

Note

Safe after partial allocation or a previous close; pointers are reset to NULL. Must not run concurrently with integration or inside its callback.

Here is the caller graph for this function:

odesolve_integrate()

int odesolve_integrate	(	double	ystart[],
		int	n,
		double	x1,
		double	x2,
		double	eps,
		double	h1,
		void(*	derivs )(double, double *, double *) )

Integrate an ODE system with adaptive Cash-Karp RK45 steps.

Parameters

[in,out]	ystart	Non-null array of n initial states at x1, replaced with states at x2 on success; unchanged on an error return.
[in]	n	Positive equation count, no larger than allocated nmax.
[in]	x1	Finite start coordinate in caller-selected units.
[in]	x2	Finite end coordinate in the same units.
[in]	eps	Positive dimensionless tolerance relative to each state scale abs(y) + abs(h * dydx) + TINY.
[in]	h1	Nonzero initial step, in x units, directed from x1 toward x2.
[in]	derivs	Non-null synchronous callback filling n dy/dx values from n states; borrowed arrays must not be modified or retained except for writing the derivative output.

Returns

Zero on success, one when n exceeds nmax, two for a zero next step, or three for exhausted iterations or underflow during rejection.

Precondition

odesolve_open() succeeded. The callback must not recursively integrate because accepted states and intermediate buffers are shared globals.

Here is the call graph for this function:

Here is the caller graph for this function:

odesolve_open()

int odesolve_open

(

int

n

)

Allocate the legacy solver's shared workspace for n equations.

Parameters

[in]

n

Positive equation capacity; n * 5 must be representable as int.

Returns

One after all allocations succeed, zero when an allocation returns NULL.

Postcondition

Sets nmax only after success; partial allocations remain for odesolve_close(). Arrays are zero-initialized and shared across all calls.

Precondition

Any earlier workspace has been closed. Calling open twice without close overwrites owned pointers. The legacy solver is not thread-safe.

Here is the caller graph for this function:

rkck()

void rkck	(	double	x,
		int	n,
		double	h,
		void(*	derivs )(double, double *, double *) )

Compute one Cash-Karp trial state and embedded truncation error.

Parameters

[in]	x	Current independent-variable coordinate.
[in]	n	Positive equation count within allocated shared workspace.
[in]	h	Trial step in x units; may be negative.
[in]	derivs	Non-null synchronous callback filling n derivatives from the supplied n-state arrays; no temporary pointer may be retained.

Precondition

y[] and dydx[] contain the accepted state and its current derivative.

Postcondition

Fills ytemp[] with the fifth-order trial state and yerr[] with its embedded error estimate. Does not commit the shared accepted y[] state.

Here is the caller graph for this function:

rkqs()

int rkqs	(	double *	x,
		int	n,
		double	htry,
		double	eps,
		double *	hdid,
		double *	hnext,
		void(*	derivs )(double, double *, double *) )

Reject and shrink trial RK45 steps until scaled local error is acceptable.

Parameters

[in,out]	x	Non-null current coordinate; advanced only on acceptance.
[in]	n	Positive equation count covered by shared workspace.
[in]	htry	Initial trial step in independent-variable units.
[in]	eps	Positive dimensionless scaled-error tolerance.
[out]	hdid	Non-null slot for the accepted step; unchanged on underflow.
[out]	hnext	Non-null slot for the next step proposal; unchanged on underflow.
[in]	derivs	Non-null synchronous derivative callback; temporary state and output arrays are borrowed and must not be retained.

Returns

Zero after accepting shared y[] and updating x, hdid and hnext; two when a rejected step becomes too small to change the coordinate.

Precondition

Current y, dydx and yscal arrays are initialized by odesolve_integrate().

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

ak

double* ak

dydx

double* dydx

nmax

int nmax

y

double* y

yerr

double* yerr

yscal

double* yscal

ytemp

double* ytemp