Update Section 2 and 3 for the overview document

[MSeeker1340]

Updates section 2 to reflect https://github.com/JuliaDiffEq/PDERoadmap/issues/29#issuecomment-412912592 .

I'm not sure about the naming however, in particular what "discretized differential operator" should refer to (L or A = LQ). I think calling L the "stencil operator" is appropriate, but there might be a more suitable name?

[MSeeker1340]

This does not actually concern the implementation directly but I have a question about the boundary conditions:

Can we actually solve Q from (5) and (6) using arbitrary B, b and R? For the most common boundary conditions we will certainly implement the respective Qs directly (or somehow incorporate into the derivative operators), but it would be a great plus if it can be automated generated.

[MSeeker1340]

I modified section 3.1 to more closely match section 2. In particular, one $Q$ is provided for each $B$, which I shall reference in the rest of section 3.

Also I used the name "stencil operators" in reference to $L_2$, $L^+_1$ and $L^-_1$, which of course is subject to change.

[jlperla]

I'm not sure about the naming however, in particular what "discretized differential operator" should refer to (L or A = LQ). I think calling L the "stencil operator" is appropriate, but there might be a more suitable name?

You are right, we should be careful with that and I like your naming. You are probably right that what people have in mind for the discretized one is the L Q product, which is square and enforces the boundaries. It is what you want to throw at solvers/etc.

Also don't forget to think through #28 sooner rather than later. If we allowed for affine operators from the beginning, then we would have it as A Q in general (instead of L because it could be affine). I assume that would stay affine and be reasonable to work with, based on what @ChrisRackauckas says about the algebra, but it makes a lot of changes to the document.

If so, then I think we need to have some notation so we can easily distinguish the operator with the imposed boundary condition. Maybe some sort of thing supersripted so it is A^{?} := A Q for some ?

[jlperla]

Can we actually solve Q from (5) and (6) using arbitrary B, b and R? For the most common boundary conditions we will certainly implement the respective Qs directly (or somehow incorporate into the derivative operators), but it would be a great plus if it can be automated generated.

This is a tricky problem, and something that @ChrisRackauckas and @dlfivefifty have made some comments on for algorithms. The one thing I would say is that I am pretty sure that B, b and R are not sufficient to create the Q. You need to know the (composed) stencil L as well of size M_bar X M When I finally am able to deliver on the examples it will help.

[MSeeker1340]

@jlperla Can you take a look at the updated section 3.2? The changes are more drastic than I anticipated.

In particular, the original version tries to solve the equation on the extended domain (24) and (25). Now that we agree to have each component operator handle their own $Q$, we don't need to worry about extending $rI$ anymore. So this becomes completely unnecessary and solving the equation on the interior nodes becomes a much more straightforward option.

[jlperla]

Now that we agree to have each component operator handle their own $Q$,

To be more clear there, I think @chrisrackauckas was less saying that they would handle their own directly and more that there may be simplifications when we apply the Q distributed.

My strong preference would be to keep the composed operator with padded zeros and the whole Q in the document. Then we can see how the application of the Q simplifies when distributed across the components of the operator. In that example, my guess is that the r I component of the operator, when padded with zeros and multiplied by the Q, ends up back at r I... But those things should be verified from the full L Q we know works (with the padded and composed L stencil of the operator)

[MSeeker1340]

My strong preference would be to keep the composed operator with padded zeros and the whole Q in the document. Then we can see how the application of the Q simplifies when distributed across the components of the operator.

I get your point. When we get to more complex examples (such as section 3.5), where we need to add together two differential operators, it makes sense to do what you're suggesting and compare the applicating Q globally vs distributed across the components. But I just don't see the necessity of verification in the simpler examples for the simple scaling operator rI.

Another angle: the intuition behind writing a discretized derivative operator as LQ is that Q handles the creation of the ghost nodes and L is just a plain stencil operation. This is simply not necessary for the scaling operator rI, as it is neither a derivative operator nor is in need of ghost nodes. Forcing it under the LQ framework is just too convoluted.

It's just my opinion nevertheless. If you think this is necessary then I will add back the verification section.

[MSeeker1340]

Updated section 2 to handle general, inhomogeneous bc with the corresponding B and Q operators. These should cover all section 3.

