Operators

Linear operators between two spaces in ApproxFun are represented by subtypes of Operator. Every operator has a domainspace and rangespace. That is, if a Fun f has the space domainspace(op), thenop*f is a Fun with space rangespace(op).

Note that the size of an operator is specified by the dimension of the domain and range space.

Calculus operators

Differential and integral operators are perhaps the most useful type of operators in mathematics. Consider the derivative operator on CosSpace:

julia> D = Derivative(CosSpace())ConcreteDerivative : CosSpace(【0.0,6.283185307179586❫) → SinSpace(【0.0,6.283185307179586❫) 0.0  -1.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅   ⋅
  ⋅    0.0  -2.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅   ⋅
  ⋅     ⋅    0.0  -3.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅   ⋅
  ⋅     ⋅     ⋅    0.0  -4.0    ⋅     ⋅     ⋅     ⋅     ⋅   ⋅
  ⋅     ⋅     ⋅     ⋅    0.0  -5.0    ⋅     ⋅     ⋅     ⋅   ⋅
  ⋅     ⋅     ⋅     ⋅     ⋅    0.0  -6.0    ⋅     ⋅     ⋅   ⋅
  ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    0.0  -7.0    ⋅     ⋅   ⋅
  ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    0.0  -8.0    ⋅   ⋅
  ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    0.0  -9.0  ⋅
  ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    0.0  ⋱
  ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅   ⋱
julia> f = Fun(θ->cos(cos(θ)), CosSpace());
julia> fp = D*f;
julia> fp(0.1) ≈ f'(0.1) ≈ sin(cos(0.1))*sin(0.1)true

Here, we specified the domain space for the derivative operator, and it automatically determined the range space:

julia> rangespace(D) == space(fp) == SinSpace()true

Operators can be identified with infinite-dimensional matrices, whose entries are given by the canonical bases in the domain and range space. In this case, the relevant formula is

\[\mathop{D} \cos{kθ} = -k \sin{kθ}.\]

That is, the (k,k+1)th entry is as follows:

julia> k,j = 5,6;
julia> ej = Fun(domainspace(D),[zeros(j-1);1]);
julia> D[k,j] ≈ (D*ej).coefficients[k] ≈ -ktrue

The Chebyshev space has the property that its derivatives are given by ultraspherical spaces:

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

Functionals

A particularly useful class of operators are functionals, which map from functions to scalar numbers. These are represented by operators of size 1 × ∞: that is, infinite-dimensional analogues of row vectors.

As an example, the evaluation functional f(0) on CosSpace has the form:

julia> B = Evaluation(CosSpace(),0)ConcreteEvaluation : CosSpace(【0.0,6.283185307179586❫) → ConstantSpace(Point(0)) 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  ⋯
julia> B*f ≈ f(0)true

As can be seen from the output, rangespace(B) is a ConstantSpace(Point(0)), a one-dimensional space used to represent scalars whose domain is a single point, 0.

Closely related to functionals are operators with finite-dimensional range. For example, the Dirichlet operator represents the restriction of a space to its boundary. In the case, of Chebyshev(), this amounts to evaluation at the endpoints ±1:

julia> B = Dirichlet(Chebyshev())ConcreteDirichlet : Chebyshev() → 2-element ArraySpace:
ConstantSpace{Point{Float64}, Float64}[ConstantSpace(Point(-1.0)), ConstantSpace(Point(1.0))] 1.0  -1.0  1.0  -1.0  1.0  -1.0  1.0  -1.0  1.0  -1.0  ⋯
 1.0   1.0  1.0   1.0  1.0   1.0  1.0   1.0  1.0   1.0  ⋯
julia> size(B)(2, ℵ₀)
julia> B*Fun(exp)Fun(2-element ArraySpace: ConstantSpace{Point{Float64}, Float64}[ConstantSpace(Point(-1.0)), ConstantSpace(Point(1.0))], [0.3678794411714403, 2.718281828459044])
julia> B*Fun(exp) ≈ Fun([exp(-1),exp(1)])true

Multiplication

A Multiplication operator sends a Fun to a Fun in the corresponding space by multiplying a given function. The Multiplication operators are presented in matrix form in ApproxFun.

julia> x = Fun();
julia> M = Multiplication(1 + 2x + x^2, Chebyshev())ConcreteMultiplication : Chebyshev() → Chebyshev() 1.5 1.0 0.25 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 2.0 1.75 1.0 0.25 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 1.0 1.5 1.0 0.25 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.25 1.0 1.5 1.0 0.25 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.25 1.0 1.5 1.0 0.25 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.25 1.0 1.5 1.0 0.25 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.25 1.0 1.5 1.0 0.25 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.25 1.0 1.5 1.0 0.25 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.25 1.0 1.5 1.0 ⋱ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.25 1.0 1.5 ⋱ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋱ ⋱ ⋱
julia> (M * x).coefficients == ((1 + 2x + x^2) * x).coefficients == M[1:4,1:2] * x.coefficientstrue

