# Two questions related to using operators in ApproxFun.jl

**URL:** <https://discourse.julialang.org/t/two-questions-related-to-using-operators-in-approxfun-jl/79831>\
**Category:** Numerics\
**Tags:** question, approxfun\
**Created:** [April 22, 2022, 11:15am UTC](https://discourse.julialang.org/t/two-questions-related-to-using-operators-in-approxfun-jl/79831 "2022-04-22T11:15:24Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![andreasvarga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasvarga/32/11634_2.png) [@andreasvarga](https://discourse.julialang.org/u/andreasvarga)\
**Post date:** [April 22, 2022, 11:15am UTC](https://discourse.julialang.org/t/two-questions-related-to-using-operators-in-approxfun-jl/79831/1 "2022-04-22T11:15:24Z")

</div>

I am about to register [PeriodicSystems.jl](https://github.com/andreasvarga/PeriodicSystems.jl), a new package dedicated to handle linear control systems with periodic time-varying coefficient matrices. For stability analysis, the computation of characteristic exponents of the state matrix A(t) is a possible approach. I already implemented a function `pseig`, which computes the characteristic exponents as the logarithms of the eigenvalues of the monodromy matrix (i.e., the state-transition matrix over one period). The following example illustrates the associated computational difficulties and the way to overcome them.

Consider the following periodic matrix of period T = 2pi  
`A(t) = [0 1; -10*cos(t) -24-10*sin(t)]`  
which has characterisitic exponents equal to 0 and -24 (you have to believe me).

However by computing the monodromy matrix and its eigenvalues, one of the resulting characteristic exponents has no **one** exact digit using Float64 based computations with high relative tolerance set for solving the underlying differential equations :

```julia
julia> using PeriodicSystems

julia> A = PeriodicFunctionMatrix(t -> [0 1; -10*cos(t) -24-10*sin(t)],2pi);

julia> @time ev = psceig(A; reltol = 1.e-10)
  0.004462 seconds (21.75 k allocations: 924.010 KiB, 79.31% compilation time)
2-element Vector{Float64}:
 -Inf
  -1.4052750767605842e-13

```

However, full accuracy can be achieved by determining the mondromy matrix as a product of say 500 matrices and computing the eigenvalues using the periodic Schur decomposition:

```julia
julia> @time ev = psceig(A,500; reltol = 1.e-10)
  0.029436 seconds (329.04 k allocations: 21.171 MiB)
2-element Vector{Float64}:
   1.8376538159566882e-15
 -23.99999999997412

```

Note that the evaluation of the 500 factors can be done in parallel, in which case a substantial speedup of computations can be achieved.

I recently included the [ApproxFun.jl](https://github.com/JuliaApproximation/ApproxFun.jl) package as an alternative way to define periodic matrices and performed the same computation using operators based computations. I must say, I was pleasantly surprized to obtain almost the same accuracy as before:

```julia
s = Fourier(0..2π);
a = Fun(t -> [0 1; -10*cos(t) -24-10*sin(t)],s);
d = domain(a);
D = Derivative(d);
ND = [D 0I; 0I D];
A = a-ND;
N = 50; n = (2*N+1)*2;
@time ev = eigvals(Matrix(A[1:n,1:n]));
  0.044203 seconds (88.73 k allocations: 4.403 MiB)
ev[sortperm(imag(ev),by=abs)][1:4]
4-element Vector{ComplexF64}:
    -24.000000000000135 + 0.0im
    -23.999999585663453 + 0.0im
 -6.042828006006208e-14 + 0.0im
                    0.0 + 0.0im

```

I guess that the operator A is a Fredholm operator (i.e., sum of a diagonal ND with a Toeplitz operator, but I am not sure) and I only displayed the eigenvalues in the central strip (i.e., in this case, those with null imaginary parts) using a truncated 2\*N+1 blocks representation which is typical for Toeplitz operators (not sure if this is the best choice).

My questions are:

1. What is the explanation for obtaining the characteristic exponents with double multiplicities?
2. There exists a nice compact way to build the diagonal operator ND for an arbitrary order A(t)?

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [April 22, 2022, 5:55pm UTC](https://discourse.julialang.org/t/two-questions-related-to-using-operators-in-approxfun-jl/79831/2 "2022-04-22T17:55:45Z")

</div>

1. Probably you could look at the eigenvectors to figure out why they are repeated.

2. No I don’t think so.

---

<div class="post-metadata">

**Author:** ![tobydriscoll](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tobydriscoll/32/1843_2.png) [@tobydriscoll](https://discourse.julialang.org/u/tobydriscoll)\
**Post date:** [April 22, 2022, 6:02pm UTC](https://discourse.julialang.org/t/two-questions-related-to-using-operators-in-approxfun-jl/79831/3 "2022-04-22T18:02:37Z")

</div>

FWIW, I tried running the problem in Chebfun, which is ApproxFun’s favorite Matlab uncle. The eigenvalues it gives me are

```julia
 -0.000000000000067 + 0.000000000000000i
 -0.000000000000067 - 1.000000000000009i
 -0.000000000000067 + 1.000000000000009i
 -0.000000000000064 - 2.000000000000016i
 -0.000000000000064 + 2.000000000000016i
 -0.000000000000067 - 3.000000000000025i
 -0.000000000000067 + 3.000000000000025i
...
 -0.000000000000005 -23.000000000000210i
 -0.000000000000005 +23.000000000000210i
 -0.000000000000007 -24.000000000000234i
 -0.000000000000007 +24.000000000000234i
 24.000000047940063 + 0.000000000000000i
 23.999999974147293 - 1.000000041935243i
 23.999999974147293 + 1.000000041935243i
 23.999999978621233 - 1.999999919103405i
 23.999999978621233 + 1.999999919103405i
 24.000000054662237 - 3.000000104379710i
 24.000000054662237 + 3.000000104379710i
 23.999999950623149 - 3.999999901300158i
 23.999999950623149 + 3.999999901300158i
 24.000000012168623 - 5.000000066251817i
 24.000000012168623 + 5.000000066251817i

```

Is that any less surprising?

---

<div class="post-metadata">

**Author:** ![tobydriscoll](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tobydriscoll/32/1843_2.png) [@tobydriscoll](https://discourse.julialang.org/u/tobydriscoll)\
**Post date:** [April 22, 2022, 6:05pm UTC](https://discourse.julialang.org/t/two-questions-related-to-using-operators-in-approxfun-jl/79831/4 "2022-04-22T18:05:47Z")

</div>

To clarify, I actually used the negative of the operator in the original code snippet.

---

<div class="post-metadata">

**Author:** ![andreasvarga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasvarga/32/11634_2.png) [@andreasvarga](https://discourse.julialang.org/u/andreasvarga)\
**Post date:** [April 22, 2022, 8:23pm UTC](https://discourse.julialang.org/t/two-questions-related-to-using-operators-in-approxfun-jl/79831/5 "2022-04-22T20:23:00Z")

</div>

The computed values are exactly what I would expect to obtain and their accuracy is indeed completely satisfactory. I suppose there are no double eigenvalues at -24 and 0 (it is not possible to check this from your result).  
Apparently you used another value of N (probably less than 50).

For comparison, these are the actually computed eigenvalues (sorted acording to the magnitudes of imaginary parts):

```julia
202-element Vector{ComplexF64}:
     -24.000000000000135 + 0.0im
     -23.999999585663453 + 0.0im
  -6.042828006006208e-14 + 0.0im
                     0.0 + 0.0im
     -24.000000393549527 - 0.9999999987185656im
     -24.000000393549527 + 0.9999999987185656im
  1.4765966227514582e-14 - 1.000000000000013im
  1.4765966227514582e-14 + 1.000000000000013im
  1.0936130473426786e-14 - 2.0000000000000027im
  1.0936130473426786e-14 + 2.0000000000000027im
      -23.99999966468944 - 2.0000000019980897im
      -23.99999966468944 + 2.0000000019980897im
                         ⋮
     -15.976780830674567 + 46.5840563132711im
      -32.02321916932537 - 46.58405631327116im
      -32.02321916932537 + 46.58405631327116im
 -1.7763568394002505e-14 - 46.99999999999988im
 -1.7763568394002505e-14 + 46.99999999999988im
  1.2434497884199578e-14 - 48.00000000000027im
  1.2434497884199578e-14 + 48.00000000000027im
  -2.842170943040401e-14 - 49.00000000000017im
  -2.842170943040401e-14 + 49.00000000000017im
 -1.1102230246251565e-14 - 50.000000000000064im
 -1.1102230246251565e-14 + 50.000000000000064im

```

---

<div class="post-metadata">

**Author:** ![tobydriscoll](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tobydriscoll/32/1843_2.png) [@tobydriscoll](https://discourse.julialang.org/u/tobydriscoll)\
**Post date:** [April 25, 2022, 1:47pm UTC](https://discourse.julialang.org/t/two-questions-related-to-using-operators-in-approxfun-jl/79831/6 "2022-04-25T13:47:01Z")

</div>

Chebfun chooses the discretization size automatically. The eigenfunctions are apparently discretized satisfactorily at 85 points, but that’s with a Chebyshev rather than Fourier method, which might use a factor of 2/pi fewer at N=54 or so.

The eigenfunction at -24 is low-frequency and therefore well resolved.

 ![eigenfunction](https://global.discourse-cdn.com/julialang/original/3X/7/b/7b335f91f5c2fbe9ae5e71eeb23d448acf7d6feb.png)

---

<div class="post-metadata">

**Author:** ![andreasvarga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasvarga/32/11634_2.png) [@andreasvarga](https://discourse.julialang.org/u/andreasvarga)\
**Post date:** [April 26, 2022, 3:21pm UTC](https://discourse.julialang.org/t/two-questions-related-to-using-operators-in-approxfun-jl/79831/7 "2022-04-26T15:21:30Z")

</div>

I am now able to compute the expected two real characteristic exponents without multiplicities. I employed another approach by generating explicitly the truncated Toeplitz - E(j0) operator in  
_J. Zhou, T. Hagiwara, and M. Araki, “Spectral characteristics and eigenvalues computation of the harmonic state operators in continuous-time periodic systems,” Systems and Control Letters, vol. 53, pp. 141–155, March 2004._  
and computing the eigenvalues. For N = 50, the resulting 8 eigenvalues with the smallest imaginary parts are:

```julia
8-element Vector{ComplexF64}:
   9.51898692030938e-15 - 1.5747448610226719e-16im
     -24.00000020745413 + 3.0152922283089535e-8im
     -23.99999976835794 + 0.9999999341699454im
 -4.068736769884833e-15 + 0.9999999999999764im
 3.4333666581810447e-15 - 1.000000000000011im
    -23.999999877051295 - 1.0000000414655406im
     -23.99999999433544 - 1.9999999296315698im
 1.2104134588879378e-14 + 1.999999999999935im

```

which do not exhibit any multiplicities of the first two eigenvalues in the central strip.

So, for my first question I still would like to understand what happens when using ApproxFun. Could it be that some boundary condition (e.g., addressing periodicity) should be added to the definition of the operator?

---

<div class="post-metadata">

**Author:** ![andreasvarga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasvarga/32/11634_2.png) [@andreasvarga](https://discourse.julialang.org/u/andreasvarga)\
**Post date:** [April 28, 2022, 8:42am UTC](https://discourse.julialang.org/t/two-questions-related-to-using-operators-in-approxfun-jl/79831/8 "2022-04-28T08:42:46Z")

</div>

I think I got it right now using

```julia
ev, V = ApproxFun.eigs(A,600,tolerance=1E-10);
ev[sortperm(imag(ev),by=abs)][1:4]
4-element Vector{ComplexF64}:
    -24.000000021801753 + 0.0im
 -6.440239494284886e-14 + 0.0im
    -24.000000080278454 + 0.9999997750716979im
     -24.00000008027845 - 0.9999997750716979im

```

The multiplicities are as expected, so the characteristic exponents are just the two eigenvalues in the central strip (with zero imaginary parts).

I wonder if it is possible to explicitly extract the finite matrix which underlies the eigenvalue computation with `eigs`. This would be usefull to build lifted (finite dimensional and constant) representations of periodic systems.

Regarding my second question, I could imagine to build a code generation sequence for arbitrary dimension `n` of `A(t)` to be evaluated, for example for `n = 2`, as

`eval(Meta.parse("ND = [D 0I; 0I D]"))`

but I wonder if the construction of `ND` would be efficient for large values of` n`. What would be a reasonable limit for `n`?

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [April 28, 2022, 9:45am UTC](https://discourse.julialang.org/t/two-questions-related-to-using-operators-in-approxfun-jl/79831/9 "2022-04-28T09:45:13Z")

</div>

I think it actually solves a generalised eigenvalue problem:

[https://github.com/JuliaApproximation/ApproxFun.jl/blob/999a138a0b0c02f780e26043eef738d03f50f256/src/Extras/eigs.jl#L17](https://github.com/JuliaApproximation/ApproxFun.jl/blob/999a138a0b0c02f780e26043eef738d03f50f256/src/Extras/eigs.jl#L17)

though for `Fourier` it might reduce to a standard eigenvalue problem. Maybe the key is `pruneeigs` which checks the eigenvectors for convergence?

---

<div class="post-metadata">

**Author:** ![andreasvarga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasvarga/32/11634_2.png) [@andreasvarga](https://discourse.julialang.org/u/andreasvarga)\
**Post date:** [April 28, 2022, 11:51am UTC](https://discourse.julialang.org/t/two-questions-related-to-using-operators-in-approxfun-jl/79831/10 "2022-04-28T11:51:26Z")

</div>

Yes, indeed, `pruneeigs` is used to dismiss two of the multiple eigenvalues of `A1` (`C1 = I`). I also found the explanations in [#645](https://github.com/JuliaApproximation/ApproxFun.jl/issues/645), so I think this aspect is now clear for me. Thanks.

---

<div class="post-metadata">

**Author:** ![andreasvarga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasvarga/32/11634_2.png) [@andreasvarga](https://discourse.julialang.org/u/andreasvarga)\
**Post date:** [April 30, 2022, 12:32pm UTC](https://discourse.julialang.org/t/two-questions-related-to-using-operators-in-approxfun-jl/79831/11 "2022-04-30T12:32:01Z")

</div>

My last question: It is possible to obtain the same results using Chebysev functions with a kind of periodic boundary conditions?