Update the overview and appendix

[MSeeker1340]

@jlperla Opening this up for discussion. The overview section can be updated alongside the appendix, and after this is done we can move on to update the examples.

I think adding an overview section to appendix A (Gaussian Elimination and Q) is very much needed because it is not evident what it tries to achieve. I'll post my summary here so please point out anything that is wrong:

• The main text and the examples focuses on constructing $LQ$ derivative operators that act on the interior. The equation to solve is $LQu = p$. The $B$ operator does not directly come into the equation.

• The appendix, on the other hand, tries to look at the discretized problem in another angle. Instead of solving the equation on the interior, the equation to solve is

  $$ L\bar{u} = p\quad B\bar{u} = b $$

  Here $B$ plays a pivotal role and $Q$ does not show itself. And we want to show via examples that these two approaches are equivalent. In particular, the equation for the interior can be derived from the equation on the extended domain by stacking the two equations and using Gaussian elimination.

I also feel like the relationship between Gaussian elimination and the Q operator is not explained clearly in A.1.4 and A.2.1. We're still pulling out Q from thin air and just checking the $LQ$ product gives the same matrix as Gaussian elimination. I haven't done the algebra fully but I feel what should be done is this: we express the elimination steps as well as dimension reduction as applying elementary matrices. This should be able to reproduce the Q matrix somehow, and then the equivalence can be adequately explained.

[jlperla]

I am trying to remember how far I got on the gaussian elimination, but I don't think I made it all the way to showing where the Q came from... It has something to do with the product of the elementary matrices performing the gaussian elimination of the B, but I don't think I cracked it. If the section isn't helpful then leave it in an appendix and we can complete later... But if you want to take a little time to see if you can figure out how the gaussian elimination elementary matrices to eliminate the B rows of the stacked system are connected to the Q^B, I think it is helpful.

[MSeeker1340]

Updated Appendix A with an overview and updated examples. I have made many references to section 2 but it is not updated yet, so it may not seem right (for example, the L1 operators in section should actually be their extended version).

I think we should add an extra example to make things complete, which uses inhomogeneous boundary condition, and prove the equivalence between the extended equation and $LQu = p$ for affine $Q$s. Modifying the first example should be enough, as the same elimination steps can be used.

[MSeeker1340]

Also, in the original A.1.3 there's talk about the restriction operator when the matrix is singular. This seems out of place for me and doesn't really make sense, because if the matrix is singular how come we can define an inverse in (71)?

[ChrisRackauckas]

It doesn't make sense to talk about it being singular since it's not square.

[jlperla]

I am not sure exactly where it is in the document, but recall that the gaussian elimination steps start with the fully stacked system with the L and the B both in a square system. It could have been referring to that, or it could have been useless

[MSeeker1340]

I am not sure exactly where it is in the document, but recall that the gaussian elimination steps start with the fully stacked system with the L and the B both in a square system. It could have been referring to that, or it could have been useless

Yes I'm pretty sure the singularity refers to the stacked system [L;B], which is square.

I think it is a legitimate concern, but as I pointed out the original discussion doesn't make much such because the inverse of [L;B] is used anyways, which is only possible if it's nonsingular.

[ChrisRackauckas]

That's non-singular iff Q is linear.

[jlperla]

Love it. Don't lose these propositions, even if we don't prove them all at first.

[ChrisRackauckas]

Yeah I remember looking at it and then under this setup (linear BCs, etc.) the partial Gaussian elimination allows the solution to be in a family of solutions up to an additive constant. That additive constant can be determined by the BCs though by looking at the trivial case, and this is why you can/must have affine operators in some cases.

[MSeeker1340]

Yeah I remember looking at it and then under this setup (linear BCs, etc.) the partial Gaussian elimination allows the solution to be in a family of solutions up to an additive constant.

I can definitely see the possibility of this. On the other hand the fact that the equation is singular should not affect what we're trying to do in the appendix, i.e. proving that the extended equation using L and B is the same as the equation on the interior, which is LQu = p. We're only using Gaussian elimination to get the same equation, and whether the equation has an unique solution is a property of the problem itself.

So I think this part can be safely removed, but for the sake of completeness we can add a note in the overview section of Appendix A saying that the equation may be singular.

[MSeeker1340]

@jlperla Can you take a look at the updated section 2, in particular the new section regarding operator composition?

For the extension operator, the original one appearing in Appendix A at (65) is obviously wrong and we should have E_11 be applied to the right. Right multiplication also naturally results in the two equivalent interpretations, i.e. extension for the stencil/restriction for the domain, which I found helpful.