It is possible for domain space and range space to be different under Mulitplication.

julia> c = Fun(θ -> cos(θ), CosSpace());
julia> Multiplication(c, SinSpace())ConcreteMultiplication : SinSpace(【0.0,6.283185307179586❫) → SinSpace(【0.0,6.283185307179586❫) 8.974302782386682e-17 0.5 … ⋅ ⋅ 0.5 8.974302782386682e-17 ⋅ ⋅ ⋅ 0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 ⋅ ⋅ ⋅ 8.974302782386682e-17 ⋱ ⋅ ⋅ … ⋱ ⋱

If a function is given by the expansion

\[\mathop{f}(θ) = \sum_{n=1}^∞ f_n \sin{nθ}.\]

Then the matrix above can be easily derived from

\[\begin{aligned} \cos{θ} \mathop{f}(θ) &= \cos{θ} \sum_{n=1}^∞ f_n \sin{nθ} \\ &= \sum_{n=1}^∞ f_n \cos{θ} \sin{nθ} \\ &= \sum_{n=1}^∞ \frac{1}{2} f_n \left(\sin{(n-1)θ} + \sin{(n+1)θ}\right) \\ &= \sum_{n=1}^∞ \frac{1}{2} \left(f_{n-1} + f_{n+1}\right) \sin{nθ}, \end{aligned}\]

where $f_0 = 0$.

Algebraic manipulation of operators

Operators can be algebraically manipulated, provided that the domain and range spaces are compatible, or can be made compatible. As a simple example, we can add the second derivative of a Fourier space to the identity operator:

julia> D2 = Derivative(Fourier(),2)DerivativeWrapper : Fourier(【0.0,6.283185307179586❫) → Fourier(【0.0,6.283185307179586❫) 0.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅      ⋅      ⋅      ⋅   ⋅
  ⋅   -1.0    ⋅     ⋅     ⋅     ⋅     ⋅      ⋅      ⋅      ⋅   ⋅
  ⋅     ⋅   -1.0    ⋅     ⋅     ⋅     ⋅      ⋅      ⋅      ⋅   ⋅
  ⋅     ⋅     ⋅   -4.0    ⋅     ⋅     ⋅      ⋅      ⋅      ⋅   ⋅
  ⋅     ⋅     ⋅     ⋅   -4.0    ⋅     ⋅      ⋅      ⋅      ⋅   ⋅
  ⋅     ⋅     ⋅     ⋅     ⋅   -9.0    ⋅      ⋅      ⋅      ⋅   ⋅
  ⋅     ⋅     ⋅     ⋅     ⋅     ⋅   -9.0     ⋅      ⋅      ⋅   ⋅
  ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅   -16.0     ⋅      ⋅   ⋅
  ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅      ⋅   -16.0     ⋅   ⋅
  ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅      ⋅      ⋅   -25.0  ⋅
  ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅      ⋅      ⋅      ⋅   ⋱
julia> D2 + IPlusOperator : Fourier(【0.0,6.283185307179586❫) → Fourier(【0.0,6.283185307179586❫) 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -3.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -3.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -8.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -8.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -15.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -15.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -24.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋱

When the domain and range space are not the same, the identity operator becomes a conversion operator. That is, to represent D+I acting on the Chebyshev space, we would do the following:

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> C = Conversion(Chebyshev(),Ultraspherical(1))ConcreteConversion : Chebyshev() → Ultraspherical(1) 1.0 0.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 0.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 0.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 0.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 0.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 0.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 0.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 0.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 0.0 ⋱ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 ⋱ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋱
julia> D + CPlusOperator : Chebyshev() → Ultraspherical(1) 1.0 1.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 2.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 3.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 4.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 5.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 6.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 7.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 8.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 9.0 ⋱ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 ⋱ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋱

ApproxFun can automatically determine the spaces, so if one writes D + I it will translate it to D + C.

Now consider the Fredholm integral operator of the second kind:

\[\mathop{L} \mathop{u} = \mathop{u} + \mathop{e}^x \int_{-1}^1 \mathop{u}(x) \mathop{dx}.\]

We can construct this using

