# DiffEqOperators and DifferentialEquations: solving a coupled system of PDEs
**URL:** https://discourse.julialang.org/t/diffeqoperators-and-differentialequations-solving-a-coupled-system-of-pdes/74722
**Category:** Numerics
**Tags:** differentialequation
**Created:** [January 17, 2022, 12:49am UTC](https://discourse.julialang.org/t/diffeqoperators-and-differentialequations-solving-a-coupled-system-of-pdes/74722 "2022-01-17T00:49:51Z")
**Posts on this page:** 15
**Page:** 1
### Author:  [@ash](https://discourse.julialang.org/u/ash)
#### Post date: [January 17, 2022, 12:49am UTC](https://discourse.julialang.org/t/diffeqoperators-and-differentialequations-solving-a-coupled-system-of-pdes/74722/1 "2022-01-17T00:49:51Z")
I am trying to solve a system of coupled PDEs shown below:
and
.
I have written some simple code using DifferentialEquations.jl and DiffEqOperators.jl hoping to solve this system:
```julia
using DifferentialEquations
using DiffEqOperators
# 300x300 grid, x \in 0,30 and t \in 0,0.3
# Define FD operator
Nₓ = 300
Δx = 30.0/(Nₓ)
x = range(0, step=Δx, length=Nₓ)
ord_deriv = 2
ord_approx = 2
Δ = CenteredDifference(ord_deriv, ord_approx, Δx, Nₓ)
bc = Dirichlet0BC(Float64)
#Parameters
f = 1.2
m = 190
q = 0.001
ϵ₁ = 0.01
ϵ = 0.01
K = 1000
Θ = 0.0025
d₀ = 0.01
D₁ = 1
p = [f, m, q, ϵ, ϵ₁, D₁, d₀, K, Θ]
function myModel!(du,u,p,t)
f, m, q, ϵ, ϵ₁, D₁, d₀, K, Θ = p
c₁, c₂ = u
du[1] = 1/ϵ*(f*c₂.*(q .- c₁)./(q .+ c₁) .+ c₁ .* (1 .- m*c₂)./(1 .- m*c₂ .+ ϵ₁) .- c₁.^2) .+
D₁*Δ*bc*c₁
du[2] = c₁ .* (1 .- m*c₂)./(1 .- m*c₂ .+ ϵ₁) .- c₂ .+ D₁/d₀*(1 + (K-1)/(1+3/Θ))^(-3/2).*Δ*bc*c₂
end
# Prepare the ODE intnegrator
tspan = (0.0,0.3)
u0 = rand(Float64, (1, Nₓ))
prob = ODEProblem(myModel!, u0, tspan, p)
dt = 0.001
solve(prob, Tsit5(), dt = dt)
```
However, I am greeted with this error:
```julia
ERROR: LoadError: MethodError: no method matching iterate(::DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Vector{Float64}, Int64})
```
which makes me think I’m doing something wrong with DiffEqOperators, but it seems to work when you just have one equation. Any thoughts?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 17, 2022, 12:47pm UTC](https://discourse.julialang.org/t/diffeqoperators-and-differentialequations-solving-a-coupled-system-of-pdes/74722/2 "2022-01-17T12:47:23Z")
That looks like a broadcasting issue on your side. It shouldn’t be itereating the DiffEqOperator but it is, so somewhere you have a `.*` where it should be `*`.
---
### Author:  [@ash](https://discourse.julialang.org/u/ash)
#### Post date: [January 17, 2022, 4:29pm UTC](https://discourse.julialang.org/t/diffeqoperators-and-differentialequations-solving-a-coupled-system-of-pdes/74722/3 "2022-01-17T16:29:51Z")
Thanks! Managed to hunt this down and fix it. To further debug the code I’ve simplified it as much as possible to the following:
```julia
using DiffEqOperators, OrdinaryDiffEq
# Define FD operator
Nₓ = 300
Δx = 30.0/(Nₓ)
x = range(0, step=Δx, length=Nₓ)
ord_deriv = 2
ord_approx = 2
const Δ = CenteredDifference(ord_deriv, ord_approx, Δx, Nₓ)
const bc = Dirichlet0BC(Float64)
function myModel!(du,u,p,t)
du[1] = Δ*bc*u[1]
#du[2] .= D₁/d₀*(1 + (K-1)/(1+3/Θ))^(-3/2)*Δ₂*bc*u[2]
du[2] = Δ*bc*u[2]
end
# Prepare the ODE intnegrator
tspan = (0.0,0.3)
u0 = rand(Float64, (1, Nₓ))
prob = ODEProblem(myModel!, u0, tspan, p)
dt = 0.001
solve(prob, Tsit5(), dt = dt)
```
However I am getting the error:
```julia
ERROR: LoadError: MethodError: no method matching *(::GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Vector{Float64}, Int64}, RobinBC{Float64, Vector{Float64}}}, ::Float64)
```
It seems to work perfectly fine in the single equation case though. I’ve already tried everything I can think of to debug this but will keep working on it.
---
### Author:  [@CollinsPeriodLee](https://discourse.julialang.org/u/CollinsPeriodLee)
#### Post date: [January 17, 2022, 6:14pm UTC](https://discourse.julialang.org/t/diffeqoperators-and-differentialequations-solving-a-coupled-system-of-pdes/74722/4 "2022-01-17T18:14:09Z")
You could try this, it works
```julia
using DiffEqOperators, OrdinaryDiffEq
# Define FD operator
Nₓ = 300
Δx = 30.0/(Nₓ)
x = range(0, step=Δx, length=Nₓ)
ord_deriv = 2
ord_approx = 2
const Δ = CenteredDifference(ord_deriv, ord_approx, Δx, Nₓ)
const bc = Dirichlet0BC(Float64)
function myModel!(du, u,p,t)
du[1,:] .= Δ* bc * u[1,:]
#du[2] .= D₁/d₀*(1 + (K-1)/(1+3/Θ))^(-3/2)*Δ₂*bc*u[2]
du[2,:] .= Δ* bc * u[2,:]
end
# Prepare the ODE intnegrator
tspan = (0.0,0.3)
u0 = rand(Float64, (2, Nₓ))
prob = ODEProblem(myModel!, u0, tspan)
dt = 0.001
solve(prob, Tsit5(), dt = dt)
```
Since you have 2 PDEs, the initial condition should be `u0 = rand(Float64, (2, Nₓ))` and you also need to match `du` & `u` in the function `mymodel!`
```julia
function myModel!(du, u,p,t)
du[1,:] .= Δ* bc * u[1,:]
#du[2] .= D₁/d₀*(1 + (K-1)/(1+3/Θ))^(-3/2)*Δ₂*bc*u[2]
du[2,:] .= Δ* bc * u[2,:]
end
```
---
### Author:  [@ash](https://discourse.julialang.org/u/ash)
#### Post date: [January 18, 2022, 2:47am UTC](https://discourse.julialang.org/t/diffeqoperators-and-differentialequations-solving-a-coupled-system-of-pdes/74722/5 "2022-01-18T02:47:43Z")
That works, thank you!
The full code now is:
```julia
using OrdinaryDiffEq
using DiffEqOperators, LinearAlgebra
using Plots
#Parameters
f = 1.2
m = 190
q = 0.001
ϵ₁ = 0.01
ϵ = 0.01
K = 1000
Θ = 0.0001
d₀ = 0.01
D₁ = 1
p = [f, m, q, ϵ, ϵ₁, D₁, d₀, K, Θ]
# Define FD operator
Nₓ = 300
Δx = 15.0/(Nₓ)
x = range(0, step=Δx, length=Nₓ)
ord_deriv = 2
ord_approx = 2
Δ = CenteredDifference(ord_deriv, ord_approx, Δx, Nₓ)
bc = Neumann0BC(Δx, 1)
function myModel!(du,u,p,t)
f, m, q, ϵ, ϵ₁, D₁, d₀, K, Θ = p
c₁ = u[1,:]
c₂ = u[2,:]
du[1,:] .= 1/ϵ*(f*c₂.*(q .- c₁)./(q .+ c₁) .+ c₁ .* (1 .- m*c₂)./(1 .- m*c₂ .+ ϵ₁) .- c₁.^2) .+ D₁*Δ*bc*c₁
du[2,:] .= c₁ .* (1 .- m*c₂)./(1 .- m*c₂ .+ ϵ₁) .- c₂ .+ D₁/d₀*(1 + (K-1)/(1+3/Θ))^(-3/2)*Δ*bc*c₂
end
# Grid
# 300x300
# Spatial step: 0.1 s.u.
# Temporal step: 0.001 t.u.
# i.e.,
# Time runs from 0 to 0.3 and space from 0 to 30
# Prepare the ODE intnegrator
tspan = (0.0,0.3)
u0 = rand(Float64, (2, Nₓ))
prob = ODEProblem(myModel!, u0, tspan, p)
dt = 0.001
sol = solve(prob, Tsit5(), dt = dt, reltol=1e-8, abstol=1e-8) # Should probably use KenCarp4 but it makes no difference
# Plot
c1 = []
c2 = []
for i in 1:length(sol.t)
append!(c1, sol.u[i][1,:])
append!(c2, sol.u[i][2,:])
end
c1r = reshape(c1, (Nₓ, length(sol.t)))
c2r = reshape(c2, (Nₓ,length(sol.t)))
p1 = heatmap(x, sol.t, c1r')
p2 = heatmap(x, sol.t, c2r')
plot(p1, p2, layout = (2,1), fmt=:png)
```
which runs perfectly fine. Now I just have to figure out why it produces no patterns
Expected solution (for c\_1, presumably, they don’t say):
(depending on the value of \Theta, see [here](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056205#fulltext)).
I get:
I tried playing around with dt and dx, approximation orders, etc. I have double-checked the model from two papers. I guess this will be an adventure. 😅
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 18, 2022, 4:35am UTC](https://discourse.julialang.org/t/diffeqoperators-and-differentialequations-solving-a-coupled-system-of-pdes/74722/6 "2022-01-18T04:35:01Z")
> [@ash](#):
>
> ```julia-auto
> sol = solve(prob, Tsit5(), dt = dt, reltol=1e-8, abstol=1e-8) # Should probably use KenCarp4 but it makes no difference
>
> ```
You are still using adaptive time stepping here since you did not set `adaptive=false` and instead only set the initial dt.
---
### Author:  [@ash](https://discourse.julialang.org/u/ash)
#### Post date: [January 18, 2022, 4:47am UTC](https://discourse.julialang.org/t/diffeqoperators-and-differentialequations-solving-a-coupled-system-of-pdes/74722/7 "2022-01-18T04:47:03Z")
Thank you. I had already tried that as well, it requires a much finer `dt` (~0.00001) to avoid instability and gives the same solution. I am not quite sure what is going on.
---
### Author:  [@ash](https://discourse.julialang.org/u/ash)
#### Post date: [January 18, 2022, 5:04am UTC](https://discourse.julialang.org/t/diffeqoperators-and-differentialequations-solving-a-coupled-system-of-pdes/74722/8 "2022-01-18T05:04:06Z")
I also tried enforcing the [no flux](https://sites.me.ucsb.edu/~moehlis/APC591/tutorials/tutorial5/node6.html) (Neumann 0) boundary condition on the initial condition but the result didn’t change much.
```julia
u0[:,2] .= u0[:,1]
u0[:, Nₓ-1] .= u0[:, Nₓ]
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 21, 2022, 9:35pm UTC](https://discourse.julialang.org/t/diffeqoperators-and-differentialequations-solving-a-coupled-system-of-pdes/74722/9 "2022-01-21T21:35:14Z")
Are you sure they are discretizing the same way? Did you look through their code?
---
### Author:  [@ash](https://discourse.julialang.org/u/ash)
#### Post date: [January 21, 2022, 9:41pm UTC](https://discourse.julialang.org/t/diffeqoperators-and-differentialequations-solving-a-coupled-system-of-pdes/74722/10 "2022-01-21T21:41:20Z")
They were discretizing quite differently but the issue here specifically is that I was using (1+1)D and they were using (2+1)D, so a very silly mistake.
Even after using (2+1)D however I haven’t been able to fully reproduce the results of this model. I am not sure what the issue is (could be discretization, they use a fully explicit scheme?). With a little effort, I moved my code over to the Gray-Scott model and it works flawlessly; the results perfectly agree with the literature, so I think I’m done with this specific model for now.
Thanks again for the help!
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 22, 2022, 2:04am UTC](https://discourse.julialang.org/t/diffeqoperators-and-differentialequations-solving-a-coupled-system-of-pdes/74722/11 "2022-01-22T02:04:32Z")
> [@ash](#):
>
> Even after using (2+1)D however I haven’t been able to fully reproduce the results of this model. I am not sure what the issue is (could be discretization, they use a fully explicit scheme?)
How accurate is their method? Euler with a big `dt`? If someone uses a very inaccurate method you can still get qualitatively okay results with quantitatively wrong boundaries for the different behaviors. Then no matter how well you try to match you might not, because being more accurate may mean matching their results less 😅. I’ve seen that in bio PDE results before with some silly "run once with a big Euler and no error estimation.
---
### Author:  [@ash](https://discourse.julialang.org/u/ash)
#### Post date: [January 23, 2022, 6:06am UTC](https://discourse.julialang.org/t/diffeqoperators-and-differentialequations-solving-a-coupled-system-of-pdes/74722/12 "2022-01-23T06:06:44Z")
They use a fully explicit three-level Du Fort-Frankel scheme with dt = 0.001, which is too large for the explicit integrators I tried in DifferentialEquations.jl (with adaptivity off of course), so I think something might be going on there. The other unknown is the exact nature of the initial condition, they just say “random noise” without giving more details.
Thanks for the advice, this is important to know!
---
### Author:  [@ziolai](https://discourse.julialang.org/u/ziolai)
#### Post date: [February 2, 2022, 6:22pm UTC](https://discourse.julialang.org/t/diffeqoperators-and-differentialequations-solving-a-coupled-system-of-pdes/74722/13 "2022-02-02T18:22:01Z")
@ash: thank you so much for sharing a wonderful example! Could you please share a reference for the D\_1 / D\_2 expression? The Guiu-Souto-2010 paper does not seem to contain it. Thx!
[edit]: I found the expression for D\_1 / D\_2 is a more recent (2012) paper by the same authors. Apologies for the noice.
---
### Author:  [@ziolai](https://discourse.julialang.org/u/ziolai)
#### Post date: [February 3, 2022, 9:43am UTC](https://discourse.julialang.org/t/diffeqoperators-and-differentialequations-solving-a-coupled-system-of-pdes/74722/14 "2022-02-03T09:43:24Z")
@ash: I am developing a notebook for my student on Turing patterns using your implementation as point of departure. I would like to acknowledge your contribution. Your contact details would be good to have. Thank you.
---
### Author:  [@ash](https://discourse.julialang.org/u/ash)
#### Post date: [February 3, 2022, 2:14pm UTC](https://discourse.julialang.org/t/diffeqoperators-and-differentialequations-solving-a-coupled-system-of-pdes/74722/15 "2022-02-03T14:14:55Z")
@ziolai My (slightly outdated) website is in my discourse bio, it has all my contact information, thanks!
Several things to note about this model
1. It’s (2+1) D, i was treating it as (1+1)D in this thread which is wrong. I was just confused by the paper.
2. I haven’t been able to reproduce the patterns in the paper exactly. They use an _explicit_ method with a large dt and I was using a far more accurate method but still couldn’t reproduce the thresholds for the phase transitions between dots/labyrinth etc.
3. The Gray-Scott model is much more fun for learning about Turing patterns in my opinion. See [here](http://mrob.com/pub/comp/xmorphia/index.html). I have implemented it in [this](https://github.com/oashour/PatternFormation.jl) unlisted and very preliminary package, but it is very heavily optimized (see my post history). There is an example in the root of the repository. Use the :NoisePatches `init_cond` to reproduce mrob’s results.
The code is a little messy and is not documented at the moment, so let me know if you have any questions!