Numerical Tools

Version 2.0 (Apr 6, 2022)

The Numerical Tools Unity asset is a collection of numerical methods functions that have been adapted from license-permissive sources. It currently focuses on solving ordinary differential equations and root-finding methods. In general the code is ported from either scientific python (SciPy) or methods from the book “Numerical Methods in Physics with Python” by Alex Gezerlis (NumPhysPy).

The principle methods provided are:

SciPy’s solve_ivp: Solve initial value problems using Runge-Kutta 23, RK45 or DOP853. Adaptations are provided for real-time evolution. Events and dense output are supported.
SciPy’s solve_bvp: Solve an ODE boundary value problem
SciPy’s fsolve: a root finder
SciPy special functions for elliptic integrals of the second kind: ellipeinc/ellipe
NumPhysPy’s Secant root finder
NumPhysPy’s LU decomposition and solver with pivot

The package has been verified by direct comparison with the output from the original python methods and this has been incorporated into a number of NUnit tests that are provided with the package.

The choice of methods and degrees to which the various modes and optional arguments are supported is largely based on what I find useful in my development of other code. It is not intended to be a “full port” of SciPy or NumPhysPy.

Code Documentation

Support: nbodyphysics@gmail.com

This is a free asset. If you value it please buy me a coffee and/or leave a review.

Contents

1 Numerical Integration Key Routines
2 Boundary Value Solver
3 Root Finding Key Routines
- 3.1 NumericalTools.Optimize.FsolveSimple
- 3.2 NumPhysPy.Secant
4 Matrix Solver
5 Example Scenes
6 Specifying a Differential Equation
- 6.1 Simple Harmonic Oscillator
- 6.2 Bessel Function

Numerical Integration Key Routines

NumericalTools.Ivp.SolveIvp

SolveIvp is the most direct interface to the SciPy solve_ivp differential equation solver. It is a port of the python code with some modes not implemented (complex numbers and vectorized mode are not supported). The available integration methods are RK45, RK23 and DOP853. Numerical Tools supports the use of events, t_eval and dense output. The result of invoking solve_ivp is a Ivp.OdeResult object.

The simplest use case for solve_ivp is to ask it to solve/numerically integrate a specified system of differential equations, given initial values from a time t0 to a time t1. The result of the function call is an Ivp.OdeResult class that holds the values of the system at time t1. An example of this use can be found in the QA unit tests: IvpBasicExample for each of the integrator types.

More sophisticated uses of solve_ivp may include the use of:

events (record the state when certain conditions occur, potentially terminating the integration)
t_eval (evaluate the solution at specific times)
dense_output (keep a record of intermediate values)

Events: The event system provides a mechanism to record the system state when a specified event function returns zero e.g. when a projectile reaches max height. The events that are triggered are available in the Ivp.OdeResult. The example SciPyBased/Scenes/Solve_Ivp/CannonWithTriggers provides demonstration code for events.

t_eval: The optional t_eval parameter can be used to pass in an array of times at which the system state should be recorded.

dense_output: In many cases the continuous evolution of the system between t0 and t1 is of interest. The dense output mode in solve_ivp can be used for this. Note that in general the DE solver will not take uniform steps; it will adjust step size to keep the error within specified bounds. Data from the solution at specific intervals is retrieved from the Ivp.OdeResult by providing a vector of times for which output values are needed. The data in the OdeResult is used to interpolate values for each of the entries in the vector of times. For an example see the LotkaVoletraPlot script and the associated demo scene.

NumericalTools.IvpRealtime

In many game/simulation scenarios the evolution of a system of DEs is best done in a continuous fashion aligned with the evolution of time in the Unity scene. Numercial Tools provides an adaptation of the SciPy solve_ivp method that allows this mode of operation: IvpRealtime. Typically an IvpRealtime object is initialized and then on each update a request is made to continue the integration to a new end time by using the IvpRealtime.EvolveTo() method. This returns the time evolved to and the state of the system. Note that due to the adaptive time step of SciPy integration methods the result may have evolved beyond the requested time and some interpolation may be required.

Examples of the direct use of IvpRealtime (as opposed to RealtimeIntegratorFactory) can be found in SciPyBased/Scenes/IvpRealtime. Events in the realtime use-case make use of a callback function that is called when the event function triggers.

RealtimeIntegratorFactory

The RealtimeIntegratorFactory is a wrapper that allow a choice of one of the SciPy realtime integration algorithms as well as two simpler, fixed timestep integration routines: Leapfrog and RK4_BASIC. Examples of its use can be found in RealtimeIntegrators/Scenes.

Boundary Value Solver

The SciPy solve_bvp method finds the solution of an ordinary differential equation given boundary values (e.g. an initial and a final position instead of an initial position and velocity). This is a more subtle and tricky part of numerical analysis and in some cases there may be more than one solution depending on seed values and in other cases there may be no solution. Most standard numerical books provide an overview of this topic.

