How to use ApproxFun to multiply differential operators without running into incompatible space issues?

Specific Domains Numerics
1

I’m trying to compute (D + x)^2, where D is the derivative operator.

julia> x = Fun(identity)
Fun(Chebyshev(),[0.0, 1.0])

julia> D = Derivative(Chebyshev())
ConcreteDerivative : Chebyshev() → Ultraspherical(1)
 ⋅  1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   ⋅
 ⋅   ⋅   2.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   ⋅
 ⋅   ⋅    ⋅   3.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   ⋅
 ⋅   ⋅    ⋅    ⋅   4.0   ⋅    ⋅    ⋅    ⋅    ⋅   ⋅
 ⋅   ⋅    ⋅    ⋅    ⋅   5.0   ⋅    ⋅    ⋅    ⋅   ⋅
 ⋅   ⋅    ⋅    ⋅    ⋅    ⋅   6.0   ⋅    ⋅    ⋅   ⋅
 ⋅   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   7.0   ⋅    ⋅   ⋅
 ⋅   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   8.0   ⋅   ⋅
 ⋅   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   9.0  ⋅
 ⋅   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   ⋱
 ⋅   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   ⋱

julia> D2 = D + x;

julia> D2^2
ERROR: AssertionError: spacescompatible(domainspace(ops[k]), rangespace(ops[k + 1]))

Directly evaluating D2^2 errors, seemingly because the range and domain spaces of D2 differ. However, the conversion from the range space (Ultraspherical(1), or ChebyshevU) to the domain space (Chebyshev) doesn’t appear to not be defined, so this doesn’t work either:

julia> Conversion(rangespace(D2), domainspace(D2))
ERROR: Not implemented
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:33
 [2] Conversion(A::Ultraspherical{Int64, ChebyshevInterval{Float64}, Float64}, B::Chebyshev{ChebyshevInterval{Float64}, Float64})
   @ ApproxFunOrthogonalPolynomials ~/.julia/packages/ApproxFunOrthogonalPolynomials/QcyFy/src/Spaces/Ultraspherical/UltrasphericalOperators.jl:132
 [3] top-level scope
   @ REPL[127]:1

Is there a way to get around this, possibly using some other basis, for which the conversion would be defined? I tried using Ultraspherical(1/2) (Legendre) as the domain space instead of Chebyshev, but that doesn’t seem to work either:

julia> x = Fun(identity, S)
Fun(Ultraspherical(0.5),[-0.0, 1.0, -0.0, -0.0])

julia> D = Derivative(S)
ConcreteDerivative : Ultraspherical(0.5) → Ultraspherical(1.5)
 ⋅  1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   ⋅
 ⋅   ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   ⋅
 ⋅   ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   ⋅
 ⋅   ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅   ⋅
 ⋅   ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅   ⋅
 ⋅   ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅   ⋅
 ⋅   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅   ⋅
 ⋅   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅   ⋅
 ⋅   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0  ⋅
 ⋅   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   ⋱
 ⋅   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   ⋱

julia> Conversion(rangespace(D2), domainspace(D2))
ERROR: Not implemented
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:33
 [2] Conversion(A::Ultraspherical{Int64, ChebyshevInterval{Float64}, Float64}, B::Chebyshev{ChebyshevInterval{Float64}, Float64})
   @ ApproxFunOrthogonalPolynomials ~/.julia/packages/ApproxFunOrthogonalPolynomials/QcyFy/src/Spaces/Ultraspherical/UltrasphericalOperators.jl:132
 [3] top-level scope
   @ REPL[132]:1

Maybe I’m doing something wrong here, but my eventual goal is to obtain the operator

D2^2 = D^2 + 2xD + x^2 + 1
1 Like
2

Create the operator with no spaces and then set the space at the end:

julia> D = Derivative();

julia> x = Fun();

julia> D^2 + 2x*D + x^2 + 1 : Chebyshev()
PlusOperator : Chebyshev() → Ultraspherical(2)
 1.1666666666666667   0.0      4.625  …    ⋅                    ⋅
 0.0                  0.84375  0.0         ⋅                    ⋅
 0.08333333333333333  0.0      0.85        ⋅                    ⋅
  ⋅                   0.03125  0.0        0.020833333333333332  ⋅
  ⋅                    ⋅       0.025      0.0                   ⋱
  ⋅                    ⋅        ⋅     …  -1.0677083333333333    ⋱
  ⋅                    ⋅        ⋅         0.0                   ⋱
  ⋅                    ⋅        ⋅        18.084375              ⋱
  ⋅                    ⋅        ⋅         0.0                   ⋱
  ⋅                    ⋅        ⋅         0.9520833333333334    ⋱
  ⋅                    ⋅        ⋅     …    ⋱                    ⋱
2 Likes
3

Thanks a lot for the swift reply! That’s very helpful. Just wanted to post that the following works as well, with the derivatives computed automatically:

julia> D = Derivative()
ConcreteDerivative : ApproxFunBase.UnsetSpace() → ApproxFunBase.UnsetSpace()

julia> D + x
PlusOperator : ApproxFunBase.UnsetSpace() → ApproxFunBase.UnsetSpace()


julia> (D + x)^2
TimesOperator : ApproxFunBase.UnsetSpace() → ApproxFunBase.UnsetSpace()


julia> (D + x)^2 : Chebyshev()
TimesOperator : Chebyshev() → Ultraspherical(2)
 1.1666666666666667   0.0      4.625               0.0                   …    ⋅                   ⋅                    ⋅                    ⋅
 0.0                  0.84375  0.0                 6.28125                   0.03125              ⋅                    ⋅                    ⋅
 0.08333333333333333  0.0      0.8499999999999999  0.0                       0.0                 0.025                 ⋅                    ⋅
  ⋅                   0.03125  0.0                 0.8854166666666666       -1.09375             0.0                  0.020833333333333332  ⋅
  ⋅                    ⋅       0.025               0.0                       0.0                -1.0785714285714285   0.0                   ⋱
  ⋅                    ⋅        ⋅                  0.020833333333333332  …  14.109375000000002   0.0                 -1.0677083333333335    ⋱
  ⋅                    ⋅        ⋅                   ⋅                        0.0                16.095238095238095    0.0                   ⋱
  ⋅                    ⋅        ⋅                   ⋅                        0.940625            0.0                 18.084375              ⋱
  ⋅                    ⋅        ⋅                   ⋅                        0.0                 0.9469696969696969   0.0                   ⋱
  ⋅                    ⋅        ⋅                   ⋅                        0.0125              0.0                  0.9520833333333334    ⋱
  ⋅                    ⋅        ⋅                   ⋅                    …    ⋅                   ⋱                    ⋱                    ⋱
1 Like