Strange eigenvalue

[rveltz]

Hi,

I am trying to get the spectrum of the following operator, it should be real negative but...

I am wondering if I am doing something wrong

using ApproxFun, LinearAlgebra
sol = Fun( X -> 0.0, Interval(-10.0, 10.0))
D = Derivative(sol.space, 2);
L1 = -I - 2D - D^2;
ApproxFun.eigs(Dirichlet(sol.space),L1-2I,100; tolerance = 1e-5)

gives

92-element Array{Complex{Float64},1}:
  -2.032800250240272 + 0.0im
  -2.004053699045951 - 0.4041140799101564im
  -2.004053699045951 + 0.40411407991015635im
  -2.366727744894575 + 0.0im
  -2.607891053151814 + 0.0im
  -2.791576558378409 + 0.0im
  -2.961489018882946 + 0.0im
 -3.1232270743384105 + 0.0im
   4.161842199302888 + 0.0im
    -4.5594165976354 + 0.0im
  -5.507851033473713 + 0.0im
 -12.381506648671865 + 0.0im
  13.298458875669132 + 0.0im
  -15.92506428457156 + 0.0im
 -28.008488707690706 + 0.0im
  -42.20133585945493 - 7.230886999195003im
 -42.201335859454936 + 7.230886999195003im
  -49.50513919128081 + 0.0im
  -92.45994637684655 + 0.0im
                     ⋮
  3067.6684441749544 + 709.7378828741976im
   4287.123993972118 + 0.0im
   4656.318955035678 + 0.0im
   6746.749676896203 + 0.0im
    7386.02040543896 + 0.0im
   10850.82308917364 + 0.0im
  11844.180136316863 + 0.0im
  17796.130834092044 + 0.0im
  19380.105512971324 + 0.0im
  30081.964227412103 + 0.0im
   32698.06214298154 + 0.0im
  53134.160137843304 + 0.0im
   57669.60736441608 + 0.0im
    99922.0088392208 + 0.0im
  108327.00882217327 + 0.0im
  205571.55994053124 + 0.0im
  222669.12343122496 + 0.0im
   482877.7424391819 + 0.0im
   522698.2906441737 + 0.0im

[dlfivefifty]

Probably your tolerance is too weak?

[rveltz]

Same issue with

ApproxFun.eigs(Dirichlet(sol.space),L1,100; tolerance = 1e-15)

[dlfivefifty]

You should be able to inspect the eigenvectors and check. Perhaps the pseudospectrum of your operator is very different from the spectrum, and therefore you have almost-eigenvectors.

[dlfivefifty]

Also I don't buy that it should be real-negative: -D^2 has positive eigenvalues and the other terms are relatively-compact perturbations

[rveltz]

D is the Laplacian

L1 is the operator for the Swift-Hohenberg equation, roughly speaking, the spectrum of L is -(1-z^2)^2

[dlfivefifty]

ah sorry very confusing to call D^2 D... would be clearer with Δ

[rveltz]

but then the problem remains?

[dlfivefifty]

like I said, maybe pseudospectra. check the eigenvector

[rveltz]

Seems like no:

λ, V = ApproxFun.eigs(Dirichlet(sol.space),L1,100; tolerance = 1e-12)
for ii in 1:length(λ)
	@show norm(L1*V[ind]-λ[ind]*V[ind])/norm(V[ind])
end

gives 1e-7

[dlfivefifty]

I would say that’s a yes, that all points are in the 10^-7 pseudospectra

note the tolerance only checks the tail of the coefficients of the eigenvectors. I’d suggest refiltering based on error in residual

[rveltz]

How do I do this?

[dlfivefifty]

Basically what you did with the for loop:

λ_keep = Vector{eltype(λ)}()
V_keep = Vector{eltype(V)}()
for (λ,V) in zip(λ,V)
    if norm(L1*V-λ*V)/norm(V) ≤ tol
        push!(λ_keep, λ)
        push!(V_keep, V)
    end
end
end

[rveltz]

