Eigenvalue solution of a 2nd order ODE of a geophysical fluid dynamic problem

community. I am working on a project with a professor of mine and he suggested me to numerically solve a geophysical fluid dynamic instability ODE equation of a paper of Boccaleti et al. (2007) (equation 12).

This equation is written as follows: \begin{equation} \{1-[\sigma+kU(z)]^2\}w_{zz} - 2\left[\frac{k}{\sigma+kU(z)}-il\right]w_z - \left\{R_i|\mathbf{K}|^2[1-\delta^2(\sigma+kU(z))^2]+\frac{2ilk}{\sigma+kU(z)}\right\}w=0, \end{equation} where $U(z)$ is the basic state zonal velocity of the fluid, $R_i$ is the Richardson Number, $\mathbf{K} = k\mathbf{\hat{i}} + l\mathbf{\hat{j}}$ is the wavenumber vector of the instability (and $|K|= (k^2+l^2)^{1/2}$), $\sigma$ is the frequency of the instability and $\delta$ is the Prandtl Ratio.

Ok, the authors argue that this is an eigenvalue problem. In this case, $\sigma$ is the eigenvalue and $w$ is the eigenvector. So my approach to solve this equation is to manipulate this equation utilizing finite difference methods available and use an eigenvalue solution algorithm. However, I arrive at the following equation:

\begin{equation} (\mathbf{A}\sigma^3 + \mathbf{B}\sigma^2 + \mathbf{C}\sigma + D)w = 0, \end{equation}

which is a nonlinear eigenvalue equation. Does it makes sense? Was my reasoning accurate?

My second problem is with the lower boundary condition (where $z=-1$): \begin{equation} w = \frac{1}{\Delta B R_i |\mathbf{K}|^2}[(1-\sigma^2)w_z + (k\sigma+il)w]+ \frac{H_y}{|\mathbf{K}|^2}[(il-k\sigma^{-1})w_z - ikl\sigma^{-1}w]. \end{equation} If I substitute this BC into my main equation and follow the same reasoning I will have to deal with $\sigma^6$. The authors, however, state that this boundary condition is less tractable than what a previous author used. Therefore, they stated: "[...], so [equations] (5)-(9) are integrated forward in time to find dominant eigensolutions. [...]".

These equations they mentioned are a linearized set of QG equations:

\begin{equation} [\partial_t + ikU(z)]u + w -v + ikR_ip = 0 \\ \end{equation}

\begin{equation} [\partial_t + ikU(z)]v + u -v + R_ip_y = 0 \\ \end{equation}

\begin{equation} R_i\delta^2[\partial_t + ikU(z)]w - R_ib+R_ip_z = 0 \\ \end{equation}

\begin{equation} [\partial_t + ikU(z)]b - R_i^{-1}v + w = 0 \\ \end{equation}

\begin{equation} iku + v_y +w_z = 0 \end{equation}

Ok, so I conclude by asking two more questions:

What does it mean to integrate equations forward in time to find eigensolutions?

If I discretize (by finite differences) and then apply the BC to the main equation, arriving at something like: \begin{equation} (\mathbf{A}\sigma^6 + \mathbf{B}\sigma^5 + \mathbf{C}\sigma^4 + \mathbf{D}\sigma^3 + \mathbf{E}\sigma^2 + \mathbf{F}\sigma + \mathbf{G})w = 0, \end{equation} would it be right and feasible?

Please be patient with my writing, and sorry if I was not clear with something in this question.

Thank you very much

References:

Boccaletti, G., R. Ferrari, and B. Fox-Kemper, 2007: Mixed Layer Instabilities and Restratification. J. Phys. Oceanogr., 37, 2228–2250, https://doi.org/10.1175/JPO3101.1

Edit:

So I tried to use RK4 algorithm to solve the system of linearized QG equations I mentioned earlier. However, I am not as fluent in numerical methods as I wished so i might have made careless mistakes. So, what I did was the following. For the $u$, $v$, and $b$ variables I did use Runge Kutta. However, I have the boundary conditions for this specific system of equations being:

\begin{equation} w = -(\eta_t +H_y v) \quad z=-1 \end{equation}

\begin{equation} p = \Delta B \eta \quad z=-1 \end{equation}

\begin{equation} w = 0 \quad z=0 \end{equation}

where $\eta$ is the vertical displacement of the base layer, and $H_y$ is the meridional gradient of the base layer.

Using these conditions we can say:

\begin{equation} \int^0_{-1} w_z dz= \eta_t+H_yv = -i\int^0_{-1}ku+v_y \ dz \end{equation}

\begin{equation} \eta_t = -H_yv - i\int^0_{-1}ku+v_y \ dz \end{equation}

where I can then use RK4. For defining $w$, I did similar pattern:

\begin{equation} \int^0_{-z} w_z dz= w(0)-w(-z) = -i\int^{-z}_{-1}ku+v_y \ dz \end{equation}

\begin{equation} w(-z) = i\int^{-z}_{-1}ku+v_y \ dz \end{equation}

Finally, for $p$:

\begin{equation} \int^{-z}_{-1} p_z\ dz = p(-z) -\Delta B = -\delta^2[\partial_t + ikU(z)]w + b \\ \end{equation}

\begin{equation} p(-z) = \Delta B \eta -\delta^2[\partial_t + ikU(z)]w + b \\ \end{equation}

Is my intuition correct somewhat correct?

• Hi, I don't know if you managed to solve the nonlinear eigenvalue problem, but perhaps you will find useful to have a look at chebfun, which can deal with ODEs with parameters: mathworks.com/examples/matlab/community/… – Dmoreno Jul 3 '18 at 9:09
• Thank you for your answer. I did have tried to solve the system of linearized QG equations as I explained in my edit. However, I am unsure if it is correct. So, I will try to use this method to see if I arrive at an easier solution – J. Doe Jul 5 '18 at 18:06