The canonical examples from the SciPy documentation are provided both as demonstration scenes (Bvp_Bratu, Bvp_SturmL) and as NUnit tests in BvpTestBasic in which the output is compared point by point to data exported from Python.

Note: this is a significant coding effort and the case where a user-provided determinant is used has not been tested. The SciPy implementation used a sparse matrix representation and solver for iterating Newton’s method. The implementation here uses a full matrix implementation of LU decomposition and solving (with pivot) to avoid the need to port the Fortran based LAPACK sparse matrix package used within SciPy.

Root Finding Key Routines

NumericalTools.Optimize.FsolveSimple

Optimize.FsolveSimple finds the root of a function given an initial starting point of the function. The function delegate allows a vector valued function and attempts to determine an input vector value for which the function evaluates to zero.

NumPhysPy.Secant

An implementation of a very simple root finding algorithm for a single valued function.

Matrix Solver

The class LUBasic provides code to perform an LU decomposition with pivot and then a method to solve the system for a given right-hand side vector. See the LUTest NUnit test for a sample usage.

This solver is also used within the solve_bvp function.

Example Scenes

RealtimeIntegrators
- Oscillator: Uses OscillatorDE.cs to illustrate how to use the RealtimeIntegratorFactor to setup a simple harmonic oscillator differential equation.
- ChaosAizawa: Model the chaotic Aizawa DE demonstrating the way arguments to a system of DEs can be passed in via the factory class.
- TwoDGravity: Illustrates a Kepler orbit. Useful to see the impact of integrator and time step choice; the SciPy error control works well, fixed timestep methods do poorly.
IvpRealtime
- Pendulum: Shows an approximate and exact pendulum DE using the IvpRealtime wrapper
- RealtimeExp: Simple exponential DE evolved in realtime
- CannonRealtimeWithEvents: Shows a cannon ball in motion with events attached to the DE system to trigger callbacks at max height and when the ball the hits the ground
- SchwXYZ: Evolves the DE for the motion of photons in the presence of a Schwarzschild black hole.
Solve_Ivp
- SimpleHarmonicMotion: Solves and graphs a solution to a simple oscillator DE
- LotkaVolterra: Plots a solution to the Lotka-Voltara DE (model of populations of predator and prey)
- CannonWithTriggers: Models the vertical motion of a cannon ball with triggers to record the peak height and time when the ground is impacted
Solve_Bvp
- Bvp_Bratu: Plot two solutions for the SciPy Bratu example
- Bvp_SturmL: Plot the solution for the Sturm-Lousville example from the SciPy docs
Fsolve/Secant
- FsolveSimple: Example shows how to use Fsolve and Secant to find the roots of exp(x-sqrt(x)) – x =0

Specifying a Differential Equation

Simple Harmonic Oscillator

The differential equation solvers in Numerical Tools expect a function that evaluates the first derivatives of each entry in the state vector y[]. This requires any equation with higher derivatives to be converted to a system of first order differential equations. Information on this can be found in many online lectures and the standard textbooks on ordinary differential equations.

As a simple example consider the simple harmonic oscillator:

$\huge{ \frac{d^2x}{dt^2} + x = 0}$

This is a second order equation for one variable x(t). It can be rewritten as two equations: one for velocity $v = frac{dx}{dt}$ and one for x(t). We can now look at the equations:

$\frac{dv}{dt} = -x$

and

$\frac{dx}{dt} = v$

where we have used the defintion of v as a second equation. The required function for Ivp.SolveIvp or IvpRealtime is:

public double[] Oscillator(double t, double[] y, params object[] args)
    {
        y_prime[0] = y[1];
        y_prime[1] = -y[0];
        return y_prime;
    }

In which y[0] holds the x variable and y[1] holds v.

Bessel Function

A more complex example is the Bessel function defined as the solution to:

$x^2 \frac{d^2y}{dx^2} + x \frac{dy}{dx} + (x^2 - a^2)y = 0$

This can be expressed by letting $z=\frac{dy}{dx}$ and $\frac{dz}{dx} = \frac{d^2y}{dx^2}$ . This leads to the equations:

$\frac{dy}{dx} = z$

and rearranging Bessel’s equation for $\frac{d^2y}{dx^2}$

$\frac{dz}{dx} = \frac{1}{x^2} \left( -x * z - (x^2 -a^2) * y \right)$ .

SciPy example of Bessel’s equation.

The resulting DE function for Numerical Tools is:

private double[] Bessel(double t, double[] y, params object[] args)
    {
        double a = (double)args[0];

        y_prime[0] = y[1];
        y_prime[1] = 1.0 / (t * t) * (-t * y[1] - (t*t - a * a) * y[0]);
        return y_prime;
    }

Note that since the equation divides by the initial value t it is important to begin evolution at some small value larger than 0.

N-body Physics

and other interesting stuff…