BoundaryValueDiffEqFIRK

Fully Implicit Runge Kutta(FIRK) Methods. To be able to access the solvers in BoundaryValueDiffEqFIRK, you must first install them use the Julia package manager:

using Pkg
Pkg.add("BoundaryValueDiffEqFIRK")

solve(prob::BVProblem, alg, dt; kwargs...)
solve(prob::TwoPointBVProblem, alg, dt; kwargs...)

Nested nonlinear solving in FIRK methods

When working with large boundary value problems, especially those involving stiff systems, computational efficiency and solver robustness become critical concerns. To improve the efficiency of FIRK methods on large BVPs, we can use nested nonlinear solving to obtain the implicit FIRK step instead of solving them as part of the global residual. In BoundaryValueDiffEq.jl, we can set nested_nlsolve as true to enable FIRK methods to compute the implicit FIRK steps using nested nonlinear solving(default option in FIRK methods is nested_nlsolve=false).

Moreover, the nested nonlinear problem solver can be finely tuned to meet specific accuracy requirements by providing detailed keyword arguments through the nested_nlsolve_kwargs option in any FIRK solver, for example, RadauIIa5(; nested_nlsolve = true, nested_nlsolve_kwargs = (; abstol = 1e-6, reltol = 1e-6)), where nested_nlsolve_kwargs can be any common keyword arguments in NonlinearSolve.jl, see Common Solver Options in NonlinearSolve.jl.

Full List of Methods

Radau IIA methods

RadauIIa1: 1 stage Radau IIA method, without defect control adaptivity
RadauIIa2: 2 stage Radau IIA method, with defect control adaptivity.
RadauIIa3: 3 stage Radau IIA method, with defect control adaptivity.
RadauIIa5: 5 stage Radau IIA method, with defect control adaptivity.
RadauIIa7: 7 stage Radau IIA method, with defect control adaptivity.

Lobatto IIIA methods

LobattoIIIa2: 2 stage Lobatto IIIa method, with defect control adaptivity.
LobattoIIIa3: 3 stage Lobatto IIIa method, with defect control adaptivity.
LobattoIIIa4: 4 stage Lobatto IIIa method, with defect control adaptivity.
LobattoIIIa5: 5 stage Lobatto IIIa method, with defect control adaptivity.

Lobatto IIIB methods

LobattoIIIb2: 2 stage Lobatto IIIb method, without defect control adaptivity.
LobattoIIIb3: 3 stage Lobatto IIIb method, with defect control adaptivity.
LobattoIIIb4: 4 stage Lobatto IIIb method, with defect control adaptivity.
LobattoIIIb5: 5 stage Lobatto IIIb method, with defect control adaptivity.

Lobatto IIIC methods

LobattoIIIc2: 2 stage Lobatto IIIc method, without defect control adaptivity.
LobattoIIIc3: 3 stage Lobatto IIIc method, with defect control adaptivity.
LobattoIIIc4: 4 stage Lobatto IIIc method, with defect control adaptivity.
LobattoIIIc5: 5 stage Lobatto IIIc method, with defect control adaptivity.

Detailed Solvers Explanation

BoundaryValueDiffEqFIRK.RadauIIa1 — Type

  RadauIIa1(; nlsolve = NewtonRaphson(), jac_alg = BVPJacobianAlgorithm(), nested_nlsolve = false, nest_tol = 0.0,
          defect_threshold = 0.1, max_num_subintervals = 3000)

1th stage RadauIIa method.

Keyword Arguments

nlsolve: Internal Nonlinear solver. Any solver which conforms to the SciML NonlinearProblem interface can be used. Note that any autodiff argument for the solver will be ignored and a custom jacobian algorithm will be used.
jac_alg: Jacobian Algorithm used for the nonlinear solver. Defaults to BVPJacobianAlgorithm(), which automatically decides the best algorithm to use based on the input types and problem type.
- For TwoPointBVProblem, only diffmode is used (defaults to AutoSparse(AutoForwardDiff) if possible else AutoSparse(AutoFiniteDiff)).
- For BVProblem, bc_diffmode and nonbc_diffmode are used. For nonbc_diffmode defaults to AutoSparse(AutoForwardDiff) if possible else AutoSparse(AutoFiniteDiff). For bc_diffmode, defaults to AutoForwardDiff if possible else AutoFiniteDiff.
nested_nlsolve: Whether or not to use a nested nonlinear solve for the implicit FIRK step. Defaults to false. If set to false, the FIRK stages are solved as a part of the global residual. The general recommendation is to choose true for larger problems and false for smaller ones.
nest_tol: The tolerance for the nested solver. Default is nothing which leads to NonlinearSolve automatically selecting the tolerance.
defect_threshold: Threshold for defect control.
max_num_subintervals: Number of maximal subintervals, default as 3000.

Note

For type-stability, the chunksizes for ForwardDiff ADTypes in BVPJacobianAlgorithm must be provided.

References

Reference for Lobatto and Radau methods:

@incollection{Jay2015,
    author="Jay, Laurent O.",
    editor="Engquist, Bj{"o}rn",
    title="Lobatto Methods",
    booktitle = {Encyclopedia of {Applied} and {Computational} {Mathematics}},
    year="2015",
    publisher="Springer Berlin Heidelberg",
}
@incollection{engquist_radau_2015,
    author = {Hairer, Ernst and Wanner, Gerhard},
    editor={Engquist, Bj{"o}rn},
    title = {Radau {Methods}},
    booktitle = {Encyclopedia of {Applied} and {Computational} {Mathematics}},
    publisher = {Springer Berlin Heidelberg},
    year = {2015},
}

References for implementation of defect control, based on the bvp5c solver in MATLAB:

@article{shampine_solving_nodate,
    title = {Solving {Boundary} {Value} {Problems} for {Ordinary} {Diﬀerential} {Equations} in {Matlab} with bvp4c
    author = {Shampine, Lawrence F and Kierzenka, Jacek and Reichelt, Mark W},
    year = {2000},
}

@article{kierzenka_bvp_2008,
    title = {A {BVP} {Solver} that {Controls} {Residual} and {Error}},
    author = {Kierzenka, J and Shampine, L F},
    year = {2008},
}

@article{russell_adaptive_1978,
    title = {Adaptive {Mesh} {Selection} {Strategies} for {Solving} {Boundary} {Value} {Problems}},
    journal = {SIAM Journal on Numerical Analysis},
    author = {Russell, R. D. and Christiansen, J.},
    year = {1978},
}