GitHub

Convergence Simulations

The convergence simulation type is useful for deriving order of convergence estimates from a group of simulations. This object will automatically assemble error vectors into a more useful manner and provide plotting functionality. Convergence estimates are also given by pair-wise estimates.

One can automatically have DifferentialEquations.jl perform the error analysis by passing a ConvergenceSimulation a vector of solutions, or using one of the provided test_convergence functions. These will give order of convergence estimates and provide plotting functionality. This requires that the true solution was provided in the problem definition.

ConvergenceSimulations can either be created by passing the constructor the appropriate solution array or by using one of the provided test_convergence functions.

The ConvergenceSimulation Type

A type which holds the data from a convergence simulation.

Fields

  • solutions::Array{<:DESolution}: Holds the solutions. For the Monte-Carlo method of test_convergence these hold a single representative trajectory unless retain_solutions = true is passed — see Memory use of Monte-Carlo studies.

  • errors: Dictionary of the error calculations. Can contain:

    • h1Errors: Vector of the H1 errors.
    • l2Errors: Vector of the L2 errors.
    • maxErrors: Vector of the nodal maximum errors.
    • node2Errors: Vector of the nodal l2 errors.

  • N: The number of simulations.

  • auxdata: Auxiliary data of the convergence simulation. Entries can include:

    • dts: The dt's in the simulations.
    • dxs: The dx's in the simulations.
    • μs: The CFL μ's in the simulations.
    • νs: The CFL ν's in the simulations.

  • 𝒪est: Dictionary of order estimates. Can contain:

    • ConvEst_h1: The H1 error order of convergence estimate for the convergence simulation. Generated via log2(error[i+1]/error[i]). Thus only valid if generated by halving/doubling the dt/dx. If alternate scaling, modify by dividing of log(base,ConvEst_h1)
    • ConvEst_l2: The L2 error order of convergence estimate for the convergence simulation. Generated via log2(error[i+1]/error[i]). Thus only valid if generated by halving/doubling the dt/dx. If alternate scaling, modify by dividing of log(base,ConvEst_l2)
    • ConvEst_max: The nodal maximum error order of convergence estimate for the convergence simulation. Generated via log2(error[i+1]/error[i]). Thus only valid if generated by halving/doubling the dt/dx. If alternate scaling, modify by dividing of log(base,ConvEst_max)
    • ConvEst_node2: The nodal l2 error order of convergence estimate for the convergence simulation. Generated via log2(error[i+1]/error[i]). Thus only valid if generated by halving/doubling the dt/dx. If alternate scaling, modify by dividing of log(base,ConvEst_node2)

  • convergence_axis: The axis along which convergence is calculated. For example, if we calculate the dt convergence, convergence_axis is the dts used in the calculation.

Plot Functions

The plot functionality is provided by a Plots.jl recipe. What is plotted is a line series for each calculated error along the convergence axis. To plot a convergence simulation, simply use:

plot(sim::ConvergenceSimulation)

All of the functionality (keyword arguments) provided by Plots.jl are able to be used in this command. Please see the Plots.jl documentation for more information.

ODE

test_convergence(dts::AbstractArray,prob::AbstractODEProblem)

Tests the order of the time convergence of the given algorithm on the given problem solved over the given dts. Keyword arguments are passed to the ODE solver.

SDE

test_convergence(dts::AbstractArray,prob::AbstractSDEProblem)

Tests the strong order time convergence of the given algorithm on the given problem solved over the given dts. Keyword arguments are passed to the ODE solver. Except:

  • numMonte: The number of simulations for each dt. Default is 10000.

Order Estimation

