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.