Approximating the solution to a system of two ODEs with the KBM method?

Question

Background

I have the following system of ODEs:

$$ \begin{aligned} \dot x (t) &= \alpha - \beta x(t) y(t) \\ \dot y (t) &= \delta x(t) y(t) - \gamma y(t) \end{aligned} $$

where all parameters and variables are positive.

Setting the two equations equal to zero and solving for $x(t)$ and $y(t)$ reveals that the system exhibits an equilibrium, $(x^*,y^*)$, where $x^*>0$ and $y^*>0$. This equilibrium can exhibit damped oscillations under certain conditions.

Question

I was wondering if one could help me understand how to approximate the system's oscillatory solution using the KBM method? This post shows that, for small perturbations from $(x^*,y^*)$, such a system's solution can be well-approximated by its linearization. However, I am hoping to use the KBM method to acquire an approximation that works for larger perturbations.

Any help would be greatly appreciated!

Answer 1

The trick to using the Krylov–Bogoliubov method is to find a way to express at least one variable in your system—the one you expect to have the oscillating behavior—in terms of an uncoupled second-order differential equation.

Suppose we attempted to do so for $x$. The first thing we'd need to do is find an explicit form of $y$ in terms of $x$; which we can get right from the first equation in the system.

$$y = \frac{\alpha - \frac{dx}{dt}}{\beta x}$$

You can then substitute this into the second equation in your system for an explicit expression for $\frac{dy}{dt}$:

$$\frac{dy}{dt} = \frac{\delta \left(\alpha -\frac{dx}{dt}\right)}{\beta }-\frac{\gamma \left(\alpha -\frac{dx}{dt}\right)}{\beta x}$$

With this in hand, how do you get a second-order equation for $x$? Well, take the first equation for $\frac{dx}{dt}$, differentiate it with respect to time, substitute in the expression for $y$ and $\frac{dy}{dt}$ you pulled above, and hope you get something that fits a "standard" KBM format. If you do so, here's what you get:

$$\frac{d^2x}{dt^2} + \alpha \delta x= \alpha \gamma -\frac{\alpha \frac{dx}{dt}}{x}-\gamma \frac{dx}{dt}+\delta x \frac{dx}{dt}+\frac{(\frac{dx}{dt})^2}{x}$$

This certainly looks like a quasi-KBM-type problem; the left-hand side looks like a simple harmonic oscillator, and the right-hand side is chock-full of all the nonlinear stuff you'd expect to see in this kind of problem. The challenge is that it doesn't appear as if there is a single perturbation parameter that governs when the nonlinear content on the right-hand side becomes small; perhaps you'll get better results when attempting the same method for $y$.

Now suppose that you did finagle the equation above into the following "standard" KGB form for $x$:

$$\frac{d^2 x}{dt^2} + x = \epsilon f\left(x,\frac{dx}{dt},t\right)$$

The idea is to assume that your solution will be approximately of the following form:

$$x(t) \approx a(t) \cos(t + \psi(t)),\quad \frac{dx}{dt} \approx -a(t) \sin (t + \psi(t))$$

The motivation is that, when $\epsilon = 0$, both $a(t)$ and $\psi(t)$ will be constant in time, and the expectation is that they will vary relatively slowly in time for small values of $\epsilon$.

The KGB method then follows by inserting this ansatz for the solution in the equation described previously, and doing an involved but straightforward analysis to find differential equations for $a$ and $\psi$ in time. A full treatment of the process, and an example of its application on the Van der Pol system, can be found starting in page 19 of R. Rand's Lecture Notes on Nonlinear Vibrations.

Answer 2

(For typesetting economy I denote your system parameters by Roman letters. )

As a side comment on the previous answer note that the initial change of variables $ y = e^u$ transforms the associated 2nd-order equation into a fairly clean format:

$$\ddot u =(DA) - B \dot u e^u - BC e^u \tag{1}$$

The stationary solution you are interested in analyzing is located at $u_0= \ln(DA/BC)$, so $e^{u_0}= (DA/BC)$. Let $v = u-u_0$ be the deviation from the equilibrium solution.

After a change of time scale proportional to $\sqrt{ AD}$ chosen to normalize some of the numerical parameters, (1) becomes

$$\ddot v =(1- e^v) - Q \dot v e^v \tag{2}$$

where $Q= 1/C$.

Thus the system now really only depends on a single frictional damping quantity $Q$. At this point I suggest sketching some phase portraits to get visual intuition about the behavior of the system.

Some additional intuition. Multiply (2) by $\dot v$. Deduce that if $Q>0$ the quantity $E= (\dot v )^2 /2 + (e^v-v-1)$ satisfies $ \dot E= - \hat Q (\dot v)^2 e^v <0$, so the Liapunov "energy" $E$ declines steadily with time, each time a half-cycle semi-circular trajectory is completed, traveling from one null value of $\dot v$ to the next. This energy decay estimate implies strong stability of the equilibrium.

Next step. Suggest using $E = \frac{1}{2} [(\dot v)^2 + v^2] + v^3/3! + \ldots $ as a "radial" coordinate in the phase plane and obtaining the differential equation for it and some convenient auxiliary angular coordinate.