calc𝒪estimates(error::Vector{Number})`

Computes the pairwise convergence estimate for a convergence test done by halving/doubling stepsizes via

log2(error[i+1]/error[i])

Returns the mean of the convergence estimates.

Memory use of Monte-Carlo studies

For an SDE or RODE problem, test_convergence runs an ensemble per step size. A full solution per trajectory carries the solver cache, the noise process, the problem and its interpolation, so keeping them makes a study over many trajectories retain far more than the error estimates it reports — at hundreds of thousands of trajectories, tens of gigabytes.

They are therefore not kept by default. Each trajectory is reduced to a ConvergenceTrajectory by the ensemble's output_func as it is solved, so the full solutions never coexist, and once the errors have been computed the ensemble is stripped to a single representative trajectory. errors, weak_errors, error_means and 𝒪est come from the full ensemble and are unaffected — what changes is that solutions[i].u holds one entry rather than trajectories of them.

Pass retain_solutions = true when the trajectories themselves are the point, such as comparing two algorithms path by path. sim[i, j] forwards into solutions[i][j], so it too reaches only the representative trajectory unless the solutions are retained.

The reduction is skipped where it cannot apply, and the solutions retained: when you supply your own EnsembleProblem (set its output_func instead), when expected_value is given, and when weak_timeseries_errors or weak_dense_errors is set.

API

DiffEqDevTools.ConvergenceSimulation — Type
ConvergenceSimulation(
    solutions, convergence_axis;
    auxdata = nothing, additional_errors = nothing, expected_value = nothing
)

Collect solutions and error series from a convergence experiment. The constructor combines the error measurements from solutions with additional_errors and estimates the observed convergence order between adjacent points on convergence_axis.

Set expected_value when the error series are supplied through additional_errors, as done by ensemble convergence tests.

source
DiffEqDevTools.ConvergenceTrajectory — Type
ConvergenceTrajectory(t, u, u_analytic, errors)

The part of a trajectory's solution that a convergence study actually consumes: its error measurements and the endpoint values the weak errors are formed from.

test_convergence stores one of these per trajectory in place of the full solution unless retain_solutions = true, which is what keeps a Monte-Carlo study's memory proportional to the errors it reports rather than to the solver state it happened to allocate along the way.

source
DiffEqDevTools.test_convergence — Function
test_convergence(dts, prob, alg[, ensemblealg]; kwargs...)
test_convergence(probs, convergence_axis, alg; kwargs...)
test_convergence(setup::ConvergenceSetup, alg; kwargs...)

Solve a problem at each point on a convergence axis and return a ConvergenceSimulation. For time-step convergence tests, each value in dts is passed to the solver as dt. Remaining keyword arguments are forwarded to solve.

The computed solutions must contain error measurements, usually obtained from an analytic solution. Use analyticless_test_convergence when only a numerical reference solution is available.

Keyword Arguments

  • retain_solutions: whether ConvergenceSimulation.solutions keeps the trajectory solutions (default false). A Monte-Carlo study over many trajectories otherwise retains one complete solution object — solver cache, noise process, problem and interpolation — per trajectory per step size, which for large trajectories is orders of magnitude more memory than the error estimates it is computing, and is what put the weak-convergence test groups past a 16 GB runner.

    By default each trajectory is reduced to a ConvergenceTrajectory by the ensemble's output_func as it is solved, so the full solutions never coexist, and the reduced ensemble is stripped to one representative trajectory once the errors have been computed from it. errors, weak_errors, error_means and 𝒪est are computed from the full ensemble either way and are unaffected; what changes is that solutions[i].u holds one entry rather than trajectories of them.

    Pass true when the trajectories themselves are the point, for example to compare two algorithms path by path.

    The reduction is skipped, and the solutions retained, where it cannot apply: when you supply your own EnsembleProblem (set its output_func instead), when expected_value is given (the weak error is formed from the trajectory values directly), and when weak_timeseries_errors or weak_dense_errors is set (both need the whole timeseries).

source
DiffEqDevTools.analyticless_test_convergence — Function
analyticless_test_convergence(dts, prob, alg, reference; kwargs...)

Estimate convergence against a numerical reference solution and return a ConvergenceSimulation. For ODE problems, reference is a solver setup whose :alg entry selects the reference algorithm. For stochastic problems, reference is the finer reference step size, and each tested solution reuses the reference noise realization.

source