Linear problems

[DanielVandH]

As mentioned in #16, it is possible to make the code faster in a specific setting by exploiting linearity to define a sparse linear system. Namely, when $\boldsymbol{q}(x, y, t, u) = A_q(x, y, t)\boldsymbol{\nabla} u + \boldsymbol{B}_q(x, y, t) + \boldsymbol{C}_q(x, y, t)u$, $R(x, y, t, u) = A_R(x, y, t) + B_R(x, y, t)u$, and the boundary conditions are also linear. In this setting, we can write the PDE in the form

$$ \dfrac{\mathrm d\boldsymbol{u}}{\mathrm dt} = \boldsymbol A\boldsymbol u + \boldsymbol b $$

for some matrix $\boldsymbol A$ and vector $\boldsymbol b$, and $\boldsymbol A$ is sparse. This could be worth implementing.

One issue would be thinking about the interface for this that doesn't disrupt the existing interface. Maybe an interface like

struct LinearFlux{F1, F2, F3} 
  Aq::F1
  Bq::F2
  Cq::F3
end
(f::LinearFlux)(x, y, t, a, b, c, p) = Aq(x, y, t, p) * [a, b] + Bq(x, y, t, p) + Cq(x, y, t, p) * (a*x + b*y + c) 
# (with obvious extensions for in-place methods)

struct LinearReaction{F1, F2}
  AR::F1
  BR::F2
end
(f::LinearReaction)(x, y, t, u, p) = AR(x, y, t, p) + BR(x, y, t, p) * u

# general form needed for the boundary condition functions

The functions Aq, Bq, Cq, AR, and BR would be FunctionWrappersWrappers, although the user can just define anonymous functions that are then wrapped (as is done already for the boundary functions). The user would then provide these structs so that the solver constructs the matrices $\boldsymbol A$ and vector $\boldsymbol b$. The number of non-zeros of $\boldsymbol A$ is already known, so pre-allocating (as a dualcache) this matrix is easy, simply storing the coordinate vectors (I, J, V) that define a sparse matrix. May be some difficulty handling the boundary conditions in this case.

The same methods used for computing $\boldsymbol A\boldsymbol u + \boldsymbol b$ would be used for solving the steady state case, where $\boldsymbol A\boldsymbol u + \boldsymbol b = \boldsymbol 0$ so that $\boldsymbol u = -\boldsymbol A^{-1}\boldsymbol b$.

Other notes:

• If $\boldsymbol{B}_q = A_R = 0$, then $\mathrm d\boldsymbol u/\mathrm dt = \boldsymbol A\boldsymbol u$.
• If the ODE is autonomous, $\boldsymbol A$ and $\boldsymbol b$ need only be computed once.
• If $\mathrm d\boldsymbol u/\mathrm dt = \boldsymbol A\boldsymbol u + \boldsymbol b(t)$, then $\boldsymbol u = \mathrm{e}^{\boldsymbol At} \left[\mathrm{e}^{-\boldsymbol At_0}\boldsymbol u(t_0) + \int_{t_0}^t \mathrm{e}^{-\boldsymbol As}\boldsymbol b(s)~\mathrm ds\right]$. Could be worth working with the matrix exponential here? Integral looks a bit troublesome though.

Actually, https://docs.sciml.ai/DiffEqDocs/stable/types/nonautonomous_linear_ode/ and https://docs.sciml.ai/DiffEqDocs/stable/features/diffeq_operator/#DiffEqOperators would probably be useful here.

[DanielVandH]

I'm not so sure this is all worth the effort anymore. A PR for this would be great, but it will not be me who does it.