NDEigensystem: 1D problem with discontinuous coefficients

Question

I am trying to use NDEigensystem to solve the 1D problem

-cs[x]^2 vx''[x] = w^2 vx[x]

with vx[x] and w the eigenfunction and eigenvalue. The coefficient cs[x] is discontinuous at x = -xp, +xp and vx[x] must vanish at the boundaries x = -L, +L.

This problem has an analytical solution, but I want to solve it numerically because it is an extremely simplified version of a similar problem with three coupled second order ODEs.

csp = 1;
csc = 13.83;
xp = 1;
L = 2;
cs[x_] := Piecewise[{{csp, Abs[x] <= xp}, {csc, Abs[x] > xp}}]
(*
Δ=0.1;
cs[x_]:=Piecewise[{{(Tanh[-(x+xp)/Δ]+1)/2*(csc-csp)+csp,\
x<0},{(Tanh[(x-xp)/Δ]+1)/2*(csc-csp)+csp,x>=0}}]
*)
Plot[cs[x], {x, -L, L}, 
 AxesLabel -> {"x", "\!\(\*SubscriptBox[\(c\), \(s\)]\)"}]
ℒ := -cs[x]^2 vx''[x];
ℬ = DirichletCondition[vx[x] == 0, True];
Neigen = 2;
{vals, funs} = 
  NDEigensystem[{ℒ, ℬ}, {vx[
     x]}, {x, -L, L}, Neigen];
i = 1;
Print[Sqrt[vals[[i]]]]
Plot[funs[[i]], {x, -L, L}, PlotRange -> All, 
 PlotLabel -> Sqrt[vals[[i]]]]

Both vx[x] and its first-order derivative should be continuous at the discontinuity points (x = -1, +1), but vx'[x] is discontinuous there:

It is possible to achieve continuity of vx'[x] by using the smoothed version of cs[x] which is commented out in the code above:

But of course this is not the problem I am interested in solving.

Answer 1

Yes, it's possible, we just need to use the documented but seldom-used and (probably) confusing last syntax of NDEigensystem.

As discussed in

Position of discontinuous coefficient influences the solution of PDE

, for time dependent PDE solved with FiniteElement method, if the discontinuous coefficient is in front of $\frac{\partial u}{\partial t}$ term rather than the diffusion term, the spatial derivative of solution will be continuous, as you expect. "But there's no such term when solving eigenvalue problem!" Well, actually, there is. As mentioned in the document of NDEigensystem:

NDEigensystem[eqns, {u, …}, t, {x, y, …} ∈ Ω, n]

gives the eigenvalues and eigenfunctions in the spatial variables {x, y, …} for solutions u, … of the coupled time-dependent differential equations eqns.

What does this mean? Some explanation can be found in Details and Options section:

For a system of first-order time-dependent equations, the time derivatives D[u[t, x, y, …], t], D[v[t, x, y, …], t], … are effectively replaced with λ u[x, y, …], λ v[x, y, …], ….

OK, you may still feel puzzled, then some further explanation can be found in the comments of user21 here:

… the idea is based on the fact that the FEM discretization of D[u, t] is the same as u (without temporal derivative). In other words the reaction term is the same as the transient term D[u, t] or any higher order transient term. … note that I am not solving $\frac{\partial u}{\partial t} = ∇ a$. I am just using that to construct the separated damping and stiffness system matrices. …

So, to summarize, with the last syntax of NDEigensystem, we're solving an eigenvalue problem corresponding to the time-dependent differential equation(s), internally the $\frac{\partial u}{\partial t}$ term is effectively replaced by $\lambda u$.

Knowing all these, the problem is easy to solve:

csp = 1;
csc = 13.83;
xp = 1;
L = 2;
cs[x_] := Piecewise[{{csp, Abs[x] <= xp}, {csc, Abs[x] > xp}}]
eq = D[vx[t, x], x, x] +(* See here! *)1/cs[x]^2 D[vx[t, x], t] == 0;
ℬ = DirichletCondition[vx[t, x] == 0, True];
Neigen = 2;
{vals, funs} = 
  NDEigensystem[{eq, ℬ}, {vx[t, x]}, t, {x, -L, L}, Neigen, 
   Method -> {"PDEDiscretization" -> {"FiniteElement", 
       "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}}];
i = 1;
Plot[{funs[[i]], D[funs[[i]], x]} // Evaluate, {x, -L, L}, PlotRange -> All, 
 PlotLabel -> Sqrt[vals[[i]]]]