I was really expecting this to be straightforward :(
I must be wrong somewhere,

sol = Fun( x -> 0.1*cos(x), Interval(-10pi, 10pi))
evsin = Fun( x -> sin(x), Interval(-10pi, 10pi))

Δ = Derivative(sol.space, 2)
L1 = -(I + 2Δ + Δ^2)

parms = (p = -0.1, b = 2.0, L1 = L1)
R_SH(sol, parms)[3] |> plot

λ, V = ApproxFun.eigs(Dirichlet(sol.space),L1,120; tolerance = 1e-9)
λ_keep = Vector{ComplexF64}()
	V_keep = Vector{eltype(V)}()
	for (_λ,_V) in zip(λ,V)
		# @show norm(L1*_V-_λ*_V)/norm(_V)
		if norm(L1*_V-_λ*_V)/norm(_V) ≤ 1e-6
			@info "one"
			push!(λ_keep, _λ)
			# push!(V_keep, _V)
		end
	end
	display(λ_keep)

but it does not even find the zero eigenvalue

julia> L1*evsin |> norm
1.538606260973298e-11

[dlfivefifty]

I think you underestimate how dodgy this stuff is, particularly since it solves a generalised eigenvalue problem. But please notice you don't have enough boundary conditions as you probably want to impose derivatives vanish as werll.

Since your operator is SPD you would be better off using a p-FEM which would result in a symmetric generalised eigenvalue problem. This should prohibit complex eigenvalues. This is doable in OrthogonalPolynomialsQuasi.jl, albeit a bit buggy. This seems to work:

julia> S = JacobiWeight(2,2) .* Jacobi(2,2);

julia> D = Derivative(axes(S,1));

julia> n = 100; # truncation size

julia> M = (S'S)[1:n,1:n]; # mass matrix

julia> L1 = -M + 2 * ((D*S)'*(D*S))[1:n,1:n] - ((D^2*S)'*(D^2*S))[1:n,1:n] # weak version of your operator
100×100 BandedMatrices.BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
 -21.5365       8.88178e-16    -4.58066         -1.38778e-17  …      ⋅               ⋅               ⋅               ⋅ 
   0.0       -150.608          -5.55112e-17    -10.7413              ⋅               ⋅               ⋅               ⋅ 
  -4.58066     -5.55112e-17  -486.016            3.55271e-15         ⋅               ⋅               ⋅               ⋅ 
   0.0        -10.7413         -1.16573e-15  -1124.42                ⋅               ⋅               ⋅               ⋅ 
  -0.042624    -6.59195e-17   -17.5611          -1.11022e-16         ⋅               ⋅               ⋅               ⋅ 
    ⋅          -0.0596737      -2.42861e-17    -24.7353       …      ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅             -0.065524         2.08167e-17         ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅              -0.0665093           ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅           …      ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅           …      ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅           …      ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅           …      ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅           …      ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
   ⋮                                                          ⋱                                                  
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅           …      ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅           …      ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅           …      ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅           …      ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅           …      ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                  ⋅               ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                -0.00993586       ⋅               ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅                 3.46945e-18    -0.00983826       ⋅               ⋅ 
    ⋅            ⋅               ⋅                ⋅              -748.327           3.46945e-18    -0.00974255       ⋅ 
    ⋅            ⋅               ⋅                ⋅           …    -1.04083e-17  -756.324           3.46945e-18    -0.00964869
    ⋅            ⋅               ⋅                ⋅                -1.46769e7       8.70241e-14  -764.321           8.67362e-18
    ⋅            ⋅               ⋅                ⋅                -1.04083e-17    -1.51348e7      -6.93889e-18  -772.318
    ⋅            ⋅               ⋅                ⋅              -764.321          -8.70484e-14    -1.5602e7        6.93889e-18
    ⋅            ⋅               ⋅                ⋅                 1.73472e-18  -772.318           0.0            -1.60788e7

julia> eigvals(Symmetric(Matrix(L1)), Symmetric(Matrix(M))) # L1, M are symmetric and banded but unfortunately we don't have generalised banded solver working yet
100-element Array{Float64,1}:
      -1.7741121549044663e12
      -1.642730926177833e12
      -2.3408751910288263e11
      -2.167633260265802e11
      -6.121712823185198e10
      -5.669226566942357e10
      -2.2572673800136753e10
      -2.090737214758445e10
      -1.0218397415260227e10
      -9.466467123271519e9
      -5.305514909569924e9
      -4.916382028793124e9
      -3.039930748381637e9
      -2.817866018646582e9
      -1.8766421283573024e9
      -1.7402194693226078e9
      -1.2285145021015291e9
      -1.139720297620768e9
      -8.434572630982723e8
      -7.829045797656573e8
      -6.025248316542239e8
      -5.596087002303199e8
      -4.452082478111935e8
      -4.137875411186183e8
      -3.387747859255062e8
      -3.151212956111086e8
      -2.6458950031324056e8
      -2.463480700653357e8
      -2.115782412917421e8
      -1.9720989895601687e8
      -1.7291742478386238e8
      -1.6138664176148584e8
      -1.4426768858655038e8
       ⋮
      -6.787129605960277e6
      -5.989272199764917e6
      -5.26390772803182e6
      -4.606506668032358e6
      -4.012685610664511e6
      -3.478207260474876e6
      -2.998980435634413e6
      -2.571060067957108e6
      -2.1906472028934928e6
      -1.8540889995271298e6
      -1.5578787305822165e6
      -1.2986557824164606e6
      -1.0732056550250317e6
 -878459.9620391804
 -711496.4307264198
 -569538.901990373
 -449957.3303705901
 -350267.7840424216
 -268132.4448165245
 -201359.60813843622
 -147903.68308765467
 -105865.19237627478
  -73490.77234662283
  -49173.172967071645
  -31451.257824161847
  -19010.00410753719
  -10680.50256655247
   -5439.957566979443
   -2411.684984001942
    -865.1341529737414
    -215.65507620069224
     -26.10247976293155

[rveltz]

you don't have enough boundary conditions

That's it I guess