# 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