Eigen value and eigen functions of linear operator

Specific Domains Numerics
, , ,
1

I am doing the linear stability analysis of a system of equations representing a free surface flow, which typically ends up with a generalised eigenvalue problem of the form:

  • L_{1} X = \lambda p(z) u
  • L_{2} X = \lambda p(z) w

where L_{i} are linear differential operators involving u, w, \Gamma, \Gamma_{z}, \Omega and \Omega_{z}. X = [u, w, \Gamma, \Omega] is a vector of variables, p is a known function, and \lambda the eigenvalue. Two algebraic equations complement this system:

  • \Gamma = p(z)u_{z},
  • \Omega =p(z)w_{z},

I also have 5 boundary conditions:

  • u(0) = w(0) = 0
  • \Gamma(1) = -C_{1}h
  • \Omega(1) = -C_{2}h
  • w(1) = h(\lambda + C_{3})

where C_{i} are constants. Note that those BCs introduce a new variable h (which can be eliminated), and that the last one involves the eigenvalue. This is why those variables have been introduced.

An important precision: p(1) = 0. Then, if eliminating \Omega and \Gamma of the problem, the two related boundary conditions degenerate, so I want to allow for the product p(z)u_{z} and p(z)w_{z} to be regular.

I am new to these kinds of problems, but I understand that I have basically two dominant ways of approaching this:

  1. the spectral method, i.e Galerkin or Tau method
  2. the pseudospectral method, e.g. Chebyshev collocation method

I have several questions before committing to a way of solving this:

  • What Julia packages would provide a high-level interface for doing this? From what I have found:

    • ApproxFun.jl works in the spectral case and would provide support for the spectral method.
    • From this post, it seems that ClassicalOrthogalPolynomial.jl is the package to use for the collocation method, to obtain the discretisation of the operator on the collocation points?

    Are those the right packages for this, or are there other better-suited?

  • I am familiar with how to impose boundary conditions on the collocation method. However, I am not familiar with the spectral method and working in the coefficient space. ApproxFun.jl has examples about solving simpler eigenvalue problems with simpler boundary conditions using basis recombination. How should I handle my more complex boundary conditions using ApproxFun? I can see how I could build my matrix by stacking operators and the algebraic relations, but that’s as far as I know for now.

As a concluding remark, I am open to any suggestions for solving this problem.

Thanks!

2

It’s hard to answer since you don’t give much detail but probably ClassicalOrthogonalPolynomials.jl can solve the problem, either with collocation or with ultraspherical spectral method

1 Like
4

Thanks for confirming this!

  • I have added some information about the eigenvalue problem. What additional details would be useful for a more detailed answer?
  • It seems ApproxFun also uses the ultraspherical spectral method. Any reason you’d recommend ClassicalOrthogonalPolynomials.jl over ApproxFun.jl? Because it’s lower-level and therefore easier to adapt to a more complex problem?
  • I have read a bit more on the ultraspherical method, but I have only seen examples with a single equation/variable. In the case where I have multiple variables and equations, the operator matrix is presumably constructed by blocks from the operators corresponding to each equation. In this case, what lines do you remove? Do you still add the BCs lines at the top, to keep the matrix close to upper triangular (by convention, if I am not mistaken) ?
5

What additional details would be useful for a more detailed answer?

If there’s a simple MWE I can try to code up an example.

It seems ApproxFun also uses the ultraspherical spectral method. Any reason you’d recommend ClassicalOrthogonalPolynomials.jl over ApproxFun.jl?

ClassicalOrthogonalPolynomials.jl is designed to be more “low level” and versatile. And it’s possible to do collocation / Galerkin / ultraspherical spectral methods in a single language. It was designed specifically to avoid some weaknesses in ApproxFun’s design, and with the idea that ApproxFun will eventually sit on top of it (though that’s not any closer to reality than it was 3 years ago).

Though it’s still very experimental.

In this case, what lines do you remove? Do you still add the BCs lines at the top, to keep the matrix close to upper triangular (by convention, if I am not mistaken) ?

