# Timestepping Method Descriptions

## Common Setup

All methods start by calculating a scaled error estimate on each scalar component of $u$:

\[err^{scaled}_i = norm(err_i/(abstol_i + max(uprev_i,u_i)reltol_i))\]

On this scaled error estimate, we calculate the norm. This norm is usually the Hairer norm:

\[norm(x) = sqrt(sum(x^2)/length(x))\]

This norm works well because it does not change if we add new pieces to the differential equation: it scales our error by the number of equations so that independent equations will not step differently than a single solve.

In all cases, the step is rejected if $err^{scaled}>1$ since that means the error is larger than the tolerances, and the step is accepted if $err^{scaled}<1$.

## Integral Controller (Standard Controller)

The proportional control algorithm is the "standard algorithm" for adaptive timestepping. Note that it is not the default in DifferentialEquations.jl because it is usually awful for performance, but it is explained first because it is the most widely taught algorithm and others build off of its techniques.

The control simply changes `dt`

proportional to the error. There is an exponentiation based on the order of the algorithm which goes back to a result by Cechino for the optimal stepsize to reduce the error. The algorithm is:

```
qtmp = integrator.EEst^(1/(alg_adaptive_order(integrator.alg)+1))/integrator.opts.gamma
@fastmath q = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),qtmp))
integrator.dtnew = integrator.dt/q
```

Thus `q`

is the scaling factor for `dt`

, and it must be between `qmin`

and `qmax`

. `gamma`

is the safety factor, `0.9`

, for how much `dt`

is decreased below the theoretical "optimal" value.

Since proportional control is "jagged", i.e. can cause large changes between one step to the next, it can effect the stability of explicit methods. Thus it's only applied by default to low order implicit solvers.

`OrdinaryDiffEq.IController`

— Type`IController()`

The standard (integral) controller is the most basic step size controller. This controller is usually the first one introduced in numerical analysis classes but should only be used rarely in practice because of efficiency problems for many problems/algorithms.

Construct an integral (I) step size controller adapting the time step based on the formula

`Δtₙ₊₁ = εₙ₊₁^(1/k) * Δtₙ`

where `k = get_current_adaptive_order(alg, integrator.cache) + 1`

and `εᵢ`

is the inverse of the error estimate `integrator.EEst`

scaled by the tolerance (Hairer, Nørsett, Wanner, 2008, Section II.4). The step size factor is multiplied by the safety factor `gamma`

and clipped to the interval `[qmin, qmax]`

. A step will be accepted whenever the estimated error `integrator.EEst`

is less than or equal to unity. Otherwise, the step is rejected and re-tried with the predicted step size.

**References**

- Hairer, Nørsett, Wanner (2008) Solving Ordinary Differential Equations I Nonstiff Problems DOI: 10.1007/978-3-540-78862-1

## Proportional-Integral Controller (PI Controller)

The proportional-integral control algorithm is a standard control algorithm from control theory. It mixes proportional control with memory in order to make the timesteps more stable, which actually increases the adaptive stability region of the algorithm. This stability property means that it's well-suited for explicit solvers, and it's applied by default to the Rosenbrock methods as well. The form for the updates is:

```
EEst,beta1,q11,qold,beta2 = integrator.EEst, integrator.opts.beta1, integrator.q11,integrator.qold,integrator.opts.beta2
@fastmath q11 = EEst^beta1
@fastmath q = q11/(qold^beta2)
integrator.q11 = q11
@fastmath q = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),q/integrator.opts.gamma))
if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min
q = one(q)
end
q
```

`beta1`

is the gain on the proportional part, and `beta2`

is the gain for the history portion. `qoldinit`

is the initialized value for the gain history.

`OrdinaryDiffEq.PIController`

— Type`PIController(beta1, beta2)`

The proportional-integral (PI) controller is a widespread step size controller with improved stability properties compared to the `IController`

. This controller is the default for most algorithms in OrdinaryDiffEq.jl.

Construct a PI step size controller adapting the time step based on the formula

`Δtₙ₊₁ = εₙ₊₁^β₁ * εₙ^β₂ * Δtₙ`

where `εᵢ`

are inverses of the error estimates scaled by the tolerance (Hairer, Nørsett, Wanner, 2010, Section IV.2). The step size factor is multiplied by the safety factor `gamma`

and clipped to the interval `[qmin, qmax]`

. A step will be accepted whenever the estimated error `integrator.EEst`

is less than or equal to unity. Otherwise, the step is rejected and re-tried with the predicted step size.

The coefficients `beta1, beta2`

are not scaled by the order of the method, in contrast to the `PIDController`

. For the `PIController`

, this scaling by the order must be done when the controller is constructed.

**References**

- Hairer, Nørsett, Wanner (2010) Solving Ordinary Differential Equations II Stiff and Differential-Algebraic Problems DOI: 10.1007/978-3-642-05221-7
- Hairer, Nørsett, Wanner (2008) Solving Ordinary Differential Equations I Nonstiff Problems DOI: 10.1007/978-3-540-78862-1

## Proportional-Integral-Derivative Controller (PID Controller)

`OrdinaryDiffEq.PIDController`

— Type```
PIDController(beta1, beta2, beta3=zero(beta1);
limiter=default_dt_factor_limiter,
accept_safety=0.81)
```

The proportional-integral-derivative (PID) controller is a generalization of the `PIController`

and can have improved stability and efficiency properties.

