Design automatic sizing of the closure

[jlperla]

Right now, the Q matrix is being designed by scratch. The first step to automatically generating them is to making https://nbviewer.jupyter.org/github/MSeeker1340/PDERoadmap/blob/function_composition/Discretizing%20Linear%20Operators.ipynb#4.-Boundary-extrapolation-operator-$Q$ general is to realize that it does not always have a single point on either side, as part of the boundary.

There are a few main cases we will see:

1. Diffusion processes always add at least a point on both sides
2. FirstDifference(x, :backwards) always add a single point on the top of the interval, FirstDifference(x, :forwards) adds a point on the bottom.
3. Hence, Upwind ultimately could add a point on both sides, depending on the drift at the boundary. Might be reasonable to just assume it adds a point on both sides
4. Jump processes add multiple points, depending on how big the jumps are. You need to see how many points backwards it can jump at the bottom, and how many points forwards it can jump at the top of the interval.

My proposal is that their is a recursive function which goes through the operators to calculate the maximum number of points on both sides, where it takes the maximum.

Lets say that this returns $M_B^{\min}$ for the lower number of points, and $M_B^{\max}$ for the upper number of points. Then I think the `$Q_A$ of the $Q$ boundary extrapolation matrix is simply $$ QB = \begin{bmatrix} 0{1 \times MB^{\min}} \ I{M} \ 0_{1 \times M_B^{\max}} \end{bmatrix}' $$

[ChrisRackauckas]

Hence, Upwind ultimately could add a point on both sides, depending on the drift at the boundary. Might be reasonable to just assume it adds a point on both sides

Q is defined as extrapolating to the boundary, which is the "add a point on both sides" interpretation.

[jlperla]

But you agree that for FirstDifference it only adds a single point on one-side, depending on which direction you do the first-differences?

[ChrisRackauckas]

No, it's just a mapping from the closure to the interior. It just so happens that the first column is zero if it's forward differences, meaning that the boundary point on the left has no effect. (And this is where the consistency assumption has to come in)

[jlperla]

It works, but I don't think adding in the extra point is really necessary. You only need to extrapolate to the boundary where the differential operator actually grabs into the boundary?

More importantly, do you agree that the thickness of the boundary depends on the operator itself (which is clear for the jump-diffusion models)?

[ChrisRackauckas]

It works, but I don't think adding in the extra point is really necessary. You only need to extrapolate to the boundary where the differential operator actually grabs into the boundary?

It is necessary to make the implementation uniform. And when it's lazy it's not even an issue.

More importantly, do you agree that the thickness of the boundary depends on the operator itself (which is clear for the jump-diffusion models)?

No. Those operators are again from the closure to the interior. This is implicit in the assumptions necessary for consistency/uniqueness/existence.

[jlperla]

Then I do not understand how we are going to handle the absorbing and reflecting barrier for the jump-diffusions, if the "closure" is not thick. I thought that we had to add in extra rows to the Q and B matrices. You will need to better sketch out how the barriers are going to work for those. My understanding was from https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/260#issuecomment-368229781 and the local discussion to that.

[jlperla]

See https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/260#issuecomment-368207062 as well. Recall that the operators (prior to applying boundary conditions) are defined on the closure. The only way that you can have it defined on the closure and still deal with the sorts of boundary conditions above/below the barriers is to have the closure big enough to handle the largest jump above/below the interior.

[ChrisRackauckas]

Notice that all of the defined operators are from the closure to the interior. I was being lazy with how I wrote Q there. You can expand it out into multiple copies of the domain, and then restrict it to the real closure using the "v'd identity", i.e. diagonal 1's, then diagonal 1's doing up, repeat. The composition of those two operators is the extrapolation to (and beyond) the boundary in a closure -> interior operator.

[jlperla]

OK, I get what you are saying about all operations being from the closure to the interior. I think we will have to discuss this in more detail, as I am not sure what "v'd identity", means here and talk

But it sounds to me like you agree that if a poisson jump process can jump $3$ points above the top of a boundary, that the thickness of the closer adds 3 points? Put another way, let the interior be 0.0:1.0:10.0 and have poisson jumps with an arrival rate $\lambda$ jump 3.0 forward. Then I believe the infinitesimal generator of that stochastic process is something like $$ L u(x)\equiv \lambda (u(x + 3.0) - u(x)) $$

Then since the grid is $x=[0.0, 1.0, ... 10.0]$, the closure is $\bar{x} = [0.0, 1.0, ... 10.0, 11.0, 12.0, 13.0]$. Consequently, the $Q : R^{10} \to R^{13}$ ??? Just checking that we agree here on the dimensions.

[jlperla]

No. Those operators are again from the closure to the interior. This is implicit in the assumptions necessary for consistency/uniqueness/existence.

The other thing that I am wondering about with the negative drift example, $d X_t = -1.0 dt$ leading to the generator $\tilde{A} = -1.0 \partial_x$, with an absorbing barrier: by adding in an extra point on the top to make the "implementation uniform" (even if it does not enter the calculations, and is not strictly necessary), doesn' that mean that we have to add in second row of the $B$ matrix? If so, they aren't we effectively forcing in a second boundary condition for a 1st order ODE in the case that we solve simple ODEs with the operator? Maybe it is a never "binding" as a boundary condition, but is that kosher?

Of course, if there is a 2nd order diffusion term, then we really need to have that two boundary values. Sorry if I am missing something obvious and elementary in the differential equations here.

[jlperla]

@dlfivefifty Any thoughts on the case of a backwards drift described above? Is it kosher to add in a point (and accompanying row in the B matrix on it). Does this mean we are adding in an extra boundary condition?

[ChrisRackauckas]

If so, they aren't we effectively forcing in a second boundary condition for a 1st order ODE in the case that we solve simple ODEs with the operator?

There always is a second boundary condition which, if it satisfies the consistency requirements (i.e is compatible with the characteristics), is trivial.

[jlperla]

There always is a second boundary condition which, if it satisfies the consistency requirements (i.e is compatible with the characteristics), is trivial.

OK, great. Now back to the size of the boundary with the jump-diffusion. Thoughts on my statements in https://github.com/JuliaDiffEq/PDERoadmap/issues/11#issuecomment-385869573

[dlfivefifty]

Does this mean we are adding in an extra boundary condition?

One can decouple the number of boundary conditions from the order of the ODE: this is what is done in [Olver & Townsend 2012]. For example, if you happen to be at an eigenvalue it's reasonable to impose an extra boundary condition.

[dlfivefifty]

I don't think it was ever in a paper, but another simple example is x*u' + u = f, which is well-posed with zero boundary conditions, as the homogeneous solution is C/x.

[jlperla]

@dlfivefifty thanks so much, that is very helpful.

[jlperla]

I am closing. I think that it is becoming clear that extending every operator with 0s in order to fit the convex hull of the composed operator will make the Gaussian-elimination for the Q much cleaner. The (lazy) expansion with 0s may also make interior boundary conditions and patched together operators, like #20 work better.