# Trouble with boundary conditions in DiffEqOperators.jl

**URL:** <https://discourse.julialang.org/t/trouble-with-boundary-conditions-in-diffeqoperators-jl/21488>\
**Category:** Numerics\
**Tags:** diffeq\
**Created:** [March 5, 2019, 9:05am UTC](https://discourse.julialang.org/t/trouble-with-boundary-conditions-in-diffeqoperators-jl/21488 "2019-03-05T09:05:46Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![danielc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielc/32/220_2.png) [@danielc](https://discourse.julialang.org/u/danielc)\
**Post date:** [March 5, 2019, 9:05am UTC](https://discourse.julialang.org/t/trouble-with-boundary-conditions-in-diffeqoperators-jl/21488/1 "2019-03-05T09:05:46Z")

</div>

Hello,

I’m trying to learn how to use the [DiffEqOperators.jl](https://github.com/JuliaDiffEq/DiffEqOperators.jl) module. It’s quite new. Chris Rackauckas [gave a talk](https://www.youtube.com/watch?v=okGybBmihOE) at JuliaCon 2018. To get started, I tried to choose a very simple problem:

- Compute the second derivative of sin(x).

The module provides `DerivativeOperator` and you can use it like this:

```julia
using DiffEqOperators

deriv = 2 # 2nd order derivative.
order = 2 # 2nd order approximation.
N = 100 # Number of grid points.
xs = range(0,stop=2π,length=N)
dx = xs[2] - xs[1] 
#
# Use Delta to denote the Laplace operator.
#
Δ = DerivativeOperator{Float64}(deriv,order,dx,N,:Dirichlet0,:Dirichlet0)

```

The last two parameters are the boundary conditions. In this example I have Dirichlet boundary conditions (i.e. y(0) = y(2\pi) = 0).

Then I tried a simple plot:

```julia
using PyPlot

ys = sin.(xs)

figure(1)
plot(xs,ys,"C0-") # sin(x)
plot(xs,Δ*ys,"C3-") # Δ(sin(x)) = - sin(x)

legend([
	L"\sin(x)",
	L"\frac{\partial^2}{\partial x^2} \sin(x) = - \sin(x)"
])

```

Here is the result. As you can see, the derivative operator did a pretty good job at making the Laplacian of sin(x), except at the very ends. Literally the first and points. At the first and last values, the derivative jumps to +15.7 and -15.7; so I imagine that there must be something wrong with the way I setup the boundary conditions. I tried `:Neumann0` boundary conditions too but it didn’t help at all. Does anyone have any suggestions?

![example](https://global.discourse-cdn.com/julialang/original/3X/8/f/8f3d28b2896c857e8646b01d6814011a1e01456f.png)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 5, 2019, 10:15am UTC](https://discourse.julialang.org/t/trouble-with-boundary-conditions-in-diffeqoperators-jl/21488/2 "2019-03-05T10:15:35Z")

</div>

Yes, the boundary condition issues are known, which is why docs haven’t been added and it hasn’t been “released”. Someone needs to find the time to address it. It would be a great GSoC project.

---

<div class="post-metadata">

**Author:** ![danielc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielc/32/220_2.png) [@danielc](https://discourse.julialang.org/u/danielc)\
**Post date:** [March 5, 2019, 10:24am UTC](https://discourse.julialang.org/t/trouble-with-boundary-conditions-in-diffeqoperators-jl/21488/3 "2019-03-05T10:24:47Z")

</div>

Thanks for the response. Is there a workaround that you can suggest? To emulate `:Dirichlet0` could I just sent the first and last values to zero after each iteration? To emulate `:Neumann0` could I just compute the finite difference of suitable order by hand at the first and last cell points for every iteration? … It seems to me that this should work, but I’m not sure if I’ll just end up breaking something I don’t expect.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 5, 2019, 11:53am UTC](https://discourse.julialang.org/t/trouble-with-boundary-conditions-in-diffeqoperators-jl/21488/4 "2019-03-05T11:53:06Z")

</div>

Well the linear operators, Dirichlet0 and Neumann0, work. But the PDEs are defined in the interior of the domain, so you’d want to plot `[0;Δ*ys;0]` with `ys` defined only on the interior. This of course is quite an odd object to handle, which is the issue.

---

<div class="post-metadata">

**Author:** ![danielc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielc/32/220_2.png) [@danielc](https://discourse.julialang.org/u/danielc)\
**Post date:** [March 5, 2019, 6:33pm UTC](https://discourse.julialang.org/t/trouble-with-boundary-conditions-in-diffeqoperators-jl/21488/5 "2019-03-05T18:33:20Z")

</div>

I’m not 100% sure I understand what you’re saying. But if the PDE is just not defined at the boundary, that sounds like I can make my grid one cell larger at each end, call that the boundary, and throw away those cells later:

```julia
N = 100 # Number of grid points.
xs = range(0,stop=2π,length=N)
dx = xs[2] - xs[1] 

# Add a throw-away grid cell on either side.
xs = [xs[1]-dx ; xs ; xs[end]+dx]
Δ = DerivativeOperator{Float64}(deriv,order,dx,N+2,:Dirichlet0,:Dirichlet0)
ys = sin.(xs)

plot(xs[2:end-1], ( ys )[2:end-1], "C0-")
plot(xs[2:end-1], (Δ*ys)[2:end-1], "C3-")

```

The real problem that I want to solve is that I want to evolve a 1D diffusion equation that looks roughly like this:

\frac{\partial u}{\partial t} = \frac{\partial}{\partial x} \left[uv - K \frac{\partial}{\partial x} \left( \frac{u}{h} \right) \right]

None of the terms are constant; they all vary with x and t. My plan was to expand that into a set of first and second order derivatives,

\frac{\partial u}{\partial t} = \frac{\partial}{\partial x} \left( uv \right) - \frac{\partial}{\partial x} \left( K \right) \frac{\partial}{\partial x} \left( \frac{u}{h} \right) - K \frac{\partial^2}{\partial x^2} \left( \frac{u}{h} \right)

then on each timestep I would

1. Compute the arrays for v, K, and h.
2. Multiply and divide to get the arrays for uv and u/h.
3. Use `DiffEqOperators` to compute the next value of u.

It looks to me like I just need to make two operators at the beginning of my program: one for \partial / \partial x and one for \partial^2 / \partial x^2. Then I can just apply them in step (3) for every iteration.

Do you think this plan can still work? Can I just make my grid one cell larger on either side? Do I need to worry about the wrong values at the boundary “propagating” into the region I care about?

Thanks for the help.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 5, 2019, 8:41pm UTC](https://discourse.julialang.org/t/trouble-with-boundary-conditions-in-diffeqoperators-jl/21488/6 "2019-03-05T20:41:20Z")

</div>

Yes, but the ODE is defined on the interior. The boundary terms are defined by the interior and the boundary condition, so they don’t show up in the ODE resulting from the spatial discretization.

---

<div class="post-metadata">

**Author:** ![danielc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielc/32/220_2.png) [@danielc](https://discourse.julialang.org/u/danielc)\
**Post date:** [March 5, 2019, 9:11pm UTC](https://discourse.julialang.org/t/trouble-with-boundary-conditions-in-diffeqoperators-jl/21488/7 "2019-03-05T21:11:27Z")

</div>

I’m not trying to be obtuse, I promise. When you say that the boundary terms don’t show up in the ODE resulting from the spatial discretization, are you saying that I can literally not worry about the boundary? For example, can I do this:

```julia
x = [x[1]-dx ; x ; x[end]+dx]
u = foo.(x)
∂ = DerivativeOperator{Float64}(1,1,dx,N+2,:Dirichlet0,:Dirichlet0)
Δ = DerivativeOperator{Float64}(2,2,dx,N+2,:Dirichlet0,:Dirichlet0)

for t in t0:dt:tmax
	u += dt .* (∂*g.(u,x,t) .+ Δ*f.(u,x,t))
end

```

Notice that there’s no effort being made here to fix the bad values at the boundaries. Can I carelessly run that for loop and trust that `u[2:end-1]` will be fine at the end?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 5, 2019, 10:54pm UTC](https://discourse.julialang.org/t/trouble-with-boundary-conditions-in-diffeqoperators-jl/21488/8 "2019-03-05T22:54:22Z")

</div>

> [@danielc](#):
>
> Can I carelessly run that for loop and trust that `u[2:end-1]` will be fine at the end?

Well it’s not carelessly. The choice of a boundary condition discretization defines the boundary as a function of the interior. If you also define the boundary in the array, you have more values than degrees of freedom.

---

<div class="post-metadata">

**Author:** ![danielc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielc/32/220_2.png) [@danielc](https://discourse.julialang.org/u/danielc)\
**Post date:** [March 5, 2019, 11:43pm UTC](https://discourse.julialang.org/t/trouble-with-boundary-conditions-in-diffeqoperators-jl/21488/9 "2019-03-05T23:43:51Z")

</div>

> Well it’s not carelessly. The choice of a boundary condition discretization defines the boundary as a function of the interior. If you also define the boundary in the array, you have more values than degrees of freedom.

I think that answers the other question I wanted to ask but was not sure how to formulate. I very much care about what happens at the boundary. I have one problem where I don’t want any material “escaping” the system, and another where I want it to escape at one end but not the other.

This is what I had in mind:

```julia
ClosedTube = DerivativeOperator{T}(...,:Neuman0,:Neuman0)
OpenTube = DerivativeOperator{T}(...,:Dirichlet0,:Neuman0)

```

So I was worried that those bad values at the boundary cells meant that the behaviour at the boundaries would be wrong. Could you confirm that,

- if I choose the correct boundary conditions,
- and I add the two throw-away cells at either end…

then the behaviour at the boundaries will be correct?

I’m still uncertain about the implications evaluating my functions at two end points with meaningless values and then iterating. I have to imagine that there must be a difference between doing this:

```julia
while some_condition()
	u += dt .* (∂*g.(u) + Δ*f.(u))
end

```

and doing this,

```julia
while some_condition()
	du = ∂*g.(u) + Δ*f.(u)
	du[1] = du[end] = 0
	u += dt * du
end

```

Does it really not make any difference inside `u[2:end-1]`?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 6, 2019, 12:12am UTC](https://discourse.julialang.org/t/trouble-with-boundary-conditions-in-diffeqoperators-jl/21488/10 "2019-03-06T00:12:14Z")

</div>

The discretized Laplacian is an N x N+2 stencil. To get an NxN stencil, a truncation has to occur which essentially embeds the boundary conditions into the stencil, making it an affine operator if non-zero BC, and a linear NxN if Neumann or Dirichlet with 0. You can then write out the boundary values using the condition. For example, if Dirichlet it’s trivial, since `[bc(L);u;bc(R)]` would be the N+2 points. Likewise you can do the same for Neumann, where a polynomial interpolation gives you a linear equation for the value at the boundary. Neumann0 with the simplest discretization of u’ implies that `[u[1];u;u[end]` has to hold.

So if you think of this extrapolation operator Q:N → Nx2 as the boundary conditions just written in the inverted form (pseudo-inverse in the space of functions that satisfy the BCs), then the operation is actually AQu. Since it’s much easier for people to think square, we give `Δ = AQ` on the interior, at least in the case of Neumann0 and Dirichlet0, and that defines the ODE on the interior while the boundary conditions are completely determined by these values and the chosen discretization (assuming linear BCs, like Robin).

[When non-zero, we don’t currently add the affine term which is the issue we need to address.]