Construct a PID step size controller adapting the time step based on the formula

`Δtₙ₊₁ = εₙ₊₁^(β₁/k) * εₙ^(β₂/k) * εₙ₋₁^(β₃/ k) * Δtₙ`

where `k = min(alg_order, alg_adaptive_order) + 1`

and `εᵢ`

are inverses of the error estimates scaled by the tolerance (Söderlind, 2003). The step size factor is limited by the `limiter`

with default value

`limiter(x) = one(x) + atan(x - one(x))`

as proposed by Söderlind and Wang (2006). A step will be accepted whenever the predicted step size change is bigger than `accept_safety`

. Otherwise, the step is rejected and re-tried with the predicted step size.

Some standard controller parameters suggested in the literature are

Controller | `beta1` | `beta2` | `beta3` |
---|---|---|---|

basic | `1.00` | `0.00` | `0` |

PI42 | `0.60` | `-0.20` | `0` |

PI33 | `2//3` | `-1//3` | `0` |

PI34 | `0.70` | `-0.40` | `0` |

H211PI | `1//6` | `1//6` | `0` |

H312PID | `1//18` | `1//9` | `1//18` |

In contrast to the `PIController`

, the coefficients `beta1, beta2, beta3`

are scaled by the order of the method. Thus, standard controllers such as PI42 can use the same coefficients `beta1, beta2, beta3`

for different algorithms.

In contrast to other controllers, the `PIDController`

does not use the keyword arguments `qmin, qmax`

to limit the step size change or the safety factor `gamma`

. These common keyword arguments are replaced by the `limiter`

and `accept_safety`

to guarantee a smooth behavior (Söderlind and Wang, 2006). Because of this, a `PIDController`

behaves different from a `PIController`

, even if `beta1, beta2`

are adapted accordingly and `iszero(beta3)`

.

**References**

- Söderlind (2003) Digital Filters in Adaptive Time-Stepping DOI: 10.1145/641876.641877
- Söderlind, Wang (2006) Adaptive time-stepping and computational stability DOI: 10.1016/j.cam.2005.03.008
- Ranocha, Dalcin, Parsani, Ketcheson (2021) Optimized Runge-Kutta Methods with Automatic Step Size Control for Compressible Computational Fluid Dynamics arXiv:2104.06836

## Gustafsson Acceleration

The Gustafsson acceleration algorithm accelerates changes so that way algorithms can more swiftly change to handle quick transients. This algorithm is thus well-suited for stiff solvers where this can be expected, and is the default for algorithms like the (E)SDIRK methods.

```
gamma = integrator.opts.gamma
niters = integrator.cache.newton_iters
fac = min(gamma,(1+2*integrator.alg.max_newton_iter)*gamma/(niters+2*integrator.alg.max_newton_iter))
expo = 1/(alg_order(integrator.alg)+1)
qtmp = (integrator.EEst^expo)/fac
@fastmath q = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),qtmp))
if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min
q = one(q)
end
integrator.qold = q
q
```

In this case, `niters`

is the number of Newton iterations which was required in the most recent step of the algorithm. Note that these values are used differently depending on acceptance and rejection. When the step is accepted, the following logic is applied:

```
if integrator.success_iter > 0
expo = 1/(alg_adaptive_order(integrator.alg)+1)
qgus=(integrator.dtacc/integrator.dt)*(((integrator.EEst^2)/integrator.erracc)^expo)
qgus = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),qgus/integrator.opts.gamma))
qacc=max(q,qgus)
else
qacc = q
end
integrator.dtacc = integrator.dt
integrator.erracc = max(1e-2,integrator.EEst)
integrator.dt/qacc
```

When it rejects, its the same as the proportional control:

```
if integrator.success_iter == 0
integrator.dt *= 0.1
else
integrator.dt = integrator.dt/integrator.qold
end
```

`OrdinaryDiffEq.PredictiveController`

— Type`PredictiveController()`

The Gustafsson acceleration algorithm accelerates changes so that way algorithms can more swiftly change to handle quick transients. This algorithm is thus well-suited for stiff solvers where this can be expected, and is the default for algorithms like the (E)SDIRK methods.

```
gamma = integrator.opts.gamma
niters = integrator.cache.newton_iters
fac = min(gamma,(1+2*integrator.alg.max_newton_iter)*gamma/(niters+2*integrator.alg.max_newton_iter))
expo = 1/(alg_order(integrator.alg)+1)
qtmp = (integrator.EEst^expo)/fac
@fastmath q = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),qtmp))
if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min
q = one(q)
end
integrator.qold = q
q
```

In this case, `niters`

is the number of Newton iterations which was required in the most recent step of the algorithm. Note that these values are used differently depending on acceptance and rejectance. When the step is accepted, the following logic is applied:

```
if integrator.success_iter > 0
expo = 1/(alg_adaptive_order(integrator.alg)+1)
qgus=(integrator.dtacc/integrator.dt)*(((integrator.EEst^2)/integrator.erracc)^expo)
qgus = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),qgus/integrator.opts.gamma))
qacc=max(q,qgus)
else
qacc = q
end
integrator.dtacc = integrator.dt
integrator.erracc = max(1e-2,integrator.EEst)
integrator.dt/qacc
```

When it rejects, it's the same as the `IController`

:

```
if integrator.success_iter == 0
integrator.dt *= 0.1
else
integrator.dt = integrator.dt/integrator.qold
end
```