I am looking at the following system of ODEs:

\begin{array}{r}{\left[c_{2}(k)-\partial_{\tau}^{2}\right] \varphi_{2}\left(\tau \right)=f_{21}(\tau) \varphi_{1}\left(\tau \right)} \\ {\left[c_{1}(k)-\partial_{\tau}^{2}\right] \varphi_{1}\left(\tau \right)=f_{11}(\tau) \varphi_{1}\left(\tau \right)+f_{12}(\tau) \varphi_{2}\left(\tau \right)} \end{array}

where $c_1(k)$ and $c_2(k)$ are known scalar functions of an unknown number $k$ and $f_{21}, f_{12}, f_{11}$ are known functions. The unknowns are $k, \varphi_{1,2}$ and I'd like to find $k$ numerically.

When there's only one ODE, $c(k)$ can be obtained by discretising time and solving a matrix eigenvalue problem:

\begin{array}{r}{\left[c(k)-\partial_{\tau}^{2}\right] \varphi\left(\tau \right)=f(\tau) \varphi\left(\tau \right)} \end{array}

I am curious if I can recast the coupled problem as a (generalized)eigenvalue problem as well? Otherwise, another possibility is to try a bisection-like method to find the correct $k$ so that the matrix is singular, but this seems clunky.

Any comments/references appreciated!

  • $\begingroup$ Where did you get this problem from? $\endgroup$
    – nicoguaro
    Sep 6, 2019 at 23:37
  • 2
    $\begingroup$ It comes up in the context of studying growth rate of perturbations in certain dynamical systems. But, I have a feeling it appears more generally! $\endgroup$
    – KartMan
    Sep 7, 2019 at 3:14
  • $\begingroup$ I think that depending on the form of $c(k)$ it can be written as a polynomial or nonlinear eigenvalue problem. $\endgroup$
    – nicoguaro
    Sep 7, 2019 at 22:05
  • $\begingroup$ @KartMan, I have a Mathematica package that can solve this kind of eigenvalue problem, if it is of interest to you $\endgroup$
    – SPPearce
    Nov 12, 2019 at 8:20

1 Answer 1


Assume the equations are discretized on the $\tau$ grid, and introduce several column-vectors, using the notation where the superscript stands for the grid point index i $\in$ [1,...,n].

$\vec{\tau}=\left[ \tau^1,\tau^2,...,\tau^{n} \right] $

$\vec{\phi_1}=\left[ \phi_1^1,\phi_1^2,...,\phi_1^{n} \right] $

$\vec{\phi_2}=\left[ \phi_2^1,\phi_2^2,...,\phi_2^{n} \right] $

$\vec{f_{11}}=\left[ f_{11}^1,f_{11}^2,...,f_{11}^{n} \right] $

$\vec{f_{12}}=\left[ f_{12}^1,f_{12}^2,...,f_{12}^{n} \right] $

$\vec{f_{21}}=\left[ f_{21}^1,f_{21}^2,...,f_{21}^{n} \right] $

Next, the second derivative operator on the grid is a matrix $\hat{\Delta}$ of size $n \times n$.

Also, use $\alpha (k) = c_2(k) / c_1(k)$, and denote $c_1(k) = \lambda$ and $c_2(k) = \alpha \lambda$

Now, let's write both equations as a single linear system of size $2n \times 2n$ for the compound state column-vector $\vec{\phi}$ which is

\begin{equation} \vec{\phi} = [\phi_1^1,\phi_1^2,...,\phi_1^{n},\phi_2^1,\phi_2^2,...,\phi_2^{n}] \end{equation}

We'll use here several matrices of size $2n \times 2n$:

The left hand side matrix $L$ \begin{bmatrix} \hat{I} & \hat{0}\\ \hat{0} & \alpha \hat{I} \end{bmatrix}

the differential operator matrix $D$ \begin{bmatrix} \hat{\Delta} & \hat{0} \\ \hat{0} & \hat{\Delta} \end{bmatrix}

the cross-coupling matrix $C$ \begin{bmatrix} \hat{0} & \hat{I} \vec{f_{12}} \\ \hat{I} \vec{f_{21}} & \hat{0} \end{bmatrix}

and the forcing term matrix $F$ \begin{bmatrix} \hat{I} \vec{f_{11}} & \hat{0} \\ \hat{0} & \hat{0} \end{bmatrix}

Here $\hat{I}$ is the unit matrix of size $n \times n$, $\hat{0}$ is the null matrix of size $n \times n$.

Now our problem can be written as

\begin{equation} \lambda L \vec{\phi} = \left( D + C + F \right) \vec{\phi} \end{equation}


\begin{equation} \lambda \vec{\phi} = L^{-1} \left( D + C + F \right) \vec{\phi} \tag{*}\label{*} \end{equation}

This a linear eigenvalue problem, where the right hand side matrix contains the parameter $\alpha$, so the eigenvalues are functions of $\alpha$, which in turn is a function of $k$.

For a given value of $k$ there is a spectrum of eigenvalues, $\lambda_0, \lambda_1$ etc. Let's say we are interested in the smallest eigenvalue $\lambda_0$ which we'd just call $\lambda$. Solution of the eigenvalue problem $\eqref{*}$ by standard methods of linear algebra defines a function $\lambda(k)$, on the other hand by construction $\lambda$ must be equal to $c_1(k)$. The root of the equation $\lambda(k) = c_1(k)$ (if it exists) defines the solution of the problem. Since $\lambda(k)$ and $c_1(k)$ are two nonlinear functions the root of the equation $\lambda(k) = c_1(k)$ has to be sought by some iterative numerical technique.

  • $\begingroup$ How is this a linear eigevalue problem when both $\lambda$ and $L$ are functions of the unknown $k$ which needs to be determined? $\endgroup$
    – KartMan
    Sep 9, 2019 at 19:10
  • $\begingroup$ It is a linear eigenvalue problem for determining $\lambda$, as long as $\alpha$ is given. But if you need to find $k$ then you have to solve the nonlinear problem $c_1(k) = \lambda(k)$ using the above machinery. $\endgroup$ Sep 9, 2019 at 21:01
  • $\begingroup$ well, as I had described in the orginal post, $k$ and thus $c(k)$ are unknown and need to be determined. I don't see how this formulation casts the problem of finding $k$ into an eigenvalue/generalised eigenvalue problem. $\endgroup$
    – KartMan
    Sep 10, 2019 at 1:03
  • $\begingroup$ I added some explanations on how that could be done. To determine $k$ you have to solve a nonlinear equation; but for your example with a single ODE you'd have to do it as well, to solve $c(k)=\lambda$. $\endgroup$ Sep 10, 2019 at 2:57
  • $\begingroup$ if you can find $\lambda$ getting $k$ is quite easy in my case. There's no numerics required. Also, I had indicated in my earlier post that you can do some iterative (e.g. bisection method) procedure to get $k$ but this is not the eigenvalue formulation. I don't think what you wrote describes this as any eigenvalue problem. $\endgroup$
    – KartMan
    Sep 10, 2019 at 11:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.