# Concrete Operators in DiffEqOperators.jl for time-indepedent problems

**URL:** <https://discourse.julialang.org/t/concrete-operators-in-diffeqoperators-jl-for-time-indepedent-problems/77415>\
**Category:** Numerics\
**Tags:** question\
**Created:** [March 4, 2022, 5:29pm UTC](https://discourse.julialang.org/t/concrete-operators-in-diffeqoperators-jl-for-time-indepedent-problems/77415 "2022-03-04T17:29:02Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![ziolai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ziolai/32/23422_2.png) [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Post date:** [March 4, 2022, 5:29pm UTC](https://discourse.julialang.org/t/concrete-operators-in-diffeqoperators-jl-for-time-indepedent-problems/77415/1 "2022-03-04T17:29:02Z")

</div>

We are looking into the use of DiffEqOperators for time-independent problems. We are considering the Poisson equation on the unit interval and unit square.

**First** a 1D implementation was made to figure out how to use it.

DiffEqOperators provides several functions to construct finite difference approximations of a differential equation. The ‘matrix’ Ah is built up in two steps:

1. Construct the finite difference approximation using CenteredDifference or UpwindDifference
2. Construct the boundary conditions using one of several available functions for different types (e.g. DirichletBC for the Dirichlet boundary condition).  
The complete Ah ‘matrix’ is then constructed by multiplying the result of these two steps. For example,  
A = CenteredDifference(ord, ord\_approx, h, N)  
bc = Dirichlet0BC(Float64)  
Ah = A \* bc  
The word ‘matrix’ is between quotation marks because DiffEqOperators does not really return a matrix but an abstract operator. This operator can then be concretized when required. DiffEqOperators also provides an implementation of \ for solving the differential equation constructed using the operators.

u = (A \* bc) \ f;  
u = vcat(0, u, 0)

**Next** the same was tried in 2D

Dxx = CenteredDifference{1}(ord, ord\_approx, h, N); Dyy = CenteredDifference{2}(ord, ord\_approx, h, N); bc = Dirichlet0BC(Float64)  
Qx,Qy = MultiDimBC(bc, size(F))  
Axx = Dxx \* Qx; Ayy = Dyy \* Qy; A =Axx+Ayy;  
u = A \ f;

This, however, resulted in the error message

```
     MethodError: Cannot `convert` an object of type
       GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, Stati
     cArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Vector{Float64}, Int64}, MultiDimDirectionalBC{Float64,
     RobinBC{Float64, Vector{Float64}}, 1, 2, 1}} to an object of type
       AbstractMatrix

```

Please advice. Thx!

---

<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:** [March 6, 2022, 2:51pm UTC](https://discourse.julialang.org/t/concrete-operators-in-diffeqoperators-jl-for-time-indepedent-problems/77415/2 "2022-03-06T14:51:51Z")

</div>

Your post has just revealed a long standing issue with concretization in multiple dimensions that we will fix asap so thank you for that. In the meantime, can you share the system you are trying to solve and I’ll find a way to solve it with the current infrastructure?

---

<div class="post-metadata">

**Author:** ![ziolai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ziolai/32/23422_2.png) [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Post date:** [March 6, 2022, 4:40pm UTC](https://discourse.julialang.org/t/concrete-operators-in-diffeqoperators-jl-for-time-indepedent-problems/77415/3 "2022-03-06T16:40:46Z")

</div>

Thx! We are trying to solve Laplacian(u) = f on the domain 0 \< x,y \< 1. Nothing fancy. We have a home-brewed implementation that we can share. Please let us know in case we can contribute. We will be happy to test drive your solution. Thx again.

---

<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:** [March 7, 2022, 1:56pm UTC](https://discourse.julialang.org/t/concrete-operators-in-diffeqoperators-jl-for-time-indepedent-problems/77415/4 "2022-03-07T13:56:59Z")

</div>

Hi, in the meantime I suggest using our other package [`MethodOfLines`](https://github.com/SciML/MethodOfLines.jl), which uses Symbolics.jl and DiffEqOperators.jl under the hood and existing ODE/Nonlinear solvers. Here’s an example, please ask away if you have any questions. Note that this does not currently scale well to very large problems, but there are improvements in the works that will hopefully make this more reasonable.

```julia
using ModelingToolkit, MethodOfLines, DomainSets, NonlinearSolve

# Define the system of equations
@parameters x y
@variables u(..)

Dxx = Differential(x)^2
Dyy = Differential(y)^2

xmin = 0.
xmax = 1.
ymin = 0.
ymax = 1.

f(x,y) = x^2 + y^2 #If you have array data, define a piecewise function or use an interpolation

eq = Dxx(u(x,y)) + Dyy(u(x,y)) ~ f(x,y)

domain = [x ∈ Interval(xmin,xmax), 
          y ∈ Interval(ymin,ymax)]

# Arbitrary bcs can be defined in a natural way
# Dirichlet 0 would be u(xmin,y) ~ 0.
# Neumann would be Differential(x)(u(xmax,y)) ~ α(x, y) where α is a function, can be replaced with a constant value.
# Periodic would be u(xmin,y) ~ u(xmax,y)
bcs = [u(xmin, y) ~ u(xmax, y),
       u(x, ymin) ~ u(x, ymax)]

@named pdesys = PDESystem(eq, bcs, domain, [x, y], [u(x,y)])

# Create discretization

dx = 0.1
dy = 0.1

order = 2 # Order of the finite element approximation, must be even

discretization = MOLFiniteDifference([x => dx, y => dy], approx_order = order, grid_align = center_align) # edge_align gives better results with Neumann BCs, though uses extra grid points

prob = discretize(pdesys, discretization)

sol = NonlinearSolve.solve(prob, NewtonRaphson())

# Center aligned grid
disc_x = xmin:dx:xmax
disc_y = ymin:dy:ymax

# Edge aligned grid
# disc_x_e = xmin-dx/2:dx:xmax+dx/2
# disc_y_e = ymin-dy/2:dy:ymax+dy/2

sol_u = reshape(sol.u, (length(disc_x), length(disc_y)))

# Plot the solution
using Plots
heatmap(disc_x, disc_y, sol_u)

```

---

<div class="post-metadata">

**Author:** ![ziolai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ziolai/32/23422_2.png) [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Post date:** [March 7, 2022, 2:14pm UTC](https://discourse.julialang.org/t/concrete-operators-in-diffeqoperators-jl-for-time-indepedent-problems/77415/5 "2022-03-07T14:14:49Z")

</div>

Thx. I currently experience a difficulty installing MethodofLines. I will try to resolve it.

---

<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:** [March 7, 2022, 2:24pm UTC](https://discourse.julialang.org/t/concrete-operators-in-diffeqoperators-jl-for-time-indepedent-problems/77415/6 "2022-03-07T14:24:50Z")

</div>

What’s the issue? If it comes up for you chances are it will come up for others.

---

<div class="post-metadata">

**Author:** ![ziolai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ziolai/32/23422_2.png) [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Post date:** [March 7, 2022, 2:36pm UTC](https://discourse.julialang.org/t/concrete-operators-in-diffeqoperators-jl-for-time-indepedent-problems/77415/7 "2022-03-07T14:36:15Z")

</div>

Not sure. I fail to understand the output of Pkg (see below).

julia\> import Pkg; Pkg.add(“MethodOfLines”)  
Resolving package versions…  
ERROR: Unsatisfiable requirements detected for package Symbolics [0c5d862f]:  
Symbolics [0c5d862f] log:  
├─possible versions are: 0.1.0-4.3.0 or uninstalled  
├─restricted to versions \* by an explicit requirement, leaving only versions 0.1.0-4.3.0  
├─restricted by compatibility requirements with MethodOfLines [94925ecb] to versions: 4.0.0-4.3.0  
│ └─MethodOfLines [94925ecb] log:  
│ ├─possible versions are: 0.1.0 or uninstalled  
│ └─restricted to versions \* by an explicit requirement, leaving only versions 0.1.0  
└─restricted by compatibility requirements with StaticArrays [90137ffa] to versions: 0.1.0-0.1.32 or uninstalled — no versions left  
└─StaticArrays [90137ffa] log:  
├─possible versions are: 0.8.0-1.4.1 or uninstalled  
├─restricted to versions \* by an explicit requirement, leaving only versions 0.8.0-1.4.1  
└─restricted by compatibility requirements with ReactionMechanismSimulator [c2d78dd2] to versions: 0.8.0-0.12.5  
└─ReactionMechanismSimulator [c2d78dd2] log:  
├─possible versions are: 0.1.0-0.4.0 or uninstalled  
└─restricted to versions \* by an explicit requirement, leaving only versions 0.1.0-0.4.0

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [March 7, 2022, 2:45pm UTC](https://discourse.julialang.org/t/concrete-operators-in-diffeqoperators-jl-for-time-indepedent-problems/77415/8 "2022-03-07T14:45:47Z")

</div>

It is telling you that `StaticArrays` is upper bounded to version 0.12 by `ReactionMechanismSimulator`, but such an ancient version of `StaticArrays` means nly very old verions of `Symbolics` can be installed. `MethodOfLines` requires `Symbolics` version 4.0 or up, while the last version of `Symbolics` to work with `StaticArrays` version 0.12 is 0.1.32.

In general `ReactionMechanismSimulator` appears to depend on a lot of old package versions, which will make it tricky to use this package together with other packages. I do see some activity on their repo though, so maybe they’re in the process of updating their dependencies.

---

<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:** [March 7, 2022, 3:15pm UTC](https://discourse.julialang.org/t/concrete-operators-in-diffeqoperators-jl-for-time-indepedent-problems/77415/9 "2022-03-07T15:15:04Z")

</div>

@ziolai I noticed that you were interested in the [Brusselator equation](https://discourse.julialang.org/t/reference-for-brusselator-in-differentialequation-jl/75810/3), you should be aware that I’m currently working on allowing this equation in `MethodOfLines`.

What is your application, if I can ask?

---

<div class="post-metadata">

**Author:** ![ziolai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ziolai/32/23422_2.png) [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Post date:** [March 7, 2022, 3:31pm UTC](https://discourse.julialang.org/t/concrete-operators-in-diffeqoperators-jl-for-time-indepedent-problems/77415/10 "2022-03-07T15:31:59Z")

</div>

@nilshg : I removed ReactionMechanismSimulator and now have MethodsOfLines up and running.

@xtalax : my objective it to introduce solving reaction-diffusion equations in a master course that I am teaching.

Thx both!