Edit: thinking about it I'm not actually sure what "absorbing/reflecting" boundary specifically refers to. I'm under the impression that they specifically refer to homogenous bc so I used Dirichlet/Neumann for the inhomogeneous ones, but this could be wrong.

[jlperla]

Forcing it under the LQ framework is just too convoluted.

Why have a special case for the r I? Better to do the algebra adding things up correctly with the padded zeros, and not have to assume that the Q drops out with arbitrary padding of zeros around the identity. The operation of composing the operators would need to add up all sorts of stuff. I want to be able to write simple linear algebra as checks on all this, and that involves adding up the identity with padded zeros.

Also, are you confident when we use more complicated boundary conditions that the padded zeros always drop out returning the identity for any Q? I am not.

[jlperla]

Edit: thinking about it I'm not actually sure what "absorbing/reflecting" boundary specifically refers to.

In this case the reflecting barrier is always homogenous u'(zval) = 0.

I would hold off on making too many changes to the examples. I will write more of them later. For people reading this solving HJBE, it is helpful to use the correct language for those processes to help them map their problems. Use dirichlet and Neuman for the boundary conditions themselves is helpful for some, but completely getting rid of reflecting and absorbing is not

[jlperla]

Another angle: the intuition behind writing a discretized derivative operator as LQ is that Q handles the creation of the ghost nodes and L is just a plain stencil operation.

No, I don't think that is correct at all. The L stencil is on the full extended space (after operator composition for the underlying L components. The boundary extrapolation operator Q makes sure that any vector defined on the discretized interior is mapped to a vector on the extended space which fulfills all the boundary conditions. So the Q is tightly connected to the exact boundary conditions on the extended space. Of course, if I use (L1 + L2) Q where some of the L1 and L2 were padded with zeros, then it could be that there is a submatrix of Q times the L1 before padding which delivers the same result of that multiplication. But that is something to be verified, and it may depend on the particular Q

[MSeeker1340]

I want to be able to write simple linear algebra as checks on all this, and that involves adding up the identity with padded zeros.

Sorry but I'm a bit at loss here. Why do we need to pad the identity operator with zeros in the first place? Can you give me an example?

I'd like to not think of rI as a special case, but instead the derivative operators (and the operators for jump diffusion) as special. It is because ghost nodes are needed that we make the distinction between the interior and the extended domain in the first place (and as a result, the Q operators). Regular operators that act on the interior directly and do not need the extended domain should, by this logic, not need an associated Q operator.

Let me make things clearer using the following stationary equation as an example:

$$ u''(x) + u'(x) + f(x)u(x) = g(x) $$

The discretized left hand side can be obviously composed by adding three parts: the discretized second-order derivative, the discretized first-order derivative and the scalar field $F$, a generalization of the scaling operator. The difference between the first two and the last one, is that $F$ does not require and boundary condition, and can be simply written as a diagonal matrix, mapping the interior points to the interior. I don't see any need of zero-padding, nor any need of verification.

Use dirichlet and Neuman for the boundary conditions themselves is helpful for some, but completely getting rid of reflecting and absorbing is not

They are still there. Basically I've made sure to refer the bc as Dirichlet/Nuemann when it can be inhomogenous, and absorbing/reflecting when we know it to be homogenous.

[ChrisRackauckas]

Yeah that's pretty much the point. You only need to make the ghost nodes that the operator actually uses.

[jlperla]

Let me make things clearer using the following stationary equation as an example:

The problem with your example is that you have not specified the boundary conditions. Q can only be calculated by looking at the boundary conditions on the composed operators. It is rectangular, depending on the maximum stencil size. It also cannot, as far as we know be determined without knowing the composed L operator... so I just don't understand what you are suggesting?

I think you are skipping some steps which may or may not apply for arbitrary Q . I think this comes down to your statement that "Now that we agree to have each component operator handle their own Q" which I did not agree to at all. I agreed that there were simplifications when you apply the Q distributed, which may make things easier so that the multiplication and ultimately not require padding zeros in practice, but that is different. My feeling is that you should forget about those simplifications for now, and we can verify them down the road and see how general they are at that point.

Following your example, if you have u′′(x)+u′(x)+f(x)u(x)=g(x) then L_tilde := D_xx + D_x + f(x) and the ODE is then L_tilde u(x) = g(x). That is leaving out boundary conditions, which are crucial for determining the Q matrix.

