5

I am trying to solve the eigenvalue problem of a 1st-order ODE system using NDEigenvalue. It should be finite difference method for ODE. And I want to tune the the DifferenceOrder and also the number of mesh points.
The relevant option I found in NDEigenvalue seems to be PDEDiscretization that probably descretizes the spatial part and in my case there's only a spatial part. After some searching, I still cannot find the correct way to set the DifferenceOrder and the number of mesh points.

The following is my stupid try (the particular ODE maybe not so important for this question?)

eps0 = 1*^-10; points = 500;
f[R_, r_] := Exp[r/R] - 1(*Tan[π/2r/R]*); 
m0[R_, r_] := 1.0/R f[R, r]; A[l_, B_, r_] := l/r + B/2 r;
Fop1[y_, l_, pm_] := I (-D[y, r] + pm A[l, B, r] y);
variables = {α, β, γ, δ};
lhs = {Fop1[δ[r], l + 1, -1] + kz γ[r] + 
    m0[R, r] α[r], 
   Fop1[γ[r], l, 1] - kz δ[r] + m0[R, r] β[r], 
   Fop1[β[r], l + 1, -1] + kz α[r] - 
    m0[R, r] γ[r], 
   Fop1[α[r], l, 1] - kz β[r] - m0[R, r] δ[r]};
bc = DirichletCondition[
  Table[component@r == 0, {component, variables}], True]

l = -2; R = 2; B = 20; kz = 2; Rcutoff = 
 1.0 Min[6 R, 16/Sqrt[B]]; Nless = 50;
Re@NDEigenvalues[{lhs, bc}, variables, {r, eps0, Rcutoff}, Nless, 
  Method -> {"PDEDiscretization" -> {"MethodOfLines",  
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "MaxPoints" -> points, "MinPoints" -> points, 
        "DifferenceOrder" -> 2}}}]
xzczd
  • 65,995
  • 9
  • 163
  • 468
xiaohuamao
  • 4,728
  • 18
  • 36

1 Answers1

3

Several issues here:

  1. "MethodOfLines" is a method for solving well-posed initial-boundary value problem so we can't use it for eigenvalue problem. Despite actually being used (see the comments below for more details), "MethodOfLines" isn't an option for NDEigenvalue.

  2. Though NDSolve`FiniteDifferenceDerivative facilitates self-implementation of finite difference method (FDM), FDM hasn't been implemented in NDSolve and its friends so far, that's the reason why I wrote pdetoode/pdetoae.

  3. At the moment, NDEigenvalue together with NDEigensystem are based on "FiniteElement" method only. The option amounts to difference order of FDM is "MeshOrder", which unfortunately cannot be higher than 2 up to now.

So, in short, it's not possible to use FDM in NDEigenvalue etc., at least now.

xzczd
  • 65,995
  • 9
  • 163
  • 468
  • 2
    The first point is not quite correct. In fact we use "MethodOfLines" for every PDE when using NDEigensystem. This is shown here. – user21 Apr 13 '18 at 05:44
  • It should be possible to adapt the code here to use the FDM approach. – user21 Apr 13 '18 at 05:45
  • @user21 Er… I think you've mis-pasted the link in the first comment? (It's the same as the one in the second comment. ) – xzczd Apr 13 '18 at 06:34
  • Yes, they are the same link. The point is that every PDE given to NDEigensystem is transformed into a time dependent PDE. – user21 Apr 13 '18 at 06:46
  • Thanks! Let's see if there would be any other answer. So you mean NDEigenvalues uses difference order 2 in FEM sense. Hmm,,, It gives the same pattern as DifferenceOrder 4 in the FDM based on your code (the DifferenceOrder dependent problem I mentioned previously). Probably I will post a separate question simply addressing this strange problem itself. – xiaohuamao Apr 13 '18 at 09:11