# Solving the eigenvalue from a set of coupled second order differential equation numerically

I met a problem in solving a set of coupled differential equation, as shown below:

$$A_1\psi_1(z)+A_2\frac{d^2\psi_1(z)}{dz^2}+A_3\frac{d\psi_2(z)}{dz}=\lambda\psi_1(z)$$ $$A_4\psi_2(z)+A_5\frac{d^2\psi_2(z)}{dz^2}+A_3\frac{d\psi_1(z)}{dz}=\lambda\psi_2(z)$$ with the following 4 boundary condition:

$$\psi_1(0)=\psi_1(d)=0$$ $$\psi_2(0)=\psi_2(d)=0$$

where $$A_{i}$$ is a constant coefficient.

I've been stuck in this question for a long time. As the coupling term $$A_{3}$$ present in the problem. Without the boundary condition on the first order derivative, I was unable to determine the value of $$\lambda$$. Is there any numerical method to solve this type of eigenvalue problem?

This problem can be interpreted as a coupled convection-diffusion-reaction equation in two variables. You can use the Finite Element Method to solve it. Must be mentioned that you need to use some form of stabilisation for the convection term. You can use SUPG (Streamline Upwind Petrov Galerkin) stabilisation.

After discretisation, we can write it in the form the generalised eigenvalue problem

$$$$\mathbf{K} \mathbf{\Psi} = \lambda \mathbf{M} \mathbf{\Psi},$$$$ which can be solved for $$\lambda$$(s).

You can refer to this textbook for the details.

I would be happy to discuss further if you want.

• Thank you for you suggestion and kindly help. It seems this method is quite general but hard to implement. It will take me some time to go through it. I will learn this method for my later used. Aug 20, 2021 at 1:50

You can rewrite the problem in matrix form as

$$\begin{bmatrix} A_1 + A_2 \frac{d^2}{dz^2} &A_3\frac{d}{dz}\\ A_3\frac{d}{dz} &A_4 + A_5\frac{d^2}{dz^2} \end{bmatrix}\begin{Bmatrix} \psi_1\\ \psi_2\end{Bmatrix} = \lambda \begin{Bmatrix} \psi_1\\ \psi_2\end{Bmatrix}\, ,$$

that can be thought as

$$\mathcal{L} \Psi = \lambda \Psi\, ,$$

where $$\mathcal{L}$$ is your matrix operator and $$\Psi$$ is your vector variable.

Thus, you can think of methods for solving it. The following pop to my mind:

• Finite differences.

• Weighted residuals:

• Galerkin (such as FEM).
• Collocation.

Since you have a problem in 1D, I would try with finite differences or collocation first.

• I think you are right, but the given condition seems not sufficient to implement the finite difference method. Aug 20, 2021 at 1:55
• Sorry, I think the main problem of using finite difference method is the unknown $\lambda$ inside the equation which need to be determined. Aug 20, 2021 at 7:15
• Oh, you can use finite differences for eigenvalue problems (see scicomp.stackexchange.com/a/19679/9667). After discretization you end up with a matrix that represents the discrete version of your operator $\mathcal{L}$ and solve $Lu = \lambda u$. Aug 20, 2021 at 11:34
• Thank you for your reply. I try to write the operator into a discretized matrix form with a finite difference method. I encounter a problem when I try to use this method. As I want to study the change of $\lambda$ when I vary the boundary d, It seems that this method was unable to tell the difference when I change the boundary d. Aug 20, 2021 at 18:22
• By changing $d$ do you mean changing the size of the interval? If so, I would think that it is better to rewrite the equations in dimensionless form and use $s = z/d$ as variable. Aug 20, 2021 at 18:50

We can show how the eigensystem looks like for a given $$A_i, d$$, for example, in a case of $$A_1=A_4=-1, A_2=A_5=-1,A_3=1$$ and $$d=10$$ we have (in the pictures below $$\psi_1$$ - blue,$$\psi_2$$ - red, $$\lambda$$ is shown under pictures)