In ultraspherical spectral method you just remove the bottom lines and add the missing rows to the top. This is roughly equivalent to rectangular collocation. When the variable coefficients are “nice” you can maintain optimal complexity using SemiseparableMatrices.jl. But this is easier to explain via an example.

1 Like
7

Thanks for these answers.

I’ll cook a minimal example that retains the difficulties that I still have and come back to you.

I guess the main issues that remain for me are the multivariate aspect (most examples are given with a single variable and a single equation), and the inclusion of the boundary condition containing the eigenvalue with an additional variable.

EDIT: what about this “toy” problem?

- u''(z) - a(z)\, w(z) = \lambda b(z) u(z)
- w''(z) - u(z) = \lambda b(z) w(z)

with z \in [0,1], a(z) = 1+z and b(z) = z and the boundary conditions:

  • u(0) = 0
  • w(0) = 0
  • w'(1) + \lambda w(1) + u(1) = 0
8

There are also finite-element methods (and finite-differences etcetera)… lots of ways to discretize systems of PDEs, with different tradeoffs, and lots of books and courses about the topic.

1 Like
9

I did a quick implementation using an ultraspherical discretisation:

It seemed that you were missing a boundary condition since for a system of two 2nd order diff eqs one would expect 4 conditions. So I added (arbitrarily) u(0) = w(1). This produced a reasonable eigenfunction:

image
image1228×854 12.4 KB

But if you think it is indeed only 3 boundary conditions the code can be modified (though one would need to decide whether to drop two rows from u or w).

It’s also easy to modify this for code for collocation, Galerkin, or indeed FEM (including hp-FEM).

Though you mentioned multiple variables. PDEs is less developed. In theory it’s possible to do triangles, disk, and rectangles but there are a lot of missing features that you would need for complicated bcs.

Edit: Actually Galerkin and FEM are more complicated for adding general BCs so that would take a bit of work.

4 Likes
10

I was aiming at those spectral/pseudo-spectral methods as they seem to be the usual way for those eigenvalue problems resulting from hydrodynamic stability analysis, presumably thanks to their exponential convergence.

11

Thank you very much for this example. I should be able to work from this. I have a couple of additional questions, if you find the time at some point:

  • Indeed, it was a typo and I forgot to type in one Bcs.

But if you think it is indeed only 3 boundary conditions, the code can be modified (though one would need to decide whether to drop two rows from u or w).

Just to clarify this. Here, there are 4 BCs, so you symmetrically drop 2 rows from u and 2 rows from w. If there are only 3, you have to drop 1 row from one, and 2 rows from the other, right? Then the choice is which operator you’d drop 2 from, and which you’d drop 1 from?

12

If there are only 3, you have to drop 1 row from one, and 2 rows from the other, right? Then the choice is which operator you’d drop 2 from, and which you’d drop 1 from?

I think talking about what rows to drop in an ill-posed equation is not meaningful :sweat_smile:

Some equations naturally have only 3 bcs. For example, I think if one of the operators is degenerate:

u'' + … 
xw'' + …

Since x vanishes at zero its a singular point of the ODE. So in this case you probably only need 3 bcs, and what’s natural is to drop two rows from u and one row from w.

The τ method gives more flexible ways of “dropping rows” which become important in PDEs, where the question is quite nuanced, see

But I would probably suggest embracing FEM for PDEs as it naturally deals with corners, etc. Though you’d need to formulate your BCs in weak form which something I don’t know about (but I believe is standard and probably you can just ask AI these days).

How to reference ClassicalOrthogonalPolynomials.jl in a publication? Simply link to the GitHub repo?

You can just cite my Acta Numerica paper, I’ve added a CITATION.bib:

2 Likes
13

I believe

https://academic.oup.com/mnras/article/455/4/4274/1267403

is a good reference for why spectral methods are “better” than mesh-based methods for hydrodynamic stability.

Another relevant reference shows hp-FEM avoids pollution effects:

(I believe a consequence of this is that spectral methods also avoid pollution…)

1 Like
14

Thanks a lot for all these additional references!

