# Solving linear system with summed DiffEqOperators

**URL:** <https://discourse.julialang.org/t/solving-linear-system-with-summed-diffeqoperators/60998>\
**Category:** General Usage\
**Tags:** diffeq\
**Created:** [May 12, 2021, 4:39am UTC](https://discourse.julialang.org/t/solving-linear-system-with-summed-diffeqoperators/60998 "2021-05-12T04:39:27Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![bclyons12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bclyons12/32/33734_2.png) [@bclyons12](https://discourse.julialang.org/u/bclyons12)\
**Post date:** [May 12, 2021, 4:39am UTC](https://discourse.julialang.org/t/solving-linear-system-with-summed-diffeqoperators/60998/1 "2021-05-12T04:39:27Z")

</div>

I’m trying to solve an elliptic differential equation that has multiple derivative terms and wanted to use `DiffEqOperators.jl` to build the matrices. Is it possible to sum `DiffEqOperator`s and then solve using backslash? When I try this, I get an error. Here’s an example, based on [https://discourse.julialang.org/t/how-to-discretize-and-solve-a-3d-poisson-equation-with-dirichlet-boundary-conditions-owing-to-diffeqoperators/40796](https://discourse.julialang.org/t/how-to-discretize-and-solve-a-3d-poisson-equation-with-dirichlet-boundary-conditions-owing-to-diffeqoperators/40796):

```julia
using DiffEqOperators
Δx = 0.1
x = Δx:Δx:1-Δx # Solve only for the interior: the endpoints are known to be zero!
N = length(x)
f = sin.(2π*x)

Δ = CenteredDifference(2, 2, Δx, N)
bc = Dirichlet0BC(Float64)

# Works
u = (Δ*bc) \ f
println("This solve is successful")

#Fails
u = ((0.5*Δ+0.5*Δ)*bc) \ f
println("But this solve fails")

```

with output

```julia
This solve is successful
ERROR: LoadError: MethodError: Cannot `convert` an object of type 
  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}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}} to an object of type 
  AbstractMatrix{T} where T
Closest candidates are:
  convert(::Type{AbstractMatrix{T} where T}, ::DerivativeOperator) at /Users/Lyons/.julia/packages/DiffEqOperators/DMNmH/src/derivative_operators/concretization.jl:63
  convert(::Type{AbstractMatrix{T} where T}, ::SciMLBase.DiffEqScaledOperator) at /Users/Lyons/.julia/packages/SciMLBase/XuLdB/src/operators/diffeq_operator.jl:97
  convert(::Type{AbstractMatrix{T} where T}, ::SciMLBase.DiffEqArrayOperator) at /Users/Lyons/.julia/packages/SciMLBase/XuLdB/src/operators/basic_operators.jl:102
  ...
Stacktrace:
  [1] (::DiffEqOperators.var"#117#118")(op::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}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}})
    @ DiffEqOperators ~/.julia/packages/DiffEqOperators/DMNmH/src/composite_operators.jl:38
  [2] MappingRF
    @ ./reduce.jl:93 [inlined]
  [3] afoldl(::Base.MappingRF{DiffEqOperators.var"#117#118", Base.BottomRF{typeof(Base.add_sum)}}, ::Base._InitialValue, ::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}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}, ::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}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}})
    @ Base ./operators.jl:533
  [4] _foldl_impl(op::Base.MappingRF{DiffEqOperators.var"#117#118", Base.BottomRF{typeof(Base.add_sum)}}, init::Base._InitialValue, itr::Tuple{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}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}, 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}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}})
    @ Base ./tuple.jl:263
  [5] foldl_impl(op::Base.MappingRF{DiffEqOperators.var"#117#118", Base.BottomRF{typeof(Base.add_sum)}}, nt::Base._InitialValue, itr::Tuple{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}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}, 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}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}})
    @ Base ./reduce.jl:48
  [6] mapfoldl_impl(f::DiffEqOperators.var"#117#118", op::typeof(Base.add_sum), nt::Base._InitialValue, itr::Tuple{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}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}, 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}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}})
    @ Base ./reduce.jl:44
  [7] mapfoldl(f::Function, op::Function, itr::Tuple{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}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}, 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}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}}; init::Base._InitialValue)
    @ Base ./reduce.jl:160
  [8] mapfoldl
    @ ./reduce.jl:160 [inlined]
  [9] #mapreduce#218
    @ ./reduce.jl:287 [inlined]
 [10] mapreduce
    @ ./reduce.jl:287 [inlined]
 [11] #sum#221
    @ ./reduce.jl:501 [inlined]
 [12] sum(f::Function, a::Tuple{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}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}, 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}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}})
    @ Base ./reduce.jl:501
 [13] convert(#unused#::Type{AbstractMatrix{T} where T}, L::DiffEqOperators.DiffEqOperatorCombination{Float64, Tuple{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}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}, 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}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}}, Vector{Float64}})
    @ DiffEqOperators ~/.julia/packages/DiffEqOperators/DMNmH/src/composite_operators.jl:37
 [14] \(L::DiffEqOperators.DiffEqOperatorCombination{Float64, Tuple{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}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}, 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}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}}, Vector{Float64}}, x::Vector{Float64})
    @ SciMLBase ~/.julia/packages/SciMLBase/XuLdB/src/operators/common_defaults.jl:14

```

---

<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:** [May 13, 2021, 1:33am UTC](https://discourse.julialang.org/t/solving-linear-system-with-summed-diffeqoperators/60998/2 "2021-05-13T01:33:57Z")

</div>

`u = SparseMatrixCSC((0.5*Δ+0.5*Δ)*bc) \ f` Hmm that’s missing a concretization. Can you open an issue?

---

<div class="post-metadata">

**Author:** ![briochemc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/briochemc/32/4209_2.png) [@briochemc](https://discourse.julialang.org/u/briochemc)\
**Post date:** [May 13, 2021, 1:34am UTC](https://discourse.julialang.org/t/solving-linear-system-with-summed-diffeqoperators/60998/3 "2021-05-13T01:34:59Z")

</div>

I’m not sure what the solution is, but FWIW, this happens because the ghost operators have their own backslash defined [here](https://github.com/SciML/DiffEqOperators.jl/blob/6fabb464146f2760abc284f8769ecb15fdb3cd98/src/derivative_operators/ghost_derivative_operator.jl#L27-L40), which I think is necessary since the ghost operator is not strictly a linear operator but an affine one AFAIU.

Maybe you need to add a `*(operator combination, bc)`, similar to [this](https://github.com/SciML/DiffEqOperators.jl/blob/6fabb464146f2760abc284f8769ecb15fdb3cd98/src/derivative_operators/ghost_derivative_operator.jl#L10-L12) for operator composition?

---

<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:** [May 13, 2021, 1:36am UTC](https://discourse.julialang.org/t/solving-linear-system-with-summed-diffeqoperators/60998/4 "2021-05-13T01:36:29Z")

</div>

Yeah exactly. It seems like it’s just a missing method.

---

<div class="post-metadata">

**Author:** ![bclyons12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bclyons12/32/33734_2.png) [@bclyons12](https://discourse.julialang.org/u/bclyons12)\
**Post date:** [May 13, 2021, 1:49am UTC](https://discourse.julialang.org/t/solving-linear-system-with-summed-diffeqoperators/60998/5 "2021-05-13T01:49:47Z")

</div>

Submitted as  
[https://github.com/SciML/DiffEqOperators.jl/issues/391](https://github.com/SciML/DiffEqOperators.jl/issues/391)

---

<div class="post-metadata">

**Author:** ![bclyons12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bclyons12/32/33734_2.png) [@bclyons12](https://discourse.julialang.org/u/bclyons12)\
**Post date:** [May 13, 2021, 4:32am UTC](https://discourse.julialang.org/t/solving-linear-system-with-summed-diffeqoperators/60998/6 "2021-05-13T04:32:58Z")

</div>

This does work:

```julia
A1, b1 = sparse(0.5*Δ*bc)
A2, b2 = sparse(0.5*Δ*bc)
As = A1 + A2
u = As \ f

```

so I can combine linear operators this way, though it’s less elegant than adding a method for `DiffEqOperatorCombination` as @briochemc suggested.

Separately but related, calling `sparse` on a `DiffEqOperator` is much slower and more memory intensive than building the same `SparseMatrixCSC` from three vectors. Here’s some example text, but I can open a separate issue/thread for this if it’s too off-topic:

```julia
N = 10001
Rᵢ = 20.
r = collect(range(0.,Rᵢ,length=N))
Δr = r[2]

function by_hand()

    R = zeros(3N-2)
    C = zeros(3N-2)
    V = zeros(3N-2)
    for i in 2:(N-1)

        # diagonal - 1
        j = 3*(i-1)
        R[j], C[j], V[j] = i, i-1, 1.0/Δr^2

        # diagonal
        j = 3*(i-1) + 1
        R[j], C[j], V[j] = i, i, -2.0/Δr^2

        # diagonal - 1
        j = 3*(i-1) + 2
        R[j], C[j], V[j] = i, i+1, 1.0/Δr^2

    end
    # Zero derivative at boundaries (Neumann)
    j= 1
    R[j], C[j], V[j] = 1, 1, 1.
    j = 2
    R[j], C[j], V[j] = 1, 2, -1.

    j = 3N - 2
    R[j], C[j], V[j] = N, N, 1.
    j = 3N - 3
    R[j], C[j], V[j] = N, N-1, -1.

    A = sparse(R,C,V)
end
function by_DEO()
    Δ = CenteredDifference(2,2,Δr,N)
    bc = Neumann0BC(Δr,1)
    A,b = sparse(Δ*bc)
    return A
end

println("Create by hand")
@btime Abh = by_hand();
println("Create by DiffEqOperators")
@btime Abd = by_DEO();

```

Output:

```julia
Create by hand
  6.089 ms (185603 allocations: 5.27 MiB)
Create by DiffEqOperators
  5.218 s (257528 allocations: 6.01 GiB)

```

I would have thought `sparse` on a `DerivativeOperator` or `GhostDerivativeOperator` would have done something like my `by_hand` function. Any thoughts?

---

<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:** [May 13, 2021, 5:47am UTC](https://discourse.julialang.org/t/solving-linear-system-with-summed-diffeqoperators/60998/7 "2021-05-13T05:47:28Z")

</div>

> [@bclyons12](#):
>
> I would have thought `sparse` on a `DerivativeOperator` or `GhostDerivativeOperator` would have done something like my `by_hand` function. Any thoughts?

At it to the issue. It should be doing it like you did by hand, it may be falling back to an iteration on the indices method.

---

<div class="post-metadata">

**Author:** ![bclyons12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bclyons12/32/33734_2.png) [@bclyons12](https://discourse.julialang.org/u/bclyons12)\
**Post date:** [May 13, 2021, 4:20pm UTC](https://discourse.julialang.org/t/solving-linear-system-with-summed-diffeqoperators/60998/8 "2021-05-13T16:20:23Z")

</div>

I opened a separate issue on git

[https://github.com/SciML/DiffEqOperators.jl/issues/392](https://github.com/SciML/DiffEqOperators.jl/issues/392)