julia> x = Fun();
julia> Σ = DefiniteIntegral(Chebyshev())ConcreteDefiniteIntegral : Chebyshev() → ConstantSpace 2.0 0.0 -0.6666666666666666 0.0 … 0.0 -0.031746031746031744 0.0 ⋯
julia> L = I + exp(x)*ΣLowRankPertOperator : Chebyshev() → Chebyshev() 3.5321317555040164 0.0 … -0.0401925675476828 0.0 ⋯ 2.2606364159699406 1.0 -0.03588311771380858 0.0 ⋱ 0.5429906790681532 0.0 -0.008618899667748463 0.0 ⋱ 0.08867369969732763 0.0 -0.0014075190428147243 0.0 ⋱ 0.010948480884187472 0.0 -0.00017378541086011858 0.0 ⋱ 0.0010858526238278882 0.0 … -1.7235755933776e-5 0.0 ⋱ 8.995464590859037e-5 0.0 -1.4278515223585772e-6 0.0 ⋱ 6.396872924803985e-6 0.0 -1.0153766547307912e-7 0.0 ⋱ 3.9842496133455947e-7 0.0 0.9999999936757943 0.0 ⋱ 2.2073543451034698e-8 0.0 -3.5037370557197933e-10 1.0 ⋱ ⋮ ⋱ … ⋱ ⋱ ⋱
julia> u = cos(10x^2);
julia> (L*u)(0.1)1.3777980523127369
julia> u(0.1) + exp(0.1)*sum(u)1.377798052312737

Note that DefiniteIntegral is a functional, i.e., a 1 × ∞ operator. when multiplied on the left by a function, it automatically constructs the operator $\mathop{e}^x \int_{-1}^1 \mathop{f}(x) \mathop{dx}$ via

julia> M = Multiplication(exp(x),ConstantSpace())ConcreteMultiplication : ConstantSpace → Chebyshev() 1.2660658777520082
 1.1303182079849703
 0.2714953395340766
 0.04433684984866382
 0.005474240442093736
 0.0005429263119139441
 4.497732295429518e-5
 3.1984364624019926e-6
 1.9921248066727974e-7
 1.1036771725517349e-8
  ⋮
julia> M*ΣTimesOperator : Chebyshev() → Chebyshev() 2.5321317555040164 0.0 … -0.0401925675476828 0.0 ⋯ 2.2606364159699406 0.0 -0.03588311771380858 0.0 ⋱ 0.5429906790681532 0.0 -0.008618899667748463 0.0 ⋱ 0.08867369969732763 0.0 -0.0014075190428147243 0.0 ⋱ 0.010948480884187472 0.0 -0.00017378541086011858 0.0 ⋱ 0.0010858526238278882 0.0 … -1.7235755933776e-5 0.0 ⋱ 8.995464590859037e-5 0.0 -1.4278515223585772e-6 0.0 ⋱ 6.396872924803985e-6 0.0 -1.0153766547307912e-7 0.0 ⋱ 3.9842496133455947e-7 0.0 -6.324205735469197e-9 0.0 ⋱ 2.2073543451034698e-8 0.0 -3.5037370557197933e-10 0.0 ⋱ ⋮ ⋱ … ⋱ ⋱ ⋱

Note that Σ*exp(x) applies the operator to a function. To construct the operator that first multiplies by exp(x), use Σ[exp(x)]. This is equivalent to Σ*Multiplication(exp(x),Chebyshev()).

Operators and space promotion

It is often more convenient to not specify a space explicitly, but rather infer it when the operator is used. For example, we can construct Derivative(), which has the alias 𝒟, and represents the first derivative on any space:

julia> f = Fun(cos,Chebyshev(0..1)); (𝒟*f)(0.1)-0.09983341664681751
julia> f = Fun(cos,Fourier()); (𝒟*f)(0.1)-0.09983341664682804

Behind the scenes, Derivative() is equivalent to Derivative(UnsetSpace(),1). When multiplying a function f, the domain space is promoted before multiplying, that is, Derivative()*f is equivalent to Derivative(space(f))*f.

This promotion of the domain space happens even when operators have spaces attached. This facilitates the following construction:

julia> D = Derivative(Chebyshev());
julia> D^2ConcreteDerivative : Chebyshev() → Ultraspherical(2) ⋅ ⋅ 4.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 6.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 8.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 10.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 12.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 14.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 16.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 18.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋱ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋱ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋱

Note that rangespace(D) ≠ Chebyshev(), hence the operators are not compatible. Therefore, it has thrown away its domain space, and thus this is equivalent to Derivative(rangespace(D))*D.

Concatenating operators

The concatenation functions vcat, hcat and hvcat are overriden for operators to represent the resulting combined operator, now with a rangespace or domainspace that is an ArraySpace.