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?