• Thank you for your help. May I ask which method did you use to calculate the result ? Your result look consistent to the boundary condition of the problem. It looks correct. My target is to solve the spectrum of E for different d given a constant for all $A_{i }$ Aug 20, 2021 at 1:53
• I have used NDEigensystem implemented in Mathematica. If you know $A_i$ then we can make a research about function $\lambda (d)$. Aug 20, 2021 at 2:02

The differential equations can also be written as the following system of first order differential equations

$$\frac{d}{dz}x(z) = M\,x(z), \tag{1}$$

with $$x(z)=\begin{bmatrix}\psi_1(z) & d/dz\,\psi_1(z) & \psi_2(z) & d/dz\,\psi_2(z)\end{bmatrix}^\top$$ and

$$M = \begin{bmatrix} 0 & 1 & 0 & 0 \\ \frac{\lambda-A_1}{A_2} & 0 & 0 & -\frac{A_3}{A_2} \\ 0 & 0 & 0 & 1 \\ 0 & -\frac{A_3}{A_5} & \frac{\lambda-A_4}{A_5} & 0 \\ \end{bmatrix}.$$

The solution to these differential equations can be obtained using

$$x(z) = e^{M\,z} x(0), \tag{2}$$

using the matrix exponential, defined as

$$e^X = \sum_{n=0}^\infty \frac{1}{n!}X^n = I + X + \frac{1}{2} X^2 + \cdots.$$

The first constraints can also be written as that $$x(0)$$ has to lie in the following span

$$\left\{\begin{pmatrix}0\\1\\0\\0\end{pmatrix},\begin{pmatrix}0\\0\\0\\1\end{pmatrix}\right\}.$$

The last constraints requires that $$x(d)$$ lies in that same span, which is equivalent that it lies orthogonal to the its null space. This is equivalent to

$$\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix} e^{M\,d} \begin{bmatrix} 0 & 0 \\ 1 & 0 \\ 0 & 0 \\ 0 & 1 \end{bmatrix} = \begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix}. \tag{3}$$

By taking some norm, such as the Frobenius norm, of the left hand side of $$(3)$$ the eigenvalue problem becomes a root finding problem. Though, I am not sure if this would be the easiest way of solving this, but I thought it might be worth mentioning (but would be too long for a comment).

For example when also using $$A_1 = A_2 = A_4 = A_5 = -1$$, $$A_3=1$$ and $$d = 10$$ and plotting the resulting values for the norm of the left hand side of $$(3)$$ yields the following plot

Which seems to roughly match the results from Alex Trounev. Namely, there seems to be a root near $$-0.65137$$, $$-0.35523$$, $$0.13826$$, $$0.82910$$, $$1.7173$$, $$2.8031$$, $$4.0862$$, $$5.5667$$, $$7.2445$$ and $$9.1196$$.

• But $\lambda$ is not known. Aug 21, 2021 at 23:39
• @nicoguaro I know, that it why one should look for the roots of the norm of the LHS of $(3)$ (as a function of $\lambda$). Aug 22, 2021 at 3:56
• But, turning to a root solving problem is commonly a bad idea. Indeed, that problem is sometimes solved by eigenvalues. Aug 22, 2021 at 4:36
• @nicoguaro Can't one combine this with discretizing the problem? Discretizing does allow one to solve it using eigenvalues, but discretizing always introduces some errors (especially for the bigger magnitude eigenvalues). For example using a uniform grid of 50 points in $z$ for the same example as in my answer yields $-0.65120$, $-0.35716$, $0.13156$, $0.81296$, $1.6842$, $2.7418$, $3.9813$, $5.3977$, $6.9851$ and $8.7370$. But one could use these obtained values as initial guesses for the root finding, which I think would have better convergence compared to increasing grid size. Aug 22, 2021 at 22:12
• Yes, that's a good approach. Aug 22, 2021 at 22:22