If you think this is ok, I will go back to the appendix and reflect the changes. I also need to update the reference operators in section 3.1, after which the document update can be finished except for the examples.

[MSeeker1340]

Forgot to mention a small nitpick: I don't think the continuous boundary condition (2) is correct, as B_tilde*u_tilde(x) should result in a vector of boundary values instead of some field b_tilde(x). Now that we're writing the discrete B as affine, perhaps we should also do the same and rewrite (2) as B_tilde*u_tilde(x) = 0?

[jlperla]

ow that we're writing the discrete B as affine, perhaps we should also do the same and rewrite (2) as B_tilde*u_tilde(x) = 0?

I agree, it would be much more consistent, and probable closer to how it would be implemented in the higher level interface. We just need to be a little careful with examples for people not used to working with affine algebra and functions (e.g. me) to take that extra step making from a typical non-homogeneous boundary condition to the affine version.

@jlperla Can you take a look at the updated section 2, in particular the new section regarding operator composition? If you think this is ok, I will go back to the appendix and reflect the changes.

Looks great as far as I can tell, and exactly what we were talking about for the difference between the composition after using the distributed law. You leave it open how the B_k is derived from B and how Q^B_k is derived from Q^B, which I think is the right thing to do. Examples will help a lot there. Perhaps @ChrisRackauckas can do a sanity check. I think that Appendix B on the affine relationships is very important to get right (as it was part of the original discussion with Chris). I am not sure, though, whether it belongs in a separate appendix or merged into sections 1 and 2. I defer to you.

For the extension operator, the original one appearing in Appendix A at (65) is obviously wrong and we should have E_11 be applied to the right. Right multiplication also naturally results in the two equivalent interpretations, i.e. extension for the stencil/restriction for the domain, which I found helpful.

Not sure exactly what you are referrring to, but I trust your judgement and suspected lots of mistakes are in there.

[ChrisRackauckas]

Can you post a current build of the PDF?

[MSeeker1340]

@ChrisRackauckas Here:

linear_operators_overview.pdf

The focus is on the overview as well as appendix A, but I still need to update the outdated references.

[ChrisRackauckas]

While we could add in arbitrary points in the extension, the algebra will be easier if we add 0s.

It's not an arbitrary choice. It's because the operator is local and the PDE is only defined on the interior, so there must be 0 influence off the diagonal in that discretization in order for it to be stable.

[ChrisRackauckas]

On the other hand, for the reacting boundaries, we know the boundary extrapolation operator is

You might want to demonstrate this from the definition of the BC, solving for the boundary values.

[ChrisRackauckas]

85 is dimensionally incorrect (operators on one side, vector on the other).

Therefore, we actually want to solve u = AQu.

no, u' = AQu, not ubar on the left side. A maps back to the interior. 86 is wrong in the same way.

I'm getting the impression B hasn't been touched yet. Part A is good.

[MSeeker1340]

85 is dimensionally incorrect (operators on one side, vector on the other). Therefore, we actually want to solve �u = AQu. no, u' = AQu, not ubar on the left side. A maps back to the interior. 86 is wrong in the same way.

Yeah I haven't touched appendix B yet. The first example Jesse wants me to implement (https://github.com/JuliaDiffEq/PDERoadmap/issues/33) does not use affine operators so I want to leave this for later.

I'm getting the impression B hasn't been touched yet.

What do you mean by this?

[ChrisRackauckas]

What do you mean by this?

Yeah I haven't touched appendix B yet.

It looked like you didn't fix Appendix B yet.

[MSeeker1340]

Oh I interpreted it as B the boundary operator ;)

[ChrisRackauckas]

Thanks for writing these examples out. It is a very clear demonstration of the general phenomena. Okay is there anything left here? Does this put us on the same page?

[MSeeker1340]

I think this should be enough, though it's only a theoretical framework and we probably need more discussions about the implementation/interface, especially for irregular grid/domain (https://github.com/JuliaDiffEq/PDERoadmap/issues/4).

@jlperla I'll go ahead and update the outdated references in appendix A & operators in section 3.1, after which I think the document should be in a good state.

[MSeeker1340]

@jlperla You can merge this now if everything is ok. To summarize, section 2 and 3.1 is updated to use affine B and Q^B notation and a section about composing stencils/zero padding is added to section 2. Appendix A is also updated to reflect the changes along with an additional example. The examples in section 3 and appendix B are untouched for now.

[jlperla]

I defer to you guys to merge it. No need to have me involved in some of those details now that we are on the same page

[MSeeker1340]

Oh I don't have write access to PDERoadmap yet.

[jlperla]

@chrisrackauckas Can you add him?

[ChrisRackauckas]

There just did.