15

I have played around with ClassicalOrthogonalPolynomials.jl and the ultraspherical method. I now want to compare the result with rectangular collocation. Is there a way to automatically deal with the mapping from [-1, 1] to [0, 1]:

T = chebyshevt(0 .. 1) # solution basis is T_n

when building the grid;

x1 = reverse(ChebyshevGrid{1}(n - 2)) # 1st kind points (Chebyshev-Gauss), sorted
x2 = reverse(ChebyshevGrid{2}(n)) # 2nd kind points (Chebyshev-Gauss-Lobatto), sorted

which are currently between -1 and 1 ?

16

For 1st kind points there’s a grid routine:

julia> using ClassicalOrthogonalPolynomials: grid

julia> T = chebyshevt(0..1);

julia> n = 5; x1 = grid(T, n-2) # 1stkind grid
3-element Vector{Float64}:
 0.9330127018922193
 0.5
 0.0669872981077807

I haven’t thought about 2nd kind points. It would be possible to add support for grid(T, n; kind=2)…

17

Thanks, in the meantime using

x1_ref = reverse(ChebyshevGrid{1}(n - 2))   # (n-2) first-kind nodes in [-1,1]
x2_ref = reverse(ChebyshevGrid{2}(n))       # n second-kind (Lobatto) nodes in [-1,1]
x1 = (x1_ref .+ 1) ./ 2            # map to [0,1]
x2 = (x2_ref .+ 1) ./ 2            # map to [0,1]

is good enough.

I have been trying to build more complex operators using ComplexOrthogonalPolynomials.jl, and I am getting a couple errors. Consider

T = chebyshevt(0 .. 1) # solution basis is T_n 
C = cp.ultraspherical(2, 0 .. 1) # RHS basis is C_n^(2)
s = axes(T, 1) 
p0 = s .^2 
u0 = 1 .- s .^ 3 
F_2 = 1 
mu_0 = 0.4 
C_mu = -14.6 

a = (3 * im * (F_2 * sqrt_p0 .* u0 .- 1) + 4 * k * p0 * (C_mu+ 9/mu_0)) .* T

1. Inexact Error

Now, trying to get the coefficients by doing (C\a)[1:15, 1:15] results in

julia> (C\a)[1:15, 1:15]
ERROR: InexactError: Float64(9.216666666666669 - 1.88671875im)
Stacktrace:
  [1] Real
    @ ./complex.jl:44 [inlined]
  [2] convert
    @ ./number.jl:7 [inlined]
  [3] setindex!
    @ ./array.jl:996 [inlined]
  [4] inbands_setindex!
    @ ~/.julia/packages/BandedMatrices/GLf5L/src/banded/BandedMatrix.jl:613 [inlined]
  [5] banded_setindex!
    @ ~/.julia/packages/BandedMatrices/GLf5L/src/banded/BandedMatrix.jl:623 [inlined]
  [6] setindex!
    @ ~/.julia/packages/BandedMatrices/GLf5L/src/banded/BandedMatrix.jl:634 [inlined]
  [7] _default_blasmul_loop!
    @ ~/.julia/packages/ArrayLayouts/cFEh2/src/muladd.jl:169 [inlined]
  [8] default_blasmul!(α::Float64, A::BandedMatrices.BandedMatrix{…}, B::BandedMatrices.BandedMatrix{…}, β::Bool, C::BandedMatrices.BandedMatrix{…})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/cFEh2/src/muladd.jl:191
  [9] materialize!(M::ArrayLayouts.MulAdd{…})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/cFEh2/src/muladd.jl:260
 [10] muladd!
    @ ~/.julia/packages/ArrayLayouts/cFEh2/src/muladd.jl:75 [inlined]
 [11] _mulbanded_copyto!
    @ ~/.julia/packages/LazyArrays/gHXdm/ext/LazyArraysBandedMatricesExt.jl:321 [inlined]
 [12] copyto!_layout(::BandedMatrices.BandedColumns{…}, srclay::LazyArrays.ApplyBandedLayout{…}, dest::BandedMatrices.BandedMatrix{…}, src::SubArray{…})
    @ LazyArraysBandedMatricesExt ~/.julia/packages/LazyArrays/gHXdm/ext/LazyArraysBandedMatricesExt.jl:328
 [13] copyto!_layout
    @ ~/.julia/packages/ArrayLayouts/cFEh2/src/ArrayLayouts.jl:273 [inlined]
 [14] copyto!
    @ ~/.julia/packages/ArrayLayouts/cFEh2/src/ArrayLayouts.jl:280 [inlined]
 [15] _BandedMatrix
    @ ~/.julia/packages/LazyArrays/gHXdm/ext/LazyArraysBandedMatricesExt.jl:303 [inlined]
 [16] BandedMatrix
    @ ~/.julia/packages/BandedMatrices/GLf5L/src/banded/BandedMatrix.jl:261 [inlined]
 [17] sub_materialize
    @ ~/.julia/packages/BandedMatrices/GLf5L/src/generic/indexing.jl:28 [inlined]
 [18] sub_materialize(L::LazyArrays.ApplyBandedLayout{…}, V::SubArray{…})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/cFEh2/src/ArrayLayouts.jl:136
 [19] sub_materialize(V::SubArray{Float64, 2, LazyArrays.ApplyArray{…}, Tuple{…}, false})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/cFEh2/src/ArrayLayouts.jl:137
 [20] layout_getindex
    @ ~/.julia/packages/ArrayLayouts/cFEh2/src/ArrayLayouts.jl:143 [inlined]
 [21] getindex(A::LazyArrays.ApplyArray{Float64, 2, typeof(*), Tuple{…}}, kr::UnitRange{Int64}, jr::UnitRange{Int64})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/cFEh2/src/ArrayLayouts.jl:158
 [22] top-level scope
    @ REPL[163]:1
