# How to use DEigensystem with periodic boundary conditions on the derivative?

Since I found I can use DEigensystem to obtain eigenvalues and eigefunctions for differential operators, I am using it to verify solution to some of my HW's.

In this problem, from the textbook

I am not able to figure how the syntax to use. The eigenvalues are $n^2$ for $n=0,1,2,\dots$, and the eigenfunctions are $\{1,\cos(n x),\sin(n x)\}$

I do not know how to write the periodic boundary conditions for $y'(0)=y'(2\pi)$. This is what I tried (this below is still not complete)

ClearAll[y,x];
op={-y''[x],PeriodicBoundaryCondition[y[x],x==2 Pi],Function[x,x-2 Pi]]]};
eigf= NDEigenvalues[op,y[x],{x,0,2 Pi},5]


To add the derivatives part of the problem, I tried

op={-y''[x]+
PeriodicBoundaryCondition[y'[x],x==2 Pi,Function[x,x-2 Pi]],
PeriodicBoundaryCondition[y[x],x==2 Pi,Function[x,x-2 Pi]]};

eigf= NDEigenvalues[op,y[x],{x,0,2 Pi},5]


op={-y''[x],PeriodicBoundaryCondition[y'[x],x==2 Pi,Function[x,x-2 Pi]],
PeriodicBoundaryCondition[y[x],x==2 Pi,Function[x,x-2 Pi]]};

eigf= NDEigenvalues[op,y[x],{x,0,2 Pi},5]


I also noticed NDEigenvalues gives any result, while DEigenvalues does nothing at all.

Any suggestions how to use Mathematica to verify the above analytical results? Does PeriodicBoundaryCondition support derivatives? I can post the hand derivation if needed.

## Appendix

Hand derivation follows

Solution using transformation.

Let $\tau=x-\pi$, hence the above system becomes \begin{align*} \frac{d\phi^{2}}{d\tau^{2}}+\lambda\phi & =0\\ \phi\left( -\pi\right) & =\phi\left( \pi\right) \\ \frac{d\phi}{d\tau}\left( -\pi\right) & =\frac{d\phi}{d\tau}\left( \pi\right) \end{align*} The characteristic equation is $r^{2}+\lambda=0$ or $r=\pm\sqrt{-\lambda}$. Assuming $\lambda$ is real. There are three cases to consider.

Case $\lambda<0$

Let $s=\sqrt{-\lambda}>0$ \begin{align*} \phi\left( \tau\right) & =c_{1}\cosh\left( s\tau\right) +c_{2} \sinh\left( s\tau\right) \\ \phi^{\prime}\left( \tau\right) & =sc_{1}\sinh\left( s\tau\right) +sc_{2}\cosh\left( s\tau\right) \end{align*} Applying first B.C. gives \begin{align} \phi\left( -\pi\right) & =\phi\left( \pi\right) \nonumber\\ c_{1}\cosh\left( s\pi\right) -c_{2}\sinh\left( s\pi\right) & =c_{1} \cosh\left( s\pi\right) +c_{2}\sinh\left( s\pi\right) \nonumber\\ 2c_{2}\sinh\left( s\pi\right) & =0\nonumber\\ c_{2}\sinh\left( s\pi\right) & =0 \tag{1} \end{align} Applying second B.C. gives \begin{align} \phi^{\prime}\left( -\pi\right) & =\phi^{\prime}\left( \pi\right) \nonumber\\ -sc_{1}\sinh\left( s\pi\right) +sc_{2}\cosh\left( s\pi\right) & =sc_{1}\sinh\left( s\pi\right) +sc_{2}\cosh\left( s\pi\right) \nonumber\\ 2c_{1}\sinh\left( s\pi\right) & =0\nonumber\\ c_{1}\sinh\left( s\pi\right) & =0 \tag{2} \end{align} Since $\sinh\left( s\pi\right)$ is zero only for $s\pi=0$ and $s\pi$ is not zero because $s>0$. Then the only other option is that both $c_{1}=0$ and $c_{2}=0$ in order to satisfy equations (1)(2). Hence trivial solution. Hence $\lambda<0$ is not an eigenvalue.

Case $\lambda=0$

The space equation becomes $\frac{d\phi^{2}}{d\tau^{2}}=0$ with the solution $\phi\left( \tau\right) =A\tau+B$. Applying the first B.C. gives \begin{align*} \phi\left( -\pi\right) & =\phi\left( \pi\right) \\ -A\pi+B & =A\pi+B\\ 0 & =2A\pi \end{align*} Hence $A=0$. The solution becomes $\phi\left( \tau\right) =B$. And $\phi^{\prime}\left( \tau\right) =0$. The second B.C. just gives $0=0$. Therefore the solution is $$\phi\left( \tau\right) =C$$ Where $C$ is any constant. Hence $\lambda=0$ is an eigenvalue.

Case $\lambda>0$

\begin{align*} \phi\left( \tau\right) & =c_{1}\cos\left( \sqrt{\lambda}\tau\right) +c_{2}\sin\left( \sqrt{\lambda}\tau\right) \\ \phi^{\prime}\left( \tau\right) & =-c_{1}\sqrt{\lambda}\sin\left( \sqrt{\lambda}\tau\right) +c_{2}\sqrt{\lambda}\cos\left( \sqrt{\lambda} \tau\right) \end{align*} Applying first B.C. gives \begin{align} \phi\left( -\pi\right) & =\phi\left( \pi\right) \nonumber\\ c_{1}\cos\left( \sqrt{\lambda}\pi\right) -c_{2}\sin\left( \sqrt{\lambda}% \pi\right) & =c_{1}\cos\left( \sqrt{\lambda}\pi\right) +c_{2}\sin\left( \sqrt{\lambda}\pi\right) \nonumber\\ 2c_{2}\sin\left( \sqrt{\lambda}\pi\right) & =0\nonumber\\ c_{2}\sin\left( \sqrt{\lambda}\pi\right) & =0 \tag{3} \end{align} Applying second B.C. gives \begin{align} \phi^{\prime}\left( -\pi\right) & =\phi^{\prime}\left( \pi\right) \nonumber\\ c_{1}\sqrt{\lambda}\sin\left( \sqrt{\lambda}\pi\right) +c_{2}\sqrt{\lambda }\cos\left( \sqrt{\lambda}\pi\right) & =-c_{1}\sqrt{\lambda}\sin\left( \sqrt{\lambda}\pi\right) +c_{2}\sqrt{\lambda}\cos\left( \sqrt{\lambda} \pi\right) \nonumber\\ 2c_{1}\sqrt{\lambda}\sin\left( \sqrt{\lambda}\pi\right) & =0\nonumber\\ c_{1}\sin\left( \sqrt{\lambda}\pi\right) & =0 \tag{2} \end{align} Both (3) and (2) can be satisfied for non-zero $\sqrt{\lambda}\pi.$ The trivial solution is avoided. Therefore the eigenvalues are \begin{align*} \sin\left( \sqrt{\lambda}\pi\right) & =0\\ \sqrt{\lambda_{n}}\pi & =n\pi\qquad n=1,2,3,\cdots\\ \lambda_{n} & =n^{2}\qquad n=1,2,3,\cdots \end{align*}

Hence the corresponding eigenfunctions are

$$\left\{ \cos\left( \sqrt{\lambda_{n}}\tau\right) ,\sin\left( \sqrt {\lambda_{n}}\tau\right) \right\} =\left\{ \cos\left( n\tau\right) ,\sin\left( n\tau\right) \right\}$$

Transforming back to $x$ using $\tau=x-\pi$

$$\left\{ \cos\left( n\left( x-\pi\right) \right) ,\sin\left( n\left( x-\pi\right) \right) \right\} =\left\{ \cos\left( nx-n\pi\right) ,\sin\left( nx-n\pi\right) \right\}$$

But $\cos\left( x-\pi\right) =-\cos x$ and $\sin\left( x-\pi\right) =-\sin x$, hence the eigenfunctions are

$$\left\{ -\cos\left( nx\right) ,-\sin\left( nx\right) \right\}$$

The signs of negative on an eigenfunction (or eigenvector) do not affect it being such as this is just a multiplication by $-1$. Hence the above is the same as saying the eigenfunctions are

$$\left\{ \cos\left( nx\right) ,\sin\left( nx\right) \right\}$$

Summary

$\lambda=0$ eigenfunctions arbitrary constant

$\lambda>0$ eigenfunctions $\left\{ \cos\left( nx\right) ,\sin\left( nx\right) \right\}$ for $n=1,2,3\cdots$

• I believe that specifying both y and y' overspecifies the problem. However, specifying only y' fails with the 'Initial Condition Creation Failed" message. Apparently, NDEigensystem is not yet capable of handling this boundary condition. Aug 28, 2017 at 8:15

You may notice in this case, at least, you do find periodic boundary conditions on the derivatives with your original setup with PeriodicBoundaryCondition handed an argument only in terms of the function and not its derivatives.

This is perhaps not so surprising given the fact that your differential operator ${\cal O}\equiv -\partial_x^2$ commutes with differentiation (not always the case, consider e.g. a distinct differential operator ${\cal O'}\equiv-\partial_x^2 + x + 5$) which means that $y'(x)$ will be eigenfunctions of your system with exactly the same eigenvalues as $y(x)$. To see this differentiate your differential equation on both sides, and pass the differentiation through the operator. $[{\cal O,\partial_x}]=0 \,\, \& \,\, {\cal O} y = \lambda y \implies {\cal O}(\partial_x y) = \lambda (\partial_x y)$.

So $[{\cal O},\partial_x]=0$ means that any solution to ${\cal O} u = \lambda u \,\, \& \,\, u(x_1)=u(x_2)$ is simultaneously a solution to both $u=y$ and $u=\partial_x y$ with associated boundary conditions.

Here we'll use NDEigensystem so we can verify the properties of the eigenfunctions do indeed satisfy differential periodic boundary conditions associated with the eigenvalue solution you calculated explicitly with NDEigenvalues.

ClearAll[y, x];
op = {-y''[x],
PeriodicBoundaryCondition[y[x], x == 2 Pi,
Function[x, x - 2 Pi]]};
{ev, ef} = NDEigensystem[op, y[x], {x, 0, 2 Pi}, 5];


### Direct check of derivatives

To directly address whether or not your solution's derivatives match at the boundary we can of course just check:

Table[Abs[(D[ef[[i]], x] /. x -> 0) - (D[ef[[i]], x] /.
x -> 2 \[Pi])], {i, 1, 5}]


yielding:

{4.94753*10^-15, 0.00114434, 0.000181246, 0.0172912, 0.00561824}

OK these are smallish numbers which is good, but to make sense of whether they are close enough to zero to count as satisfying boundary conditions, we need to make some sense of the actual precision of the solution.

### Understanding solution by relating to your analytics

My goal here is to show that we are indeed landing on a numerical approximation of periodic boundary conditions imposed for both $y(x)$ and $y'(x)$. I'll do this by relating to the actual solutions.

Plot[ef /. x -> xx // Evaluate, {xx, 0, 2 \[Pi]},
PlotPoints -> 50]


Of course since it is always possible to superimpose eigenfunctions with the same eigenvalues that's exactly what happens. We'll want to decompose into your familiar basis.

We can decompose either with a fit:

fitSol = Table[
FindFit[Table[{y, (ef[[i]] /. x -> xx )}, {xx, 0, 2 \[Pi],
2 \[Pi]/100}],
a Cos[xx Floor[i/2] ] + b Sin[xx Floor[i/2]], {a, b}, xx], {i, 1, 5}]


{{a -> -0.398942, b -> 0.}, {a -> 0.557245, b -> -0.0882589}, {a -> 0.0882589, b -> 0.557245}, {a -> -0.536595, b -> 0.17435}, {a -> -0.17435, b -> -0.536593}}

Or, of course, we could just use the orthogonality of the solutions:

orthoSol =
Table[Module[{csnorm =
Integrate[Cos[xx Floor[i/2] ]^2, {xx, 0, 2 \[Pi]}],
snorm = If[i > 1,
Integrate[Sin[xx Floor[i/2] ]^2, {xx, 0, 2 \[Pi]}],
1]}, {a ->
NIntegrate[
ef[[i]] Cos[xx Floor[i/2] ] /. x -> y, {xx, 0, 2 \[Pi]}]/
csnorm, b ->
NIntegrate[
ef[[i]] Sin[xx Floor[i/2] ] /. x -> y, {xx, 0, 2 \[Pi]}]/
snorm
}], {i, 1, 5}]


{{a -> -0.398942, b -> 0.}, {a -> 0.557264, b -> -0.0882552}, {a -> 0.0882619, b -> 0.557221}, {a -> -0.536576, b -> 0.174281}, {a -> -0.174344, b -> -0.536381}}

In any case, we can see that we do indeed have a nice superposition per eigenvalue, by using these to map back to the familiar basis and look at the difference functions:

Plot[ Table[(If[Abs[a] > Abs[b],
ef[[i]]/a - b/a Sin[ x Floor[ i/2]] - Cos[ x Floor[ i/2]],
funs[[i]]/b - a/b Cos[x Floor[i/2]] - Sin[ x Floor[ i/2]]] /.
orthoSol[[i]]), {i, 1, 5}] /. x -> xx // Evaluate, {xx, 0,
2 \[Pi]},
PlotLegends -> {"vs 1", "vs cos(x)", "vs sin(x)", "vs cos(2x)",
"vs sin(2x)"}, PlotRange -> All]


Can compare magnitude variance from exact solution, with above explicit derivatives case by case to see that we're of the right order of accuracy for our derivatives to be considered matching at the boundary. Not so bad for machine precision via FEM.

As you already know from working it out analytically, these solutions do indeed satisfy both periodic boundary conditions ($y(2\pi)=y(0)$ and $y'(2\pi)=y'(0)$).

### Numeric differential eigensystems only?

Re: NDEigen* vs DEigen* and PeriodicBoundaryConditions, it seems like it's only been promoted for numerical solution. (c.f. blurb here). You may learn more browsing user21's answers on this site.