```@meta
CollapsedDocStrings = true
```

# OrdinaryDiffEqFIRK

Fully Implicit Runge-Kutta (FIRK) methods for stiff differential equations requiring very high accuracy. These methods solve a fully coupled implicit system at each timestep, providing superior accuracy and stability compared to diagonally implicit methods.

!!! warning "Real Numbers Only"
    
    FIRK methods should only be used for problems defined on real numbers, not complex numbers.

## Key Properties

FIRK methods provide:

  - **Highest-order implicit methods** (excluding extrapolation)
  - **Superior accuracy** for very low tolerance requirements (≤ 1e-9)
  - **A-stable and L-stable** behavior for stiff problems
  - **Higher order per stage** than SDIRK methods (order 2s or 2s-1 for s stages)
  - **Geometric properties** can be displayed by some methods in certain cases
  - **Excellent for small to medium systems** with high accuracy requirements

## When to Use FIRK Methods

These methods are recommended for:

  - **Very low tolerance problems** (1e-9 and below) where accuracy is paramount
  - **Small to medium stiff systems** (< 200 equations)
  - **Problems requiring highest possible accuracy** for implicit methods
  - **Stiff problems** where SDIRK order limitations (max order 5) are insufficient
  - **Applications where computational cost is acceptable** for maximum accuracy

## Mathematical Background

RadauIIA methods are based on Gaussian collocation and achieve order 2s-1 for s stages, making them among the highest-order implicit methods available. They represent the ODE analog of Gaussian quadrature with a fixed right node (Gauss-Radau). They are L-stable, meaning they are stiffly accurate and handle large negative eigenvalues/stiffness ratios. For more details on recent advances in FIRK methods, see our paper: [High-Order Adaptive Time Stepping for the Incompressible Navier-Stokes Equations](https://arxiv.org/abs/2412.14362).

Gauss–Legendre methods use the roots of the Legendre polynomial directly as collocation points and achieve an order of 2s, the highest order IRK available. They are also A-stable, symmetric, and can be symplectic in certain usages. They conserve energy within the system but are not L-Stable like Radau, which means they preserve oscillations instead of damping. 

## Computational Considerations

### Advantages

  - **Higher accuracy per stage** than diagonal methods
  - **Better multithreading** for small systems due to larger linear algebra operations
  - **No order restrictions** like SDIRK methods (which max out at order 5)

### Disadvantages

  - **Limited to real-valued problems** - cannot be used for complex number systems
  - **Higher implementation complexity** compared to SDIRK methods

## Solver Selection Guide

### High accuracy requirements

  - **`AdaptiveRadau`**: **Recommended** - adaptive order method that automatically selects optimal order
  - **`RadauIIA5`**: 5th-order method, good balance of accuracy and efficiency
  - **`RadauIIA9`**: 9th-order method for extremely high accuracy requirements
  - **`RadauIIA3`**: 3rd-order method for moderate accuracy needs
  - **`GaussLegendre`**: 2s order method for high accuracy in nonstiff or conservative problems

### System size considerations

  - **Systems < 200**: FIRK methods are competitive due to better multithreading
  - **Systems > 200**: Consider SDIRK or BDF methods instead, or keep a Radau method and
    switch to a matrix-free linear solver (see below)

## Matrix-free (Newton-Krylov) usage

The RadauIIA methods and `AdaptiveRadau` accept a matrix-free `linsolve`, which avoids
forming a Jacobian at all:

```julia
solve(prob, RadauIIA5(linsolve = KrylovJL_GMRES()))
```

The W-transform decouples each step into one real shifted solve and one shifted solve per
complex-conjugate stage pair, and every one of them is applied through the same real
Jacobian-vector product. `concrete_jac = true` keeps a concrete Jacobian alongside the
matrix-free application, which is what a preconditioner needs.

`GaussLegendre` has no matrix-free path: it solves a single coupled
`num_stages * length(u)` block system rather than decoupled shifted systems.

## Performance Guidelines

  - **Best for tolerances ≤ 1e-9** where high accuracy justifies the cost
  - **Most efficient on small to medium systems** where linear algebra cost is manageable
  - **Should be tested against** parallel implicit extrapolation methods which specialize in similar regimes
  - **Compare with** high-order SDIRK methods for borderline cases

```@eval
first_steps = evalfile("./common_first_steps.jl")
first_steps("OrdinaryDiffEqFIRK", "RadauIIA5")
```

## Full list of solvers

```@docs
RadauIIA3
RadauIIA5
RadauIIA9
AdaptiveRadau
GaussLegendre
```
