# 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.3777980523127369julia> 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.09983341664681751julia> 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.