Lets say we discretize this stencil as L_tilde to L which we know has to be of size M X (M+2) by seeing that D_xx operator is second order. Whatever the affine boundary conditions would be, and they could be very general matrix when discretized of size (M+2) X 2 we know that the size of the Q is (M+2) X M. Then if we just naively add up matrices (which we can do as a check) then A = L * Q is M X M which is what we want to use for solving the A u = g equation as it embeds the boundary conditions.

Now lets think about the composition. First, if we were to just pad zeros, then lets define the components as L_2 discretizing the D_xx term, and of size (M+2) X M, the L_1 discretizing the D_x term, and of size (M+1) X M if we didn't padd it., and finally I as the constant term, of size M X M.

Obviously, we can't add these up to form the L directly, and padding with 0s is arbitrary (though convenient). If we padded with things other than 0s it would still work, but the Q matrix would change.

So after padding to make everything (M+2) X M lets call them L_2_pad, L_1_pad, and I_pad.

Note that L = L_2_pad + L_1_pad + diag(f(x)) * I_pad

and the whole discretization is, by linearity

A = L * Q = L_2_pad * Q + L_1_pad * Q + diag(f(x)) * Q * I_pad

here is where the simplifications come in. It is possible that a subset of Q, call it Q_1 of size(M+1) X Mgives usL_1_pdf Q = L_1 Q_1. and it is possible that a subset ofQ, call itQ_Iof sizeM X Mgives usQ I_pad = Q_I I`.

If so, then we have

A = L_2_pad * Q + L_1 * Q_1 + diag(f(x)) * Q_I * I

The hope in the implementation is that these subsets are easy enough to implement and strip out of the Q, but I don't think we can skip those steps in the math. We know that they would not be true for arbitrary paddings, for example. Once we have examples, we can see if the Q_1 and Q_I are trivial to determine... but I don't assume they are a priori.

[MSeeker1340]

The problem with your example is that you have not specified the boundary conditions. Q can only be calculated by looking at the boundary conditions on the composed operators. It is rectangular, depending on the maximum stencil size. It also cannot, as far as we know be determined without knowing the composed L operator... so I just don't understand what you are suggesting?

That's actually my point. We don't need to know the boundary condition to write out the discretized version of f(x), which is simply just diag(f(x_1),f(x_2),...,f(x_M)).

My feeling is that you should forget about those simplifications for now, and we can verify them down the road and see how general they are at that point.

Isn't having each component handling their own Q (or in the case of f(x), don't have a Q at all) more general than enforcing a single Q for the whole composed operator?

Following your example, if you have u′′(x)+u′(x)+f(x)u(x)=g(x) then L_tilde := D_xx + D_x + f(x) and the ODE is then L_tilde u(x) = g(x). That is leaving out boundary conditions, which are crucial for determining the Q matrix.

That's correct. Without the boundary condition/Q we can't define D_xx and D_x, but I thought the focus of the discussion is on f(x)?

[MSeeker1340]

Lets say we discretize this stencil as L_tilde to L...

Ok I see where our presumptions differ. If I understand this correctly, you seem to think we should treat the whole of L_tilde as a derivative operator, with a stencil and Q right? On the other hand, I'd like to think of it as

L_tilde = (D_xx + D_x) + f(x)

D_xx + D_x is a derivative operator and f(x) is just a regular matrix. Wouldn't you say this make things significantly simpler? So instead of

A = L * Q = L_2_pad * Q + L_1_pad * Q + diag(f(x)) * Q * I_pad

We have

A = A_diff + A_other = (L_2_pad + L_1_pad) * Q + diag(f(x))

Or more generally, if we use different Q for the derivative operators

A = L_2_pad * Q_2 + L_1_pad * Q_1 + diag(f(x))

[jlperla]

Ok I see where our presumptions differ.

But you are doing all of that without the boundary conditions, which might interact with the f(x). I think this is the crux of the issue. We can't discretize an operator which imposes boundary conditions on the whole equation without using the full boundary conditions.

Your way of separating out the non-derivative term as a special case may or may not work generally, but we do know that the algebra works without the special case... So let's start there and simplify if it turns out the special case works for a all affine boundary conditions. Which I am not sure is true.

[ChrisRackauckas]

Have you done the algebra yet?

[MSeeker1340]

After discussing this in detail with Jesse I agree to first leave the examples untouched and focus on getting the notations and derivations in the appendices right. Archiving this.