SciML / DiffEqOperators.jl

Linear operators for discretizations of differential equations and scientific machine learning (SciML)
https://docs.sciml.ai/DiffEqOperators/stable/
Other
284 stars 74 forks source link

Operators defined on parts of the domain #59

Open jlperla opened 6 years ago

jlperla commented 6 years ago

Right now, my understanding of the domains and the boundary values is that if the operators coincide in the domain, we can take just add in as many extra ghost nodes outside of the domain in order to determine boundary values. I have a good sense of how that would work.

But just to make sure we are thinking through the boundary values and domains generally: There are a lot of problems in economics (including one of the ones I need to solve) where different regions of the state space are governed by different stochastic processes. When this is a choice and/or a free-boundary, the problem is written as a variational equality - which works perfectly fine since the underlying functions are still defined on the entire space.

But what if the differential operators are defined on fixed but connected regions of the domain, with boundary values determining the interface? For example, L_1 = mu_1(x) d/dx + sigma^2/2 d^2/dx^2 defined on x = [0, .5] and L_2 = mu_2(x) d/dx + sigma^2/2 d^2/dx^2 on x = [0.5,1].

Normally, you might do this with two functions and boundary values, for a function v_1(x) defined on [0,0.5] and v_2(x) on [0.5,1]:

How would we do this sort of setup in the general setup? Do we separate two functions or is this a single function with "interior" boundary values?

ChrisRackauckas commented 6 years ago

This is more like FEM where you don't necessarily have the boundaries on the ends. If you write it out in operator form it's clear. v is the interior, which is [0,0.5) union (0.5,1.0]. So then u = Q*v gives the full space, and you can work out what Q is from the BCs. The other 0.0 and 1.0 are just like before. For the point at 0.5 you get two linear relations and a two constraints. To first order, you put them together to get u[mid] - u[mid-1] = u[mid] - u[mid+1] or u[mid] = 0.5*(u[mid-1] + u[mid+1]), and so Q is I along with those three BC calculations. Then the operators are well defined on [0,0.5] -> (0,0.5) and [0.5,1.0] -> (0.5,1.0) and when you piece that back together you get v_t.

jlperla commented 6 years ago

OK. Lets keep this example in our minds down the road when we think about automatically expanding the domains when we compose operators. For example, if L_1 is defined above (and we associate the domain [0, 0.5] witth it, and L_2 with [0.5, 1] then it seems like the notation for a typical composition of the operators L = L_1 + L_2 should be possible . When we do the discretizaiton of the operator, I imagine there would simply be 0's in the interior of (0.5, 1] for L_1 and [0, 0.5) for L_2.

The boundary values could be done manually, as always, for now.

ChrisRackauckas commented 6 years ago

The boundary values could be done manually, as always, for now.

Yeah, we'll need to find a good way to relate the (Q,R) to a mesh and things like that, but what we can do here is make the generation of the A stencils easy, and do the standard (Q,R) operators.

jlperla commented 6 years ago

Yes. It also tells me that we need to be able to associate the domain with individual operators so that we can stitch together the union of them for the domain when we lazily compose.

ChrisRackauckas commented 6 years ago

Well technically this operator still acts on the whole interior, it's just with a constant of zero some places.

But wait a second, why do you need mu_1 and mu_2? If it's dependent on x, isn't this just one spatially-dependent coefficient?

jlperla commented 6 years ago

Yes, the operator acts on the whole space, but we need the domain to know which ones we should zero out

This was intended as a MWE. You are right that in this case we could just have a mu(x) on its own, but there might be much more complicated dynamics (and even different stochastic processes) in different regions.