Solving the Helmholtz eq. with the Neumann condition generates `DimensionMismatch' error

[BoundaryValueProblems]

The following code segment worked sometime ago, but now with the current version of ApproxFun with Julia v1.0.1, I got the DimensionMismatch error as follows. Something must have been changed.

julia> Ω=Interval()^2  #The domain is the square [-1,1]^2.
ProductDomain{Tuple{Segment{Float64},Segment{Float64}},StaticArrays.SArray{Tuple{2},Float64,1,2}}((【-1.0,1.0】, 【-1.0,1.0】))
julia> Δ=Laplacian(Ω)  #This defines the Laplace operator on Ω.
LaplacianWrapper:Chebyshev(【-1.0,1.0】)⊗Chebyshev(【-1.0,1.0】)→Ultraspherical(2,【-1.0,1.0】)⊗Ultraspherical(2,【-1.0,1.0】)
 0.0  0.0  0.0  4.0  0.0  4.0  0.0  0.0  0.0  0.0  ⋯
 0.0  0.0  0.0  0.0  0.0  0.0  6.0  0.0  1.0  0.0  ⋱
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  6.0  ⋱
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
  ⋮    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱   ⋱
julia> B=Neumann(Ω)    #This defines the Neumann boundary operator on Ω.
DirichletWrapper:Chebyshev(【-1.0,1.0】)⊗Chebyshev(【-1.0,1.0】)→2-element ArraySpace:
TensorSpace{SV,D,Array{Float64,1}} where D where SV[2-element ArraySpace:
ConstantSpace{ApproxFun.Point{Float64},Float64}[ConstantSpace(Point(-1.0)), ConstantSpace(Point(1.0))]⊗Chebyshev(【-1.0,1.0】), Chebyshev(【-1.0,1.0】)⊗2-element ArraySpace:
ConstantSpace{ApproxFun.Point{Float64},Float64}[ConstantSpace(Point(-1.0)), ConstantSpace(Point(1.0))]]
 0.0  0.0  1.0   0.0  0.0  -4.0  0.0   0.0   0.0  9.0  ⋯
 0.0  0.0  1.0   0.0  0.0   4.0  0.0   0.0   0.0  9.0  ⋱
 0.0  1.0  0.0  -4.0  0.0   0.0  9.0   0.0   0.0  0.0  ⋱
 0.0  1.0  0.0   4.0  0.0   0.0  9.0   0.0   0.0  0.0  ⋱
 0.0  0.0  0.0   0.0  1.0   0.0  0.0   0.0  -4.0  0.0  ⋱
 0.0  0.0  0.0   0.0  1.0   0.0  0.0   0.0   4.0  0.0  ⋱
 0.0  0.0  0.0   0.0  1.0   0.0  0.0  -4.0   0.0  0.0  ⋱
 0.0  0.0  0.0   0.0  1.0   0.0  0.0   4.0   0.0  0.0  ⋱
 0.0  0.0  0.0   0.0  0.0   0.0  0.0   1.0   0.0  0.0  ⋱
 0.0  0.0  0.0   0.0  0.0   0.0  0.0   1.0   0.0  0.0  ⋱
  ⋮    ⋱    ⋱     ⋱    ⋱     ⋱    ⋱     ⋱     ⋱    ⋱   ⋱
julia> g=Fun((x,y)->sin(5π*x-7π*y),Ω) #A trial 2D sinusoid
Fun(Chebyshev(【-1.0,1.0】)⊗Chebyshev(【-1.0,1.0】),[-1.29382e-17, 0.0334004, -0.0332574, -2.18927e-18, -3.50038e-17, 4.329e-18, 0.0267048, -0.0724972, 0.0751763, -0.0235788  …  9.5227e-15, -1.01484e-14, 1.01389e-14, -9.48102e-15, 8.29904e-15, -6.84613e-15, 5.26785e-15, -3.77486e-15, 2.54284e-15, -1.61381e-15])
julia> f=B*g           #This computes the Neumann boundary condition from g.
Fun(2-element ArraySpace:
TensorSpace{SV,D,Array{Float64,1}} where D where SV[2-element ArraySpace:
ConstantSpace{ApproxFun.Point{Float64},Float64}[ConstantSpace(Point(-1.0)), ConstantSpace(Point(1.0))]⊗Chebyshev(【-1.0,1.0】), Chebyshev(【-1.0,1.0】)⊗2-element ArraySpace:
ConstantSpace{ApproxFun.Point{Float64},Float64}[ConstantSpace(Point(-1.0)), ConstantSpace(Point(1.0))]],[1.87882, 1.87882, -3.10476, -3.10476, 9.3604e-13, 5.07142e-13, -2.34628e-13, 7.83046e-13, 4.09561, 4.09561  …  1.77235e-16, -2.60956e-16, 4.62124e-17, -6.00936e-17, 7.71524e-17, 7.71524e-17, 4.06824e-17, 4.06824e-17, 0.0, 0.0])
julia> k2=(5π)^2+(7π)^2-0.001  #Slightly perturbed k^2
730.3497256806126
julia> u=\([B;Δ+k2*I],[f;0.0];tolerance=1E-12)  #The solution!
ERROR: LoadError: DimensionMismatch("dimensions must match")
Stacktrace:
 [1] coefficients(::Array{Float64,1}, ::ApproxFun.ArraySpace{Space,1,ProductDomain{Tuple{ApproxFun.Point{Float64},Segment{Float64}},StaticArrays.SArray{Tuple{2},Float64,1,2}},Any}, ::ApproxFun.ArraySpace{Space,1,ProductDomain{Tuple{ApproxFun.Point{Float64},Segment{Float64}},StaticArrays.SArray{Tuple{2},Float64,1,2}},Any}) at /Users/xxx/.julia/packages/ApproxFun/wFAIT/src/Spaces/Modifier/ArraySpace.jl:202
 [2] coefficients at /Users/xxx/.julia/packages/ApproxFun/wFAIT/src/Fun/Fun.jl:45 [inlined]
 [3] #\#183(::Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol},NamedTuple{(:tolerance,),Tuple{Float64}}}, ::Function, ::ApproxFun.InterlaceOperator{Float64,1,TensorSpace{Tuple{Chebyshev{Segment{Float64},Float64},Chebyshev{Segment{Float64},Float64}},ProductDomain{Tuple{Segment{Float64},Segment{Float64}},StaticArrays.SArray{Tuple{2},Float64,1,2}},Float64},ApproxFun.ArraySpace{Space,1,ProductDomain{Tuple{ApproxFun.Point{Float64},Segment{Float64}},StaticArrays.SArray{Tuple{2},Float64,1,2}},Any},ApproxFun.CachedIterator{Tuple{Int64,Int64},ApproxFun.BlockInterlacer{Tuple{ApproxFun.UnitCount{Int64}}}},ApproxFun.CachedIterator{Tuple{Int64,Int64},ApproxFun.BlockInterlacer{Tuple{ApproxFun.Flatten{Tuple{Array{Int64,1},ApproxFun.Repeated{Int64}}},ApproxFun.UnitCount{Int64}}}},Tuple{ApproxFun.Infinity{Bool},ApproxFun.Infinity{Bool}}}, ::Fun{ApproxFun.ArraySpace{Space,1,ProductDomain{Tuple{ApproxFun.Point{Float64},Segment{Float64}},StaticArrays.SArray{Tuple{2},Float64,1,2}},Any},Float64,Array{Float64,1}}) at /Users/xxx/.julia/packages/ApproxFun/wFAIT/src/Operators/ldiv.jl:11
 [4] (::getfield(Base, Symbol("#kw##\")))(::NamedTuple{(:tolerance,),Tuple{Float64}}, ::typeof(\), ::ApproxFun.InterlaceOperator{Float64,1,TensorSpace{Tuple{Chebyshev{Segment{Float64},Float64},Chebyshev{Segment{Float64},Float64}},ProductDomain{Tuple{Segment{Float64},Segment{Float64}},StaticArrays.SArray{Tuple{2},Float64,1,2}},Float64},ApproxFun.ArraySpace{Space,1,ProductDomain{Tuple{ApproxFun.Point{Float64},Segment{Float64}},StaticArrays.SArray{Tuple{2},Float64,1,2}},Any},ApproxFun.CachedIterator{Tuple{Int64,Int64},ApproxFun.BlockInterlacer{Tuple{ApproxFun.UnitCount{Int64}}}},ApproxFun.CachedIterator{Tuple{Int64,Int64},ApproxFun.BlockInterlacer{Tuple{ApproxFun.Flatten{Tuple{Array{Int64,1},ApproxFun.Repeated{Int64}}},ApproxFun.UnitCount{Int64}}}},Tuple{ApproxFun.Infinity{Bool},ApproxFun.Infinity{Bool}}}, ::Fun{ApproxFun.ArraySpace{Space,1,ProductDomain{Tuple{ApproxFun.Point{Float64},Segment{Float64}},StaticArrays.SArray{Tuple{2},Float64,1,2}},Any},Float64,Array{Float64,1}}) at ./none:0
 [5] top-level scope at none:0

[dlfivefifty]

Ah, it was because you were using concatenation (which combines vector-valued functions) when you really want

u=\([B;Δ+k2*I],[f,0.0];tolerance=1E-12) 

[BoundaryValueProblems]

@dlfivefifty: Thanks a lot! It works now. But could you explain to me the difference of this [f,0.0] for this Neumann problem as I posed earlier and the Dirichlet problem in README.md that uses [f;0.0]? Why does the concatenation ; work for the latter case? Is this because f in the latter is not a function but simply a vector? Thanks!

[dlfivefifty]

Yes that’s right: there’s an inconsistency with Dirichlet and Neumann conditions. Direchlet are treating as scalar valued functions defined on the boundary. Neumann conditions are (due to older code) treated as vector-valued functions.

[BoundaryValueProblems]

@dlfivefifty: Thanks for your clarification. Then, I have another question for specifying the boundary condition. In this case, the Dirichlet boundary condition. Consider the following codes:

Ω=Interval()^2  #The domain is the square [-1,1]^2.
Δ=Laplacian(Ω)  #This defines the Laplace operator on Ω.
Bd=Dirichlet(Ω)    #This defines the Dirichlet boundary operator on Ω.
g=Fun((x,y)->sin(5π*x-7π*y),Ω) #A trial 2D sinusoid
f=Fun((x,y)->g(x,y),∂(Ω))  #This computes the Dirichlet boundary condition from g.
f2=Bd*g # This should compute the same Dirichlet boundary condition from g.
k2=(5π)^2+(7π)^2-0.001  #Slightly perturbed k^2
u=\([Bd;Δ+k2*I],[f;0.0];tolerance=1E-12)  #The solution!
u2=\([Bd;Δ+k2*I],[f2;0.0];tolerance=1E-12)  #The solution too!

Then, let's check the relative L^2 error:

julia> norm(u-u2)/norm(u)
1.928042887351262e-12

Of course, I really shouldn't complain since these solutions have error tolerance of 1E-12 as I specified in solving them. Yet, I want to know what's the real difference between:

f=Fun((x,y)->g(x,y),∂(Ω))  #This computes the Dirichlet boundary condition from g.
f2=Bd*g # This should compute the same Dirichlet boundary condition from g.

Is this difference between function evaluation vs matrix multiplication? Thanks a lot!

[dlfivefifty]

That’s right: Bd*g creates a (sparse) matrix to multiply the coefficients while the other version uses transforms

[BoundaryValueProblems]

Thanks a lot, @dlfivefifty!