clean up affine operator

[vpuri3]

• define traits
• add tests
• mul!, *

[codecov]

Codecov Report

Merging #21 (36047b2) into master (aa02552) will not change coverage.
The diff coverage is 0.00%.

@@          Coverage Diff           @@
##           master     #21   +/-   ##
======================================
  Coverage    0.00%   0.00%           
======================================
  Files           4       4           
  Lines         459     498   +39     
======================================
- Misses        459     498   +39     

Impacted Files	Coverage Δ
src/sciml.jl	0.00% <0.00%> (ø)
src/basic.jl	0.00% <0.00%> (ø)
src/interface.jl	0.00% <0.00%> (ø)

📣 Codecov can now indicate which changes are the most critical in Pull Requests. Learn more

[vpuri3]

what is the matrix representation of

L = AffineOperator(A, b)
L(u) = A*u + b

we would also need a new name. affine is too vague. something that signifies that it is adding a vector, like vecsumoperator or muladdoperator or something

[vpuri3]

im going with MuladdOperator for now

[vpuri3]

thoughts @ChrisRackauckas ?

[ChrisRackauckas]

what is the matrix representation of

You can dig up some of the old conversations in DiffEqOperators.jl about this question (SciML/DifferentialEquations.jl#260 (comment), SciML/DiffEqOperators.jl#53). There is no matrix representation of it, since if there was it would be a linear operator instead of an affine operator. You can only represent an affine operator as a linear operator in a dimension of one larger. Basically, the operation is: [A b] * [u;1], so it requires something modified to the input as well.

While it this seems like it might doom the idea of using matrix-free affine operators, it turns out that it's okay because of an obscure fact that iterative linear solvers are actually compatible with affine operators. Everyone talks about and proves them for linear solvers, but that's actually over-restrictive.

You can show some nice properties of affine operators, like that the definition given of * will make Qx=b in an iterative solver converge to Qax = b-Qb (so if we actually define a type that does this in Julia, it will work with iterative solvers without a fuss. Just the same thing idea as writing in residual form).

I forget where I wrote the proof of that, some PDE solvers roadmap. But it's rather straightforward to prove. If Q is the affine operator Q(x) = Ax + c, then solving Qx = b is equivalent to solving Ax + c = b or Ax = b-c. You do the same thing of plugging the affine operator into the GMRES/CG/etc. proof, just move the affine part to the rhs residual, and show it converges to solving Ax = b-c, and so GMRES/CG/etc. solves Q(x) = c for an affine operator properly. That same trick then can be used pretty much anywhere you would've had a linear operator to extend the proof to affine operators, so then exp(A*t)*v operations via Krylov methods work for A being affine as well, and all sorts of things. So affine operators have no matrix representation but they are still compatible with pretty much any iterative or Krylov method.

And many boundary conditions give not linear operators but affine operators, hence the existence of all of this.