Some type information was truncated. Use `show(err)` to see complete types.

2. MethodError

I tried to do this:

b = -(k/2) * (5 .+ p0).*T

but it results in:

julia> b = -(k/2) * (5 .+ p0).*T
ERROR: MethodError: no method matching _broadcasted_layout_broadcasted_mul(::Tuple{…}, ::QuasiArrays.BroadcastQuasiVector{…}, ::QuasiArrays.SubQuasiArray{…})
The function `_broadcasted_layout_broadcasted_mul` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  _broadcasted_layout_broadcasted_mul(::Tuple{ContinuumArrays.AbstractWeightLayout, QuasiArrays.PolynomialLayout}, ::Any, ::Any)
   @ ClassicalOrthogonalPolynomials ~/.julia/packages/ClassicalOrthogonalPolynomials/UM9Rj/src/clenshaw.jl:128

Stacktrace:
 [1] layout_broadcasted(::Tuple{…}, ::typeof(*), a::QuasiArrays.BroadcastQuasiVector{…}, P::QuasiArrays.SubQuasiArray{…})
   @ ClassicalOrthogonalPolynomials ~/.julia/packages/ClassicalOrthogonalPolynomials/UM9Rj/src/clenshaw.jl:125
 [2] broadcasted(::QuasiArrays.LazyQuasiArrayStyle{…}, ::Function, ::QuasiArrays.BroadcastQuasiVector{…}, ::QuasiArrays.SubQuasiArray{…})
   @ QuasiArrays ~/.julia/packages/QuasiArrays/YfIp9/src/lazyquasiarrays.jl:110
 [3] broadcasted(::typeof(*), ::QuasiArrays.BroadcastQuasiVector{…}, ::QuasiArrays.SubQuasiArray{…})
   @ Base.Broadcast ./broadcast.jl:1353
 [4] top-level scope
   @ REPL[166]:1
Some type information was truncated. Use `show(err)` to see complete types.

Note that doing:

b1 = (5 .+ p0).*T
b = (-k/2)*b1b1 = (5 .+ p0).*T
b = (-k/2)*b1

works fine.

18

Can you file an issue with a MWE?

Note that sqrt_p0 is not defined

19
20

These are now fixed.