# Dirac delta in boundary conditions for MethodOfLines.jl

**URL:** <https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567>\
**Category:** Modelling & Simulations\
**Created:** [May 5, 2022, 8:11pm UTC](https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567 "2022-05-05T20:11:04Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [May 5, 2022, 8:11pm UTC](https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567/1 "2022-05-05T20:11:04Z")

</div>

As a basis for more complex models, I’m trying to get a simplified version of the [Von Foerster equation - Wikipedia](https://en.wikipedia.org/wiki/Von_Foerster_equation) discretized and running in MethodOfLines.jl.

Here’s my model:

```julia
using ModelingToolkit
using MethodOfLines
using OrdinaryDiffEq
using DomainSets

@parameters t a μ
@variables v(..) u(..)
Dt = Differential(t)
Da = Differential(a)
p = [μ => 0.04]
tmax = 100.0
amax = 100.0
da = 0.1
dt = 0.1

pde_eqns = [Dt(u(t,a)) + Da(u(t,a)) ~ -μ*u(t,a)]
domains = [t ∈ Interval(0,tmax), a ∈ Interval(0,amax)]
function a0(a)
    if a > 0
        return 0.0
    else
        return 1.0
    end
end
@register a0(a)
bcs = [u(t,0) ~ 0.0, u(0,a) ~ a0(a)]
@named pde_sys = PDESystem(pde_eqns, bcs, domains, [t,a], [u(t,a)], p)
discretization = MOLFiniteDifference([a=>da], t, approx_order = 2, grid_align = center_align)
pde_prob = discretize(pde_sys, discretization)
pde_sol = solve(pde_prob, Tsit5(), saveat = dt)

```

The boundary conditions for u(0,a) are meant to be 1 when a=0 and 0 when a\>0, which is what I’m trying to do with the registered function. However, when I look at the discretized problem, the initial conditions are 0 for all values of a:

```julia
julia> pde_prob.u0
1000-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.0
 0.0
 0.0
 0.0
 0.0

```

I think I understand what’s going on here - while the above function a0 is correct in the continuous sense, it doesn’t work in the discretized system. If I change function:

```julia
function a0(a)
    if a > da
        return 0.0
    else
        return 1.0
    end
end

```

This now does something more sensible:

```julia
julia> pde_prob.u0
1000-element Vector{Float64}:
 1.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.0
 0.0
 0.0
 0.0
 0.0

```

However, now the solution gives nonsensical results - the integral over a should be an exponential decay, but it shows exponential growth, which again points to another boundary issue problem. Any tips on how to fix this? Full code, along with the ‘correct’ exponential decay as an ODE is given below.

```julia
using ModelingToolkit
using MethodOfLines
using OrdinaryDiffEq
using DomainSets

@parameters t a μ
@variables v(..) u(..)
Dt = Differential(t)
Da = Differential(a)
p = [μ => 0.04]
tmax = 100.0
amax = 100.0
da = 0.1
dt = 0.1

# ODE version
ode_eqns = [Dt(v(t)) ~ -μ*v(t)]
v0 = [v(t) => 1.0]
@named ode_sys = ODESystem(ode_eqns)
ode_prob = ODEProblem(ode_sys, v0, (0,tmax), p)
ode_sol = solve(ode_prob, Tsit5(), saveat = dt)

# PDE version
pde_eqns = [Dt(u(t,a)) + Da(u(t,a)) ~ -μ*u(t,a)]
domains = [t ∈ Interval(0,tmax), a ∈ Interval(0,amax)]
function a0(a)
    if a > da
        return 0.0
    else
        return 1.0
    end
end
@register a0(a)
bcs = [u(t,0) ~ 0.0, u(0,a) ~ a0(a)]
@named pde_sys = PDESystem(pde_eqns, bcs, domains, [t,a], [u(t,a)], p)
discretization = MOLFiniteDifference([a=>da], t, approx_order = 2, grid_align = center_align)
pde_prob = discretize(pde_sys, discretization)
pde_sol = solve(pde_prob, Tsit5(), saveat = dt)

```

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [May 10, 2022, 1:34pm UTC](https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567/2 "2022-05-10T13:34:48Z")

</div>

I don’t know what’s happening here, but a suggestion:  
Try

```julia
using IfElse
a0(a) = IfElse.ifelse(a>0, 0.0, 1.0)

```

Let me know if that works, if not I can take a closer look.

---

<div class="post-metadata">

**Author:** ![skleinbo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skleinbo/32/36080_2.png) [@skleinbo](https://discourse.julialang.org/u/skleinbo)\
**Post date:** [May 10, 2022, 1:49pm UTC](https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567/3 "2022-05-10T13:49:24Z")

</div>

`bcs = [u(t,0) ~ 0.0, u(0,a) ~ a0(a)]`

Those seem contradictory to me 🤔

---

<div class="post-metadata">

**Author:** ![skleinbo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skleinbo/32/36080_2.png) [@skleinbo](https://discourse.julialang.org/u/skleinbo)\
**Post date:** [May 11, 2022, 11:32am UTC](https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567/4 "2022-05-11T11:32:47Z")

</div>

Also, a centered difference scheme is probably unstable for an advection equation. You should rather be using e.g. an [upwind scheme](https://en.wikipedia.org/wiki/Upwind_scheme).

That said, `MOLFiniteDifference` takes an `upwind_order` parameter, but the manual doesn’t say what its effect is.  
@xtalax It appears that upwind discretization should be performed automatically…? Can you briefly comment? I’m having a hard time understanding the source 😄

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [May 17, 2022, 2:43pm UTC](https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567/5 "2022-05-17T14:43:42Z")

</div>

Yeah the source relies heavily on Symbolics.jl and SymbolicUtils.jl, not so easy to understand if you don’t have experience with it. To see where terms are recognised and schemes applied see [this](https://github.com/SciML/MethodOfLines.jl/blob/e24d5b904a5c1049933e69f7eefa75226c9bccf4/src/discretization/generate_finite_difference_rules.jl#L158) file, also useful to understand this is the [documentation on `@rule`](https://symbolicutils.juliasymbolics.org/rewrite/). Short answer, yes advection terms are automatically intercepted and the correct upwind scheme applied :). Any more questions on functionality welcome!

`upwind_order` is the aproximation order of the upwind scheme, note that this appears to lead to instability when it is above 1, which is its default, which I suspect is due to offside stencils near the boundaries.

Current work is on implementing flux limiters to ensure stability and higher order accuracy for advection terms in a wider set of circumstances.

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [May 17, 2022, 2:51pm UTC](https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567/6 "2022-05-17T14:51:31Z")

</div>

Ah another thing @sdwfrost, you want to use `@register_symbolic` rather than `@register`, there was an error in the docs - Sorry!

---

<div class="post-metadata">

**Author:** ![skleinbo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skleinbo/32/36080_2.png) [@skleinbo](https://discourse.julialang.org/u/skleinbo)\
**Post date:** [May 18, 2022, 6:46am UTC](https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567/7 "2022-05-18T06:46:28Z")

</div>

Thank you for pointing me in the right direction!

> advection terms are automatically intercepted and the correct upwind scheme applied

Cool! That’s a neat feature, yet as far as I can tell not mentioned explicitly in the docs. Neither here [How it Works · MethodOfLines.jl](http://methodoflines.sciml.ai/dev/howitworks/) nor here [MOLFiniteDifference · MethodOfLines.jl](http://methodoflines.sciml.ai/dev/MOLFiniteDifference/).

In OPs example, the discretizer seems to pick the wrong direction though:

```julia
symbolic_discretize(pde_sys, discretization)
(ODESystem(Equation[10.0u[3](t) + Differential(t)(u[2](t)) - 10.0u[2](t) ~ -μ*u[2](t),
[...]

```

Explains the instability. Indeed, bringing the `Da` term to the rhs solves the issue, i.e. this integrates fine:

```julia
pde_eqns = [Dt(u(t,a)) ~ -μ*u(t,a) - Da(u(t,a))]
domains = [t ∈ Interval(0,tmax), a ∈ Interval(0,amax)]
function a0(a)
    if da < a <= 2da # don't start right at the boundary.
        return 1.0/da
    else
        return 0.0
    end
end
@register_symbolic a0(a)
bcs = [u(t,0) ~ 0.0, u(0,a) ~ a0(a)]
@named pde_sys = PDESystem(pde_eqns, bcs, domains, [t,a], [u(t,a)], p)
discretization = MOLFiniteDifference([a=>da], t, approx_order = 2, grid_align = center_align)
pde_prob = discretize(pde_sys, discretization)

pde_sol = solve(pde_prob, Tsit5(), saveat = dt)

```

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [May 18, 2022, 11:06am UTC](https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567/8 "2022-05-18T11:06:12Z")

</div>

> Cool! That’s a neat feature, yet as far as I can tell not mentioned explicitly in the docs. Neither here [How it Works · MethodOfLines.jl](http://methodoflines.sciml.ai/dev/howitworks/) nor here [MOLFiniteDifference · MethodOfLines.jl](http://methodoflines.sciml.ai/dev/MOLFiniteDifference/).

You’re right, I will add this.

> Explains the instability. Indeed, bringing the Da term to the rhs solves the issue, i.e. this integrates fine:

Good catch! That is indeed a big problem, pushing a hotfix as soon as possible.

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [May 18, 2022, 12:18pm UTC](https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567/9 "2022-05-18T12:18:29Z")

</div>

This appears to be related to advection terms appearing on the same side as the time derivative, this causes the ‘wrong’ scheme to be dispatched from the perspective of the case when that term is on the other side of the equation - both positive and negative coefficients fail in this case.

I will update the docs to reflect this limitation until such a time as we can automatically rearrange equations in to a set form pre discretization.

---

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [December 19, 2022, 6:53pm UTC](https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567/10 "2022-12-19T18:53:08Z")

</div>

Hi! I’m revisiting this after a while, as I need to understand more about how boundary conditions are specified. Here’s my current code - I’ve tried various ways of writing the BCs in terms of da and dt, but haven’t got a match yet. Any suggestions of where I’m going wrong? The advection term on the LHS seems to work OK now.

```julia
using ModelingToolkit
using MethodOfLines
using OrdinaryDiffEq
using DomainSets
using Plots

@parameters t a μ
@variables v(..) u(..)
Dt = Differential(t)
Da = Differential(a)
p = [μ => 0.04]
tmax = 100.0
amax = 100.0
da = 1.0
dt = 1.0

# ODE version
ode_eqns = [Dt(v(t)) ~ -μ*v(t)]
v0 = [v(t) => 1.0]
@named ode_sys = ODESystem(ode_eqns)
ode_prob = ODEProblem(ode_sys, v0, (0,tmax), p)
ode_sol = solve(ode_prob, Tsit5(), saveat = dt)

# PDE version
pde_eqns = [Dt(u(t,a)) + Da(u(t,a)) ~ -μ*u(t,a)]
domains = [t ∈ Interval(0,tmax), a ∈ Interval(0,amax)]

function get_bcs(t,a)
    if t < dt
        if a < da
            return 1.0
        else
            return 0.0
        end
    else
        return 0.0
    end
end
@register_symbolic get_bcs(t,a)

bcs = [u(t,0) ~ get_bcs(t,0.0), u(0,a) ~ get_bcs(0.0,a)]
@named pde_sys = PDESystem(pde_eqns, bcs, domains, [t,a], [u(t,a)], p)
discretization = MOLFiniteDifference([a=>da], t, approx_order = 2, grid_align = center_align)
pde_prob = discretize(pde_sys, discretization)
pde_sol = solve(pde_prob, Tsit5(), saveat = dt)

umat = pde_sol.u[u(t, a)]
usum = sum(umat, dims=1)

plot(ode_sol)
plot!(pde_sol.t, usum')

```

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 19, 2022, 7:30pm UTC](https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567/11 "2022-12-19T19:30:19Z")

</div>

I also had issues simulating those. In the end, mass conservation proved crucial and so an upwind finite volume scheme was employed. You can find some details in the Phd thesis of my student [Long time behavior of a mean-field model of interacting spiking neurons - Inria - Institut national de recherche en sciences et technologies du numérique](https://hal.inria.fr/tel-03123384)

---

<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:** [December 19, 2022, 7:44pm UTC](https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567/12 "2022-12-19T19:44:20Z")

</div>

Then maybe try the WENO scheme?

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [December 19, 2022, 8:07pm UTC](https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567/13 "2022-12-19T20:07:36Z")

</div>

So to understand, you want the corner point to be 1 and the rest of the boundary to be 0?

The way this is currently specified, this will have no effect on the solution compared with 0 bcs and ics since the corners are unused by the discretization, due to ghost cell padding.

Could you try replacing this with a gaussian pulse centered on the time space corner, perhaps one or two steps in? This will have very similar properties to the discrete delta, will reduce numerical dispersion, and have some chance of inducing the effect you are looking for in the dynamics.

EDIT: Actually, the time space corner is used, the trouble is that most of the effect of this delta will be wiped out by the zero bc at the second timestep. So centring your gaussian at t = 0 a = da should work

---

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [December 19, 2022, 8:12pm UTC](https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567/14 "2022-12-19T20:12:46Z")

</div>

@xtalax - yes, I want the corner point to be 1 and the boundary to be zero.

@ChrisRackauckas - I tried changing to a WENO scheme, but now I get `ERROR: LinearAlgebra.SingularException(98)` with the following:

```julia
using DifferentialEquations
using ModelingToolkit
using MethodOfLines
using OrdinaryDiffEq
using DomainSets
using Plots

@parameters t a μ
@variables v(..) u(..)
Dt = Differential(t)
Da = Differential(a)
p = [μ => 0.04]
tmax = 100.0
amax = 100.0
da = 1.0
dt = 1.0

# ODE version
ode_eqns = [Dt(v(t)) ~ -μ*v(t)]
v0 = [v(t) => 1.0]
@named ode_sys = ODESystem(ode_eqns)
ode_prob = ODEProblem(ode_sys, v0, (0,tmax), p)
ode_sol = solve(ode_prob, Tsit5(), saveat = dt)

# PDE version
pde_eqns = [Dt(u(t,a)) + Da(u(t,a)) ~ -μ*u(t,a)]
domains = [t ∈ Interval(0,tmax), a ∈ Interval(0,amax)]

function get_bcs(t,a)
    if t < dt
        if a < da
            return 1.0
        else
            return 0.0
        end
    else
        return 0.0
    end
end
@register_symbolic get_bcs(t,a)

bcs = [u(t,0) ~ get_bcs(t,0.0), u(0,a) ~ get_bcs(0.0,a)]
@named pde_sys = PDESystem(pde_eqns, bcs, domains, [t,a], [u(t,a)], p)
discretization = MOLFiniteDifference([a=>da], t, approx_order = 2, grid_align = center_align, advection_scheme=WENOScheme())
pde_prob = discretize(pde_sys, discretization)
pde_sol = solve(pde_prob, saveat = dt)

umat = pde_sol.u[u(t, a)]
usum = sum(umat, dims=1)

plot(ode_sol)
plot!(pde_sol.t, usum')

```

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [December 19, 2022, 8:20pm UTC](https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567/15 "2022-12-19T20:20:47Z")

</div>

See my above edit :).

Yes, I’ve seen this one before, arises in various circumstances with WENO that I haven’t yet been able to nail down. This kind of system error is notoriously difficult to debug.

From what I can tell at the moment, this shouldn’t be present if you can provide 2 bcs per spatial boundary that won’t affect your dynamics too much, which seems unlikely given that your ic starts at the boundary.

In any case please try the gaussian pulse, pure discrete dirac deltas are very rarely well handled by finite difference methods. You might also try a spectral method, though since MOL does not yet support these you wiill have to do this manually, for now.

---

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [December 19, 2022, 8:28pm UTC](https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567/16 "2022-12-19T20:28:49Z")

</div>

@xtalax - thanks! I’ll give it a go. The models I deal with are mostly solved by the method of characteristics (see [this SO post](https://math.stackexchange.com/questions/3724507/characteristic-method-for-mckendrick-von-foerster-age-structure-model) for a solution of the above) rather than finite differences, but I wanted to see how well MethodOfLines.jl could do in this case.