# Eigen value and eigen functions of linear operator
**URL:** https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859
**Category:** Numerics
**Tags:** question, approxfun, eigenvalue-problem, chebychev
**Created:** [November 14, 2025, 12:01am UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859 "2025-11-14T00:01:23Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@Liris](https://discourse.julialang.org/u/Liris)
#### Post date: [November 14, 2025, 12:01am UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/1 "2025-11-14T00:01:23Z")
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:
- 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!
---
### Author:  [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)
#### Post date: [November 14, 2025, 8:14am UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/2 "2025-11-14T08:14:57Z")
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
---
### Author:  [@Liris](https://discourse.julialang.org/u/Liris)
#### Post date: [November 14, 2025, 3:28pm UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/4 "2025-11-14T15:28:08Z")
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) ?
---
### Author:  [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)
#### Post date: [November 14, 2025, 3:57pm UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/5 "2025-11-14T15:57:46Z")
> 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](https://juliaregistries.github.io/General/packages/redirect_to_repo/ClassicalOrthogonalPolynomials) over [ApproxFun.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/ApproxFun)?
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.
---
### Author:  [@Liris](https://discourse.julialang.org/u/Liris)
#### Post date: [November 14, 2025, 4:44pm UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/7 "2025-11-14T16:44:36Z")
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
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [November 14, 2025, 6:49pm UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/8 "2025-11-14T18:49:25Z")
> [@Liris](#):
>
> I am new to these kinds of problems, but I understand that I have basically two dominant ways of approaching this:
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.
---
### Author:  [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)
#### Post date: [November 15, 2025, 12:46pm UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/9 "2025-11-15T12:46:16Z")
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:
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.
---
### Author:  [@Liris](https://discourse.julialang.org/u/Liris)
#### Post date: [November 16, 2025, 4:04pm UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/10 "2025-11-16T16:04:09Z")
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.
---
### Author:  [@Liris](https://discourse.julialang.org/u/Liris)
#### Post date: [November 16, 2025, 4:21pm UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/11 "2025-11-16T16:21:32Z")
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?
- How to reference [ClassicalOrthogonalPolynomials.jl](https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl/tree/main) in a publication? Simply link to the GitHub repo?
---
### Author:  [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)
#### Post date: [November 17, 2025, 10:19am UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/12 "2025-11-17T10:19:47Z")
> 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 😅
Some equations naturally have only 3 bcs. For example, I think if one of the operators is degenerate:
```julia
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
> **[CornerTaus.pdf](https://danfortunato.com/papers/CornerTaus.pdf)**
>
> 654.01 KB
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](https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl/tree/main) in a publication? Simply link to the GitHub repo?
You can just cite my Acta Numerica paper, I’ve added a CITATION.bib:
>
---
### Author:  [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)
#### Post date: [November 17, 2025, 10:23am UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/13 "2025-11-17T10:23:00Z")
I believe
[https://academic.oup.com/mnras/article/455/4/4274/1267403](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:
> **[A simple proof that the hp-FEM does not suffer from the pollution effect for...](https://link.springer.com/article/10.1007/s10444-023-10025-3)**
>
> In d dimensions, accurately approximating an arbitrary function oscillating with frequency ≲ k $\\lesssim k$ requires ∼ k d $\\sim k^{d}$ degrees of freedom. A numerical method for solving the Helmholtz equation (with wavenumber k) suffers from the...
(I believe a consequence of this is that spectral methods also avoid pollution…)
---
### Author:  [@Liris](https://discourse.julialang.org/u/Liris)
#### Post date: [November 18, 2025, 5:21pm UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/14 "2025-11-18T17:21:50Z")
Thanks a lot for all these additional references!
---
### Author:  [@Liris](https://discourse.julialang.org/u/Liris)
#### Post date: [December 10, 2025, 11:43am UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/15 "2025-12-10T11:43:19Z")
I have played around with [ClassicalOrthogonalPolynomials.jl](https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl/tree/main) 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]:
```julia-auto
T = chebyshevt(0 .. 1) # solution basis is T_n
```
when building the grid;
```julia-auto
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 ?
---
### Author:  [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)
#### Post date: [December 10, 2025, 1:31pm UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/16 "2025-12-10T13:31:04Z")
For 1st kind points there’s a `grid` routine:
```julia
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)`…
---
### Author:  [@Liris](https://discourse.julialang.org/u/Liris)
#### Post date: [December 11, 2025, 3:21pm UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/17 "2025-12-11T15:21:24Z")
Thanks, in the meantime using
```julia-auto
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
```julia-auto
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-auto
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:
```julia-auto
b = -(k/2) * (5 .+ p0).*T
```
but it results in:
```julia-auto
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:
```julia-auto
b1 = (5 .+ p0).*T
b = (-k/2)*b1b1 = (5 .+ p0).*T
b = (-k/2)*b1
```
works fine.
---
### Author:  [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)
#### Post date: [December 12, 2025, 6:11am UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/18 "2025-12-12T06:11:29Z")
Can you file an issue with a MWE?
Note that `sqrt_p0` is not defined
---
### Author:  [@Liris](https://discourse.julialang.org/u/Liris)
#### Post date: [December 12, 2025, 11:28am UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/19 "2025-12-12T11:28:48Z")
>
>
> Doing
>
> \`\`\`julia
> using ClassicalOrthogonalPolynomials
>
> T = chebyshevt(0 .. 1) # …solution basis is T\_n
> s = axes(T, 1)
> f = s .^ 2
>
> op = 2 \* f .\* T
> \`\`\`
>
> results in
>
> \`\`\`julia
>
> 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\[256\]:1
> Some type information was truncated. Use \`show(err)\` to see complete types.
> \`\`\`
>
> However, doing this:
>
> \`\`\`julia
> op\_tp = f .\* T
> op = 2\*op\_tp
> \`\`\`
> or this
>
> \`\`\`julia
> op = (2 \* f .- 1) .\* T
> \`\`\`
>
> or this
>
> \`\`\`julia
> f = s
> op = 2 \* f .\* T
> \`\`\`
>
> works.
>
>
> Doing
>
> \`\`\`julia
> using ClassicalOrthogonalPolynomials
>
> T = chebyshevt(0 .. 1) # s…olution basis is T\_n
> C = ultraspherical(1, 0 .. 1) # RHS basis is C\_n^(2)
> s = axes(T, 1)
>
> op = (s .+ im) .\* T
> C \\ op
> \`\`\`
>
> results in
>
> \`\`\`julia
> julia\> C \\ op
> (ℵ₀×ℵ₀ BandedMatrix{Float64} with bandwidths (0, 2) with data vcat(1×ℵ₀ FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}} with indices Base.OneTo(1)×OneToInf(), 1×ℵ₀ FillArrays.Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}} with indices Base.OneTo(1)×OneToInf(), hcat(1×1 Ones{Float64}, 1×ℵ₀ FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}} with indices Base.OneTo(1)×OneToInf()) with indices Base.OneTo(1)×OneToInf()) with indices Base.OneTo(3)×OneToInf() with indices OneToInf()×OneToInf()) \* (ℵ₀×ℵ₀ view(::RecurrenceRelationshipArrays.Clenshaw{Float64, Vector{ComplexF64}, LazyArrays.BroadcastVector{Float64, typeof(\*), Tuple{LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfUnitRange{Int64}}}, Float64}}, LazyArrays.BroadcastVector{Float64, typeof(\*), Tuple{LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfUnitRange{Int64}}}, Float64}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfUnitRange{Int64}}}, LazyBandedMatrices.Tridiagonal{Float64, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, FillArrays.Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, FillArrays.Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, :, :) with eltype Float64 with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf():
> Error showing value of type LazyArrays.ApplyArray{Float64, 2, typeof(\*), Tuple{BandedMatrices.BandedMatrix{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(vcat), Tuple{FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}}, FillArrays.Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}}, LazyArrays.ApplyArray{Float64, 2, typeof(hcat), Tuple{FillArrays.Ones{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}}}}}}, InfiniteArrays.OneToInf{Int64}}, SubArray{Float64, 2, RecurrenceRelationshipArrays.Clenshaw{Float64, Vector{ComplexF64}, LazyArrays.BroadcastVector{Float64, typeof(\*), Tuple{LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfUnitRange{Int64}}}, Float64}}, LazyArrays.BroadcastVector{Float64, typeof(\*), Tuple{LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfUnitRange{Int64}}}, Float64}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfUnitRange{Int64}}}, LazyBandedMatrices.Tridiagonal{Float64, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, FillArrays.Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, FillArrays.Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}}}:
>
> SYSTEM (REPL): showing an error caused an error
> ERROR: 1-element ExceptionStack:
> InexactError: Float64(0.5 + 1.0im)
> Stacktrace:
> \[1\] Real
> @ ./complex.jl:44 \[inlined\]
> \[2\] convert
> @ ./number.jl:7 \[inlined\]
> \[3\] setindex!
> @ ./array.jl:985 \[inlined\]
> \[4\] macro expansion
> @ ./multidimensional.jl:981 \[inlined\]
> \[5\] macro expansion
> @ ./cartesian.jl:64 \[inlined\]
> \[6\] \_unsafe\_getindex!
> @ ./multidimensional.jl:979 \[inlined\]
> \[7\] \_unsafe\_getindex(::IndexCartesian, ::SubArray{…}, ::UnitRange{…}, ::Int64)
> @ Base ./multidimensional.jl:970
> \[8\] \_getindex
> @ ./multidimensional.jl:956 \[inlined\]
> \[9\] getindex
> @ ./abstractarray.jl:1342 \[inlined\]
> \[10\] map (repeats 2 times)
> @ ./tuple.jl:409 \[inlined\]
> \[11\] \_mul\_getindex
> @ ~/.julia/packages/LazyArrays/gHXdm/src/linalg/mul.jl:219 \[inlined\]
> \[12\] getindex
> @ ~/.julia/packages/LazyArrays/gHXdm/src/linalg/mul.jl:185 \[inlined\]
> \[13\] isassigned(::LazyArrays.ApplyArray{Float64, 2, typeof(\*), Tuple{…}}, ::Int64, ::Int64)
> @ Base ./multidimensional.jl:1653
> \[14\] alignment(io::IOContext{…}, X::AbstractVecOrMat, rows::Vector{…}, cols::Vector{…}, cols\_if\_complete::Int64, cols\_otherwise::Int64, sep::Int64, ncols::Infinities.InfiniteCardinal{…})
> @ Base ./arrayshow.jl:68
> \[15\] \_print\_matrix(io::IOContext{…}, X::AbstractVecOrMat, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64, rowsA::InfiniteArrays.InfUnitRange{…}, colsA::InfiniteArrays.InfUnitRange{…})
> @ Base ./arrayshow.jl:207
> \[16\] print\_matrix(io::IOContext{…}, X::LazyArrays.ApplyArray{…}, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64)
> @ Base ./arrayshow.jl:171
> \[17\] print\_matrix
> @ ./arrayshow.jl:171 \[inlined\]
> \[18\] print\_array
> @ ./arrayshow.jl:358 \[inlined\]
> \[19\] show(io::IOContext{…}, ::MIME{…}, X::LazyArrays.ApplyArray{…})
> @ Base ./arrayshow.jl:399
> \[20\] show\_repl(io::IO, mime::MIME{Symbol("text/plain")}, x::Any)
> @ REPL ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/REPL.jl:548
> \[21\] show\_limited(io::IO, mime::MIME, x::Any)
> @ REPL ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/REPL.jl:517
> \[22\] (::REPL.var"#display##0#display##1"{REPL.REPLDisplay{…}, MIME{…}, Base.RefValue{…}})(io::Any)
> @ REPL ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/REPL.jl:540
> \[23\] with\_repl\_linfo(f::Any, repl::REPL.LineEditREPL)
> @ REPL ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/REPL.jl:833
> \[24\] display(d::REPL.REPLDisplay, mime::MIME{Symbol("text/plain")}, x::Any)
> @ REPL ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/REPL.jl:526
> \[25\] display
> @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/REPL.jl:546 \[inlined\]
> \[26\] display(x::Any)
> @ Base.Multimedia ./multimedia.jl:340
> \[27\] (::REPL.var"#print\_response##2#print\_response##3")()
> @ REPL ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/REPL.jl:590
> \[28\] repl\_backend\_loop(backend::REPL.REPLBackend, get\_module::Function)
> @ REPL ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/REPL.jl:445
> \[29\] start\_repl\_backend(backend::REPL.REPLBackend, consumer::Any; get\_module::Function)
> @ REPL ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/REPL.jl:427
> \[30\] start\_repl\_backend
> @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/REPL.jl:424 \[inlined\]
> \[31\] run\_repl(repl::REPL.AbstractREPL, consumer::Any; backend\_on\_current\_task::Bool, backend::Any)
> @ REPL ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/REPL.jl:653
> \[32\] run\_repl(repl::REPL.AbstractREPL, consumer::Any)
> @ REPL ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/REPL.jl:639
> \[33\] run\_std\_repl(REPL::Module, quiet::Bool, banner::Symbol, history\_file::Bool)
> @ Base ./client.jl:478
> \[34\] run\_main\_repl(interactive::Bool, quiet::Bool, banner::Symbol, history\_file::Bool)
> @ Base ./client.jl:499
> \[35\] repl\_main
> @ ./client.jl:586 \[inlined\]
> \[36\] \_start()
> @ Base ./client.jl:561
> Stacktrace:
> \[1\] print\_response(errio::IO, response::Any, backend::Union{…}, show\_value::Bool, have\_color::Bool, specialdisplay::Union{…})
> @ REPL ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/REPL.jl:599
> \[2\] (::REPL.var"#print\_response##0#print\_response##1"{REPL.LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
> @ REPL ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/REPL.jl:558
> \[3\] with\_repl\_linfo(f::Any, repl::REPL.LineEditREPL)
> @ REPL ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/REPL.jl:833
> \[4\] print\_response(repl::REPL.AbstractREPL, response::Any, show\_value::Bool, have\_color::Bool)
> @ REPL ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/REPL.jl:556
> \[5\] (::REPL.var"#do\_respond#73"{…})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
> @ REPL ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/REPL.jl:1192
> \[6\] (::REPL.var"#setup\_interface##32#setup\_interface##33"{…})(::REPL.LineEdit.MIState, ::Any, ::Vararg{…})
> @ REPL ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/REPL.jl:1599
> \[7\] (::REPL.LineEdit.var"#match\_input##0#match\_input##1"{…})(s::Any, p::Any)
> @ REPL.LineEdit ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/LineEdit.jl:1810
> \[8\] macro expansion
> @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/LineEdit.jl:2960 \[inlined\]
> \[9\] macro expansion
> @ ./lock.jl:376 \[inlined\]
> \[10\] (::REPL.LineEdit.var"#prompt!##2#prompt!##3"{…})()
> @ REPL.LineEdit ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/REPL/src/LineEdit.jl:2949
> Some type information was truncated. Use \`show(err)\` to see complete types.
> \`\`\`
---
### Author:  [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)
#### Post date: [December 12, 2025, 2:10pm UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/20 "2025-12-12T14:10:05Z")
These are now fixed.
---
### Author:  [@Liris](https://discourse.julialang.org/u/Liris)
#### Post date: [January 27, 2026, 3:59pm UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/21 "2026-01-27T15:59:03Z")
@dlfivefifty
I have progressed on this issue but the real system of equations is tricky so progress is a bit slow. I have a questions about the “rectangular” part of those methods.
Consider a system of equations with two variables such that:
u'' + v = \lambda v (1)
v' + u = \lambda u (2)
Three BCs are required for this system (the first equation is 2nd order, and the other is 1st order). To make the final matrix rectangular, you need:
- in the ultraspherical method, to truncate the operator corresponding to (1) at n-2, and the operator corresponding to (2) at n-3, then resulting in a (n-3) x n matrix which can be made square by appending the BCs.
- in the collocation method, it is the number of collocation points that you reduce. For example, you use `x₁ = reverse(ChebyshevGrid{1}(n-2)) # 1st kind points, sorted` [here](https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl/blob/main/examples/collocation.jl). In the case of the system above, then two grids
```julia-auto
x₁ = reverse(ChebyshevGrid{1}(n-2)) # 1st kind points, sorted
x₁_bis = reverse(ChebyshevGrid{1}(n-1)) # 1st kind points, sorted
```
can be used for the (1) and (2) equations to end up with a square matrix. However, the two equations are now discretised on different grids. I assume that there should be an extra step compared to the ultraspherical method.
However, I am not sure how to deal with this. Could you help?
---
### Author:  [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)
#### Post date: [January 28, 2026, 9:20am UTC](https://discourse.julialang.org/t/eigen-value-and-eigen-functions-of-linear-operator/133859/22 "2026-01-28T09:20:20Z")
I don’t really understand the issue: you can just use different grids for each equation.