Dealing with numerical instability

Internally, the FRAGMENT-MNP model uses a numerical solver to solve an initial value problem for the system of ODEs that form the model (see FragmentMNP.run()).

For certain input parameter ranges, care needs to be taken to avoid the solution becoming unstable. This is particularly relevant when mass concentrations in a certain size class trend asymptotically towards zero, in which case the system can be described as “stiff”. The model provides a few ways of dealing with this:

• Changing the integration method used by specifying the solver_method config option. The available methods are those provided by SciPy. “Radau”, “BDF” and “LSODA” are particularly well suited to stiff problems.

• Changing the max step size by specifying the solver_max_step config option.

• Changing the error tolerances by specifying the solver_rtol and solver_atol config options.

Example of a numerically unstable problem

Let’s create a model setup that will be numerically unstable. We do this by setting a relatively high dissolution rate \(k_\text{diss}\). First, let’s use the Runge-Kutta order 5(4) solver (RK45), which is an explicit solver and therefore not particularly good at dealing with stiff problems. This is the default option, and so we don’t need to specify it explicitly as a config option.

import numpy as np
from fragmentmnp import FragmentMNP
from fragmentmnp.examples import minimal_config, minimal_data

# Keeping it simply with two size classes
minimal_config['n_size_classes'] = 2
minimal_config['particle_size_range'] = [-6, -5]

# Arbitrary initial concs and a high dissolution rate
minimal_data['initial_concs'] = [10, 10]
minimal_data['k_frag'] = 0.1
minimal_data['k_diss'] = 0.3

# Run the model and plot the output
fmnp = FragmentMNP(minimal_config, minimal_data)
output = fmnp.run()
_ = output.plot(plot_dissolution=True, log_yaxis=True)

Something is clearly going wrong here! Let’s try LSODA instead, which is better at dealing with stiffness:

# Set the solver method to LSODA, which is particularly suited
# to dealing with stiff problems
minimal_config['solver_method'] = 'LSODA'

fmnp = FragmentMNP(minimal_config, minimal_data)
output = fmnp.run()
_ = output.plot(plot_dissolution=True, log_yaxis=True)

Better, but not quite there yet. If we look at the SciPy solution object, accessed via output.soln, we can see the number of evaluations made (nfev):

output.soln.nfev

133

Let’s try and set a smaller solver_max_step and thereby increase this number of evaluations, to see if that helps:

# Set the maximum step size (i.e. timestep) for the ODE solver
minimal_config['solver_max_step'] = 0.3

fmnp = FragmentMNP(minimal_config, minimal_data)
output = fmnp.run()
output.plot(plot_dissolution=True, log_yaxis=True)
output.soln.nfev

568

This seems to have done the trick! Alternatively, we could have set smaller atol and rtol than the defaults (10e-6 and 10e-3, respectively):

# Don't set a max step size, but set small error tolerances instead
minimal_config['solver_max_step'] = np.inf
minimal_config['solver_atol'] = 10e-14
minimal_config['solver_rtol'] = 10e-12

fmnp = FragmentMNP(minimal_config, minimal_data)
output = fmnp.run()
output.plot(plot_dissolution=True, log_yaxis=True)
output.soln.nfev

543

Note that we have to set considerably smaller tolerances than the defaults to solve the instability.

Numerical instability due to timestep size

The model does not do any scaling of timesteps when they are input to the numerical solver. Whilst this shouldn’t normally cause issues, specifying very large timestep lengths dt may be a problem. To alleviate this, you might want to consider using a different unit for your time dimension. For example, if you want your timesteps to be 1 year long, then instead of dt=31526000 (the number of seconds in a year), use dt=1 and provide rate constants (k_frag and k_diss) with units of /year instead of /s. Rate constant scaling parameters most also follow your unit conventions.