# Intro

I’m using a semi-coupled system of wave and diffusion equations to simulate a string whose mass/speed of sound diffuses (as a toy to familiarize myself with `MethodOfLines.jl` for my master’s thesis in high-performance computational physics). The initial profiles are both normal distributions, and as such I choose the spatial domain to be of length 4*pi. But I encounter an error in the discretization.

Given a step size `dx = pi / n`, a warning is issued

Warning: dx for x does not divide domain exactly, adding grid point at x = 6.283185307179586).

followed by a `BoundsError` such as

ERROR: LoadError: BoundsError: attempt to access 401-element Vector{SymbolicUtils.BasicSymbolic{Real}} at index
[CartesianIndex{1}[CartesianIndex(0,), CartesianIndex(1,)]]

any time `n` is a multiple of 25, as far as I have tested. Additionally, no flags are raised when `n = 5`. I tested for multiple different dependencies on `n`, e.g. multiples of 5 not equal to multiples of 25, ±1 around multiples of 25, and all have worked just fine.

I’ve tried using the VSCode debugger, but in this particular instance, it quickly gets over my head.

# Code (does not work)

It appears the `BoundsError` is raised from within `discretize`. As such I have omitted the later code to actually get and visualize the solution.

Note: `dx = pi / 25.`

``````using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets

# independent and dependent variables
@parameters t x        # t - time, x - space
@variables u(..) c(..) # u - wave amplitude, c - speed of sound

# differential operators
Dt  = Differential(t)   # first  order time derivative operator
Dtt = Differential(t)^2 # second order time derivative operator
Dx  = Differential(x)   # first  order spatial derivative operator
Dxx = Differential(x)^2 # second order spatial derivative operator

# governing equations
eqwave  = Dtt(u(t, x)) ~ c(t, x)^2 * Dxx(u(t, x)) # wave equation with variable speed
eqspeed =  Dt(c(t, x)) ~ Dxx(c(t, x))             # diffusion equation for mass/speed

# boundary conditions
bcs = [
# initial conditions
u(0, x)  ~ exp(-x^2),   # initial spatial wave profile
Dt(u(0, x)) ~ 0.,          # initial wave velocity
c(0, x)  ~ exp(-x^2),   # initial spatial speed profile
# boundary conditions (optional)
Dx(u(t, -2pi)) ~ 0., # left side is rigid
u(t,  2pi)  ~ 0.  # right side held at zero
]

# Space and time domains
domains = [
t ∈ Interval(0.0, 20.0),      # arbitrary time domain
x ∈ Interval(-2 * pi, 2 * pi) # arbitrary spatial domain
]

# PDE system
@named pdesys = PDESystem([equ, eqc], bcs, domains, [t, x], [u(t, x), c(t, x)])

# Method of lines discretization
dx = pi / 25. # doesn't work when denominator is a multiple of 25
discretization = MOLFiniteDifference([x => dx], t)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys, discretization)
``````

# Versions

## Packages

`DomainSets` has been held back from updating due to compatibility constraints.

[0c46a032] DifferentialEquations v7.11.0
⌅ [5b8099bc] DomainSets v0.6.7
[7a1cc6ca] FFTW v1.7.1
[94925ecb] MethodOfLines v0.10.1
[961ee093] ModelingToolkit v8.72.2
[1dea7af3] OrdinaryDiffEq v6.59.0

## Julia

1.9.3+0.x64.w64.mingw32

ahh this is likely a floating point error thing. Can you open an issue?

Of course! I’ll do that now.

Yep, floating point, this often happens where pi is involved. I haven’t seen the bounds error though, I will investigate