|
1 |
| -# DiffEqOperators |
| 1 | +# Matrix-Free Linear Operators and Specializations on Linearity |
2 | 2 |
|
3 |
| -The `AbstractDiffEqOperator` interface is an interface for declaring parts of |
4 |
| -a differential equation as linear or affine. This then allows the solvers to |
5 |
| -exploit linearity to achieve maximal performance. |
6 |
| - |
7 |
| -## Using DiffEqOperators |
8 |
| - |
9 |
| -`AbstractDiffEqOperator`s act like functions. When defined, `A` has function |
10 |
| -calls `A(u,p,t)` and `A(du,u,p,t)` that act like `A*u`. These operators update |
11 |
| -via a function `update_coefficients!(A,u,p,t)`. |
12 |
| - |
13 |
| -## Constructors |
14 |
| - |
15 |
| -### Wrapping an Array: DiffEqArrayOperator |
16 |
| - |
17 |
| -`DiffEqArrayOperator` is for defining an operator directly from an array. The |
18 |
| -operator is of the form: |
19 |
| - |
20 |
| -```math |
21 |
| -A(u,p,t) |
22 |
| -``` |
23 |
| - |
24 |
| -for some scalar `α` and time plus possibly state dependent `A`. The |
25 |
| -constructor is: |
26 |
| - |
27 |
| -```julia |
28 |
| -DiffEqArrayOperator(A::AbstractMatrix{T}, update_func = DEFAULT_UPDATE_FUNC) |
29 |
| -``` |
30 |
| - |
31 |
| -`A` is the operator array. `update_func(A,u,p,t)` is the function called by |
32 |
| -`update_coefficients!(A,u,p,t)`. If left as its default, then `update_func` |
33 |
| -is trivial, which signifies `A` is a constant. |
34 |
| - |
35 |
| -### AffineDiffEqOperator |
36 |
| - |
37 |
| -For `As = (A1,A2,...,An)` and `Bs = (B1,B2,...,Bm)` where each of the `Ai` and |
38 |
| -`Bi` are `AbstractDiffEqOperator`s, the following constructor: |
39 |
| - |
40 |
| -```julia |
41 |
| -AffineDiffEqOperator{T}(As, Bs, u_cache = nothing) |
42 |
| -``` |
43 |
| - |
44 |
| -builds an operator `L = (A1 + A2 + ... An)*u + B1 + B2 + ... + Bm`. `u_cache` |
45 |
| -is for designating a type of internal cache for non-allocating evaluation of |
46 |
| -`L(du,u,p,t)`. If not given, the function `L(du,u,p,t)` is not available. Note |
47 |
| -that in solves which exploit this structure, this function call is not necessary. |
48 |
| -It's only used as the fallback in ODE solvers, which were not developed for this |
49 |
| -structure. |
50 |
| - |
51 |
| -## Formal Properties of DiffEqOperators |
52 |
| - |
53 |
| -These are the formal properties that an `AbstractDiffEqOperator` should obey |
54 |
| -for it to work in the solvers. |
55 |
| - |
56 |
| -### AbstractDiffEqOperator Interface Description |
57 |
| - |
58 |
| - 1. Function call and multiplication: `L(du,u,p,t)` for inplace and `du = L(u,p,t)` for |
59 |
| - out-of-place, meaning `L*u` and `mul!`. |
60 |
| - 2. If the operator is not a constant, update it with `(u,p,t)`. A mutating form, i.e., |
61 |
| - `update_coefficients!(A,u,p,t)` that changes the internal coefficients, and an |
62 |
| - out-of-place form `B = update_coefficients(A,u,p,t)`. |
63 |
| - 3. `isconstant(A)` trait for whether the operator is constant or not. |
64 |
| - |
65 |
| -### AbstractDiffEqLinearOperator Interface Description |
66 |
| - |
67 |
| - 1. `AbstractDiffEqLinearOperator <: AbstractDiffEqOperator` |
68 |
| - 2. Can absorb under multiplication by a scalar. In all algorithms, things like |
69 |
| - `dt*L` show up all the time, so the linear operator must be able to absorb |
70 |
| - such constants. |
71 |
| - 3. `isconstant(A)` trait for whether the operator is constant or not. |
72 |
| - 4. Optional: `diagonal`, `symmetric`, etc. traits from LinearMaps.jl. |
73 |
| - 5. Optional: `exp(A)`. Required for simple exponential integration. |
74 |
| - 6. Optional: `expv(A,u,t) = exp(t*A)*u` and `expv!(v,A::DiffEqOperator,u,t)` |
75 |
| - Required for sparse-saving exponential integration. |
76 |
| - 7. Optional: factorizations. `ldiv!`, `factorize` et al. This is only required |
77 |
| - for algorithms which use the factorization of the operator (Crank-Nicolson), |
78 |
| - and only for when the default linear solve is used. |
| 3 | +SciML has the [SciMLOpereators.jl](https://docs.sciml.ai/SciMLOperators/stable/) library for defining linear operators. |
| 4 | +The ODE solvers will specialize on this property, for example using the lienar operators automatically with Newton-Krylov |
| 5 | +methods for matrix-free Newton-Krylov optimzations, but also allows for methods that requires knowing that part of the |
| 6 | +equation is linear, like exponential integrators. See the SciMLOperators.jl library for more information. |
0 commit comments