Hello,
I’m trying to learn how to use the DiffEqOperators.jl module. It’s quite new. Chris Rackauckas gave a talk 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:
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:
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?
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.
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.
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.
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:
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
- Compute the arrays for v, K, and h.
- Multiply and divide to get the arrays for uv and u/h.
- 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.
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.
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:
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?
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.
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:
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:
while some_condition()
u += dt .* (∂*g.(u) + Δ*f.(u))
end
and doing this,
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]?
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.]