Jacobian-free algorithm to calculate the phase sensitivity function in the phase reduction theory and its applications to K\'arm\'an's vortex street

TL;DR

This paper introduces a Jacobian-free algorithm for efficiently computing the phase sensitivity function in systems with limit cycles, especially applicable to incompressible fluid systems like Kármán's vortex street, reducing computational costs.

Contribution

The paper presents a novel Jacobian-free numerical method for calculating the phase sensitivity function, enabling applications to complex fluid systems with reduced computational effort.

Findings

The new algorithm significantly reduces computation time.

Successful application to Kármán's vortex street demonstrates practical utility.

Method extends phase reduction theory to incompressible fluid systems.

Abstract

Phase reduction theory has been applied to many systems with limit cycles; however, it has limited applications in incompressible fluid systems. This is because the calculation of the phase sensitivity function, one of the fundamental functions in phase reduction theory, has a high computational cost for systems with a large degree of freedom. Furthermore, incompressible fluid systems have an implicit expression of the Jacobian. To address these issues, we propose a new algorithm to numerically calculate the phase sensitivity function. This algorithm does not require the explicit form of the Jacobian along the limit cycle, and the computational time is significantly reduced, compared with known methods. Along with the description of the method and characteristics, two applications of the method are demonstrated. One application is the traveling pulse in the FitzHugh Nagumo equation in a periodic domain and the other is the K\'arm\'an's vortex street. The response to the perturbation added to the K\'arm\'an's vortex street is discussed in terms of both phase reduction theory and fluid mechanics.

Equations

$$\frac{d\mathbf{x}}{dt}=\mathbf{f}(\mathbf{x}),$$

$$\mathbf{Z}(\varphi )={\frac{\partial \Phi (\mathbf{x})}{\partial \mathbf{x}}|}_{\Phi (\mathbf{x})=\varphi }.$$

$$\frac{d\tilde{\mathbf{Z}}}{dt}=-{}^{t}J(\mathbf{p}(t))\tilde{\mathbf{Z}},$$

$$\tilde{\mathbf{Z}}(t)\cdot \mathbf{f}(\mathbf{p}(t))=1.$$

$$\Delta \Phi (\mathbf{x},\Delta \mathbf{x})=\Phi (\mathbf{x}+\Delta \mathbf{x})-\Phi (\mathbf{x}).$$

$$Z_j(\varphi )=\underset{ϵ\to 0}{\lim }\frac{\Delta \Phi (\mathbf{x},ϵ{\mathbf{e}}_j)}{ϵ},$$

$$J_{ij}(\mathbf{p}(t))=\frac{\partial f_i}{\partial x_j}(\mathbf{p}(t))=\underset{ϵ\to 0}{\lim }\frac{x_i(\mathbf{p}(t)+ϵ{\mathbf{e}}_j)-x_i(\mathbf{p}(t))}{ϵ},$$

$$J(\mathbf{p}(t))\mathbf{v}=\underset{ϵ\to 0}{\lim }\frac{\mathbf{x}(\mathbf{p}(t)+ϵ\mathbf{v})-\mathbf{x}(\mathbf{p}(t))}{ϵ}$$

$$\frac{d\mathbf{y}}{dt}=A(t)\mathbf{y},A(t)=J(\mathbf{p}(t)).$$

$$G_p(t)=\left({\mathbf{y}}_1(t)\cdots {\mathbf{y}}_n(t)\right).$$

$$\frac{d{\hat{G}}_p(t)}{dt}=-{}^{t}A{\hat{G}}_p(t),$$

$${}^t{\mathbf{z}}_p(t)(G_p(t+T)-G_p(t))=0.$$

$${}^t{\mathbf{z}}_p(t)(G_p(t)G_p(t+T{)}^{-1}-I)=0.$$

$$G_p(t+T)-G_p(t)=\left({\mathbf{g}}_1(t),{\mathbf{g}}_2(t),\cdots {\mathbf{g}}_n(t)\right).$$

$${}^t{\mathbf{z}}_p(t){\mathbf{g}}_k(t)=0.(k=1,\cdots n).$$

$$\mathbf{x}(t_0+T;t_0,\mathbf{p}(t_0)+{\mathbf{y}}_0)-\mathbf{x}(t_0+T;t_0,\mathbf{p}(t_0))-{\mathbf{y}}_0$$

$$\mathbf{z}=\mathbf{y}-\sum _{k=1}^{n-1}({\mathbf{a}}_k,\mathbf{y}){\mathbf{a}}_k.$$

$$\frac{\partial u}{\partial t}$$

$$\frac{\partial v}{\partial t}$$

$$\mathbf{U}$$

$$\mathbf{F}(\mathbf{U})$$

$$(g_k^u,g_k^v)$$

$$(Q_u(x;\varphi ),Q_v(x;\varphi ))={\left(\frac{\delta \Phi (u,v)}{\delta u},\frac{\delta \Phi (u,v)}{\delta v}\right)|}_{(u,v)=(u_p(x,\frac{\varphi }{\omega }),v_p(x,\frac{\varphi }{\omega }))},$$

$$\omega$$

$$\simeq$$

$$\mathbf{Z}(\varphi )=(Z_1^u,\cdots ,Z_N^u,Z_1^v,\cdots ,Z_N^v),Z_k^u=\frac{\partial \varphi }{\partial u_k},Z_k^v=\frac{\partial \varphi }{\partial v_k}(k=1,\cdots ,N).$$

$$(Q_u(k\Delta x;\varphi ),Q_v(k\Delta x;\varphi ))=\frac{\omega }{\Delta x}(Z_k^u(\varphi ),Z_k^v(\varphi )).$$

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}=-\nabla p+\frac{1}{Re}\Delta \mathbf{u},\nabla \cdot \mathbf{u}=0,$$

$$\mathbf{Q}(\mathbf{x};\varphi )=(Q_u,Q_v),Q_u(\mathbf{x};\varphi )={\frac{\delta \Phi (\mathbf{u})}{\delta u}|}_{\mathbf{u}={\mathbf{u}}_p(\mathbf{x},\varphi /\omega )},Q_v(\mathbf{x};\varphi )={\frac{\delta \Phi (\mathbf{u})}{\delta v}|}_{\mathbf{u}={\mathbf{u}}_p(\mathbf{x},\varphi /\omega )},$$

$$\Phi (\mathbf{u}(\mathbf{x},t)+{\mathbf{u}}^{\prime }\delta (\mathbf{x}-{\mathbf{x}}_0))-\Phi (\mathbf{u}(\mathbf{x},t))\simeq \int \frac{\delta \Phi (\mathbf{u})}{\delta \mathbf{u}}\cdot {\mathbf{u}}^{\prime }\delta (\mathbf{x}-{\mathbf{x}}_0)d\mathbf{x}=\mathbf{Q}({\mathbf{x}}_0;\varphi )\cdot {\mathbf{u}}^{\prime },$$

$$\omega$$

$$\simeq$$

$$\mathbf{Q}({\mathbf{x}}_{i,j};\varphi )=\frac{\omega }{\Delta S_{i,j}}(Z_{i,j}^u,Z_{i,j}^v).$$

Full text

aaa

††preprint: aaa/bbb

I Introduction

Periodic flow is the simplest type of unsteady flow. A well-known example is the Kármán’s vortex street observed in the downstream of a cylinder in a uniform flow, when $Re\gtrsim 50$, where $Re=UD/\nu$; $U$ is the velocity of the uniform flow, $D$ is the diameter of the cylinder, and $\nu$ is the kinematic viscosityWilliamson (1996). If the whole system is described in the phase space, the periodic flow is represented by the stable limit cycle (LC), a periodic orbit with a unique period such that all the orbits near the LC converge to LCStrogatz (1994).

The state near the LC can be described by a single variable called “phase”, $\phi$, and the dynamics near the LC can be described by the ordinary differential equation (ODE) of $\phi$, which is a significant reduction of the degree of the freedomKuramoto (1984) . The phase reduction technique has been applied to many problems, including mechanical vibrations, spiking neurons, and flashing firefliesNakao (2015). In the phase reduction theory, the phase sensitivity function, which gives the linear response coefficients of the phase to perturbations, is essential for the reduction of the original system into the ODE system of $\phi$ Kuramoto (1984); Nakao (2015). A practical method to calculate the phase sensitivity function, “the adjoint method”Ermentrout (1996), has been applied to many problems not only in the systems described by the ordinary differential equationsNakao (2015), but also the partial differential equations, including the reaction-diffusion systemNakao et al. (2014), convection in Hele-Shaw cellKawamura and Nakao (2015), and flagellum synchronization of a model equationKawamura and Tsubaki (2018).

If the phase reduction theory is applied to incompressible fluid systems, not only are new aspects of the periodic flows obtained, but also new techniques of the flow control will be developed. In the case of the Kármán’s vortex street, we will be able to design an efficient perturbation form to change the phase of the flow (or the timing of the separation of vortices). Further, the synchronization of two cylinders Williamson (1985); Peschard and Gal (1996); Akinaga and Mizushima (2005) can be analyzed.

However, the phase reduction theory has not been applied to incompressible fluids, except for the cases in which the linearized equation around the LC can be obtained explicitlyKawamura and Nakao (2015) and the case of direct measurements of the phase shift by adding perturbation to the flow to obtain the information of the prescribed regionsTaira and Nakao (2018). A practical problem exists, if the phase reduction theory is applied to address incompressible fluids, which is, a computational source. The description of the states near the LC requires the linearized matrix (the Jacobian) of any point along the LC; further, the explicit form of the matrix cannot be obtained because the Poisson equation must be solved to obtain the pressure. In such a case, we need to numerically calculate all the elements of the Jacobian and store the values, which requires a lot of memory. Recently, Taira and NakaoTaira and Nakao (2018) calculated the phase sensitivity function at the top of the separation point of the cylinder in a uniform flow by a direct method. However, a computational source is needed to obtain the phase sensitivity function of the whole region. Further, it requires a lot of periods for the perturbed system to converge to LC. They will be discussed later.

In this paper, we propose a new numerical algorithm to calculate the phase sensitivity function, which is applicable for the systems with large degree of freedom and without explicit expression of the Jacobian. In particular, we applied this method to analyze the Kármán’s vortex street. The obtained phase sensitivity function revealed that the downstream region comprises narrow bands of the effective area for the phase shift, and the distribution and effective directions are time-dependent. Moreover, the distribution in the upstream region is less time-dependent, which suggests that controlling the phase is more convenient. Furthermore, a comprehensive interpretation based on fluid mechanics is discussed.

The remainder of this paper is organized as follows. In Sec. II, we describe the details of the proposed method and discuss the characteristics. In Sec. III, we demonstrate the proposed method. In Sec. III.1, the phase sensitivity function of the traveling pulse is compared with known results. In Sec. III.2, we exhibit the phase sensitivity function of the Kármán’s vortex street and discuss the detailed characteristics in terms of fluid mechanics. In Sec. IV, we summarize the results.

II Methods

In this section, the proposed method is described in detail to numerically calculate the phase sensitivity function of a dynamic system with finite dimensions. In the case of the partial differential equations (PDE), the discretized system is considered.

II.1 Phase reduction

The definitions and notations of the phase reduction theory are briefly summarized in this section. Refs.Nakao (2015); Kuramoto (1984) contains more information on this theory. Let us consider the $n$-dimensional autonomous dynamical systems given by:

[TABLE]

where the vector $\bm{x}={}^{t}(x_{1},\cdots,x_{n})\in\mathbb{R}^{n}$ represents the state, and the function $\bm{f}(\bm{x})={}^{t}(f_{1}(\bm{x}),\cdots,f_{n}(\bm{x})):\mathbb{R}^{n}\mapsto\mathbb{R}^{n}$ determines the dynamics. It is assumed that eq. (1) has a stable limit cycle solution $\bm{x}(t)=\bm{p}(t)$ that satisfies $\bm{p}(t+T)=\bm{p}(t)$ for all $t$, where $T$ is the natural period.

In the phase reduction theory, a value called “phase”, $\phi$, (more precisely, the “asymptotic phase” in Ref.Nakao (2015)) is defined on and near the limit cycle (LC) as follows: on the LC, the origin of the phase ($\phi=0$) is chosen, and the orbit $\bm{p}(t)$ that starts from the origin at $t=0$ is considered; then, $\phi$ is defined as $\displaystyle\phi(t)=\omega t\;(\mod 2\pi)$, where $\omega=2\pi/T$ and $0\leq\phi(t)<2\pi$. For the phase at a particular point $\bm{x}_{0}$ near the LC, $\Phi(\bm{x}_{0})$, the orbit $\bm{x}(t)$ that starts from $\bm{x}=\bm{x}_{0}$ at $t=0$ is considered. Now, $\Phi(\bm{x}_{0})=\phi_{0}$, in which, $\phi_{0}$ satisfies $\lim_{t\to\infty}(\bm{x}(t)-\bm{p}(t+\phi_{0}/\omega))=0$.

A fundamental function of the phase reduction theory is the phase sensitivity function, $\bm{Z}(\phi)$, which is defined as:

[TABLE]

The phase sensitivity function determines the phase increment of the state near the LC due to a small perturbation $\Delta\bm{x}$, because of the relationship $\displaystyle\Phi(\bm{x}+\Delta\bm{x})=\bm{Z}(\phi)\cdot\Delta\bm{x}+O(|\Delta\bm{x}|^{2})$. Once $\bm{Z}(\phi)$ is obtained, the extent of phase shift after perturbation due to any force, and the synchronized property of the coupled oscillators, can be calculatedWinfree (1967); Kuramoto (1984); Nakao (2015).

If $\displaystyle\bm{\tilde{Z}}(t)=(1/\omega)\bm{Z}(\omega t)$ is defined, $\bm{Z}(\phi)$ is obtained as the periodic solution of the adjoint equation as follows:

[TABLE]

where $J$ is the Jacobian of the dynamical system (1), and its $(i,j)-$component is $J_{ij}=\frac{\partial f_{i}}{\partial x_{j}}(\bm{p}(t))$ Ermentrout (1996). Because the adjoint equation is linear, an additional condition to determine the amplitude of $\tilde{\bm{Z}}(t)$ is required. The following normalization relationship is imposed:

[TABLE]

Once $\tilde{\bm{Z}}(t)$ is obtained, $\bm{Z}(\phi)$ can be easily obtained by converting $t$ to $\phi$.

II.2 Problems in calculations of the phase sensitivity function in the incompressible fluid system

Let us consider the case that the dimension $n$ is large and the function $\bm{f}(\bm{x})$ is not given explicitly, which is the case of the discretized system of the incompressible fluid system for numerical calculation.

A simple way to calculate the phase shift due to the perturbation is called “direct method”Nakao (2015); Taira and Nakao (2018). A phase shift $\Delta\Phi$ due to the perturbation $\Delta\bm{x}$ added to the state $\bm{x}$ is given by:

[TABLE]

Then, the phase sensitivity function $\bm{Z}(\phi)=(Z_{1}(\phi),\cdots,Z_{n}(\phi))$ is obtained by:

[TABLE]

where $\bm{e}_{j}$ is a unit vector and its $k-$th component is $[\bm{e}_{j}]_{k}=\delta_{jk}$ ($\delta_{jk}$ is the Kronecker’s delta).

To evaluate the calculation time by the direct method, $\Delta\tau$ is assumed as the time to calculate one time step of the dynamical system (1) per one degree of freedom, i.e., time $n\Delta\tau$ is needed to calculate single time step of the whole system. It is further assumed that the time step is given by the $T/m$, where $m$ is the division number of the period $T$ and that it needs $N_{1}$ periods for the system to converge to estimate the phase shift $\Delta\Phi$. Then, $N_{1}m$ steps (equivalently, $N_{1}$ periods) are needed to calculate the phase shift due to the single perturbation $\Delta\bm{x}=\epsilon\bm{e}_{j}$. Therefore, the total calculation time is estimated as $N_{1}mn^{2}\Delta\tau$ because $n$ perturbations $\left\{\epsilon\bm{e}_{j}\mid j=1,\cdots,n\right\}$ are needed to obtain $\bm{Z}(\phi)$.

If the adjoint method is applied, an asymptotic state of the adjoint equation (3) must be obtained as $t\to-\infty$. When applying this procedure to the incompressible fluid system, several problems occur. In this case, the analytic expression of the Jacobian $J$ cannot be obtained because the pressure must be calculated by solving the Poisson equation. Thus, we need to calculate $J$, $n\times n$ matrix, every time step on $\bm{x}(t)=\bm{p}(t)$. The memory required to store the Jacobians along the LC is estimated of the order of $mn^{2}$. For example, approximately 800 GB is required to store $J(\bm{p}(t))$ with double precision, when $m=1,000$ and $n=10,000$, which is too large for the calculation.

Alternatively, the components of $J(\bm{p}(t))$ can be evaluated at each time step. If the following formula is used,

[TABLE]

the calculation time of $J$ at a particular time is of the order $n^{2}\Delta\tau$. If the adjoint system takes $N_{2}$ periods to converge, we need an order of time $N_{2}mn^{2}\Delta\tau$, which is the same order of the calculation time for the direct method if we can assume that $N_{1}$ and $N_{2}$ are of the same order. These methods require long-time integration of the original system or the adjoint system until convergence. This process can be reduced by the proposed method, which is discussed in the next subsection.

II.3 Proposed method

We propose a novel method to calculate the phase sensitivity function that can reduce the computational cost by the factor $1/N_{1}$ or $1/N_{2}$.

In this method, a Jacobian-free method is utilized to calculate the product $J\bm{v}$, where $\bm{v}$ is a vector, by the following formula:

[TABLE]

Knoll and Keyes (2004). This method can minimize the time required because the $n$ components of $J\bm{v}$ can be calculated at one time. If we need to calculate all the components of $J$ by eq. (7), the computational time of $J\bm{v}$ is of the order of $n^{2}\Delta\tau$, while it is $n\Delta\tau$ when eq. (8) is used. The formula (8) cannot be applied for the calculation of the adjoint equation (3) because no similar formula has been known for ${}^{t}J\bm{v}$ Govaerts (2000).

Now, we consider a method to obtain the periodic solution of the adjoint equation (3). The linearized equation of (1) from the limit cycle $\bm{p}(t)$ can be written as:

[TABLE]

The linearly independent set of the solution of (9), $\bm{y}_{1},\cdots\bm{y}_{n}$, can be used to define the fundamental solution matrix as follows:

[TABLE]

Then, the general solution $\bm{y}(t)$ is represented by $\bm{y}(t)=G_{p}(t)\bm{c}$, where $\bm{c}$ is a constant column vector determined by the initial condition.

Differentiating the identity $G_{p}(t)^{-1}G_{p}(t)=I$ with respect to $t$ and using the relationship $\displaystyle\frac{dG_{p}}{dt}=AG_{p}$, we obtain the following equation:

[TABLE]

where $\hat{G}_{p}(t)={}^{t}(G_{p}(t)^{-1})$. Equation (11) shows that $\hat{G}_{p}(t)$ is the fundamental solution matrix of the adjoint equation (3). Let us assume that Eq. (3) has a unique limit cycle solution $\bm{z}_{p}(t)$. Then, the periodicity condition, $\bm{z}_{p}(t+T)=\bm{z}_{p}(t)$ can be reduced to:

[TABLE]

Eq. (12) can be shown as follows: If we write $\bm{z}_{p}(t)=\hat{G}_{p}(t)\bm{c}$, the periodicity condition is, $\hat{G}_{p}(t+T)\bm{c}=\hat{G}_{p}(t)\bm{c}$. Taking the transpose of this equation and using the identity $G_{p}(t)^{-1}G_{p}(t)=I$, either ${}^{t}\bm{c}G_{p}(t)^{-1}G_{p}(t)G_{p}(t+T)^{-1}=\bm{c}{}^{t}G_{p}(t)^{-1}$, or the following equation can be obtained:

[TABLE]

This equation is equivalent to Eq. (12). Solving eq. (12), we can obtain $\bm{z}_{p}(t)$, which is proportional to $\tilde{\bm{Z}}(t)$.

The calculation procedure of $\bm{z}_{p}(t)$ using Eq. (12) is as follows. Let us rewrite $G_{p}(t+T)-G_{p}(t)$, using column vectors $\bm{g}_{k}(t)\;(k=1,\cdots,n)$ as follows:

[TABLE]

Then, the equation (12) is decomposed to $n$ orthogonality relationships between $\bm{z}_{p}$ and $\{\bm{g}_{k}\mid k=1,\cdots,n\}$, i.e.,

[TABLE]

Further, the vectors $\{\bm{g}_{k}\}$ can be obtained through time integration of the original system. Let us write $\bm{x}(t;t_{0},\bm{x}_{0})$ as the solution of Eq. (1) with $\bm{x}(t_{0})=\bm{x}_{0}$. Then, for any small vector $\bm{y}_{0}$, we obtain the following relationship, if higher order terms are negligible:

[TABLE]

This formula can be used to calculate $\bm{g}_{k}(t_{0})$ by setting $\bm{y}_{0}=\epsilon\bm{e}_{k}$, where $\epsilon$ is a small parameter. Because of the existence and uniqueness of the periodic solution, one eigenvalue of $G_{p}(t_{0}+T)-G_{p}(t_{0})$ is zero. Therefore, $\bm{g}_{1},\bm{g}_{2},\cdots\bm{g}_{n}$ are linearly dependent, and $(n-1)$ elements of $\{\bm{g}_{k}\}$ are linearly independent. Thus, there exists a non-trivial solution of (12).

The algorithm to find the solution is as follows: using the Gram-Schmidt orthonormalization and applying the relationship (16) with the linearly independent set of $\bm{y}_{0}$, we can obtain $(n-1)$ vectors $\{\bm{a}_{k}\mid k=1,\cdots,n-1\}$ which construct the orthonormal basis of the linear space $V$ spanned by $\{\bm{g}_{k}\mid k=1,\cdots,n\}$, i.e., $(\bm{a}_{i},\bm{a}_{j})=0\;(i\neq j,1\leq i,j\leq n-1)$ and $|\bm{a}_{i}|=1$. Finally, using a general vector $\bm{y}\in\mathbb{R}^{n}$ such that $\bm{y}\notin V$, the following $\bm{z}$ is obtained:

[TABLE]

By definition, $\bm{z}$ is perpendicular to any vector of $V$, i.e., $(\bm{z},\bm{a}_{i})=0\;(i=1,\cdots,n-1)$. Thus, $\bm{z}$ is the solution of eq. (15) and eq. (3).

II.4 Characteristics

The characteristics of the proposed method is as follows:

First, this method is memory-saving and Jacobian-free, i.e., an explicit expression of the Jacobian is not needed. This means that we do not need memory to store the Jacobian data over one period, which requires $mn^{2}$ variables. Second, this method is time-saving. This method enables the calculation of the phase sensitivity function by $n-1$ calculation of time evolution over one period, $T$. The total computation time is estimated as $mn(n-1)\Delta\tau\sim mn^{2}\Delta\tau$, $1/N_{1}$ or $1/N_{2}$ of the direct integration of the adjoint equation. Third, this method can be efficiently implemented by parallel computation. The bottleneck of this method is the calculation of $n-1$ vectors in the form $(G_{p}(t+T)-G_{p}(t))\bm{y}_{0}$. Each calculation needs time integration over one period, which requires more computation time. However, each process can be independently calculated in parallel. This means that the parallel computing algorithms such as MPI method work efficiently. These three characteristics are the advantages of the proposed method.

In addition, the following remarks should be considered. First, a periodic solution data $\bm{p}(t)$ must be prepared for the calculation; there are several algorithms to obtain the periodic solution numerically, e.g., Ref.Saiki (2007); Second, the proposed algorithm gives $\bm{Z}(\phi)$ for single phase, which appears a disadvantage initially. When we consider the traveling wave in a periodic domain, e.g., traveling pulse in the FitzHugh Nagumo(FHN) equationNakao et al. (2014), the phase sensitivity function $\bm{Z}(\phi)$ at one phase is sufficient because of the Galilean invariance of the solution. This example will be discussed in the next section. In such cases, the problem does not need to be considered. Generally, this problem can be amended, because the LHS of eq. (14) for different value of $t$ can be calculated by the data of $G_{p}(t)$ over two periods, which does not change the order of the calculation cost.

III Applications

The proposed method is applied to two PDE problems. First, the phase sensitivity function of a traveling pulse is calculated in the FHN equation in a periodic domain to compare with the results given in Ref.Nakao et al. (2014). Second, the phase sensitivity function of the Kármán’s vortex street is calculated.

III.1 A traveling pulse in FitzHugh Nagumo equation in a periodic domain

Nakao, Yanagita, and Kawamura developed phase reduction theory for the reaction-diffusion system, and calculated the phase sensitivity function for several solutions including the traveling pulse of FHN model (sec. IIIA in their paper). The FHN equations are described as:

[TABLE]

in one dimensional space. Here $u(x,t)$ and $v(x,t)$ are independent variables and $\alpha,\gamma,\tau,\kappa,$ and $\delta$ are constants. The details are explained in Sec. IIIA in Ref.Nakao et al. (2014).

When the system parameters are $\alpha=0,\tau^{-1}=0.018,\gamma=1,\kappa=1,$ and $\delta=0.02$, which are the same values in Ref.Nakao et al. (2014), a traveling pulse solution with wavy tails exists. We calculated the FHN equation for a one-dimensional periodic domain with the size $L=300$ with spatial grids $\Delta x=0.5$; the space was discretized to $N(=L/\Delta x)$ discrete points. The time integration of the FHN equations was calculated using the Runge-Kutta method of the second order. The periodic solution was obtained by the Newton-Raphson methodSaiki (2007), in which the GMRES(k) method with the Jacobian-free algorithmKnoll and Keyes (2004) was used for solving the linear equation. The period was divided into $m(=6,000)$ discrete states. The period of the obtained solution was $T=553.15$.

The phase sensitivity function was calculated for the discretized system with $2N$ dimensions, where:

[TABLE]

with the periodic boundary condition $u_{j+N}=u_{j},v_{j+N}=v_{j}$ for any $j$. The discretized equations for Eqs. (18) and (19) are represented as $\displaystyle\frac{d\bm{U}}{dt}=\bm{F}(\bm{U})$. The phase sensitivity function was converted to the phase sensitivity function for the reaction diffusion system, defined by (B13) in Ref.Nakao et al. (2014), given by:

[TABLE]

where $(u_{p}(x,t),v_{p}(x,t))$ is the traveling pulse solution which is time periodic. The difference between $(Q_{u}(x;\phi),Q_{v}(x;\phi))$ and $\bm{Z}(\phi)$ is just a numerical factor when the grid distance $\Delta x$ is homogeneous. Normalization of $(Q_{u}(x;\phi),Q_{v}(x;\phi))$ is given by (B9) in Ref. Nakao et al. (2014), which is:

[TABLE]

The normalization condition for $\bm{Z}(\phi)$, eq. (4) reads $\sum_{k=1}^{N}\left\{Z_{k}^{u}(\phi)g_{k}^{u}+Z_{k}^{v}(\phi)g_{k}^{v}\right\}=1$, where:

[TABLE]

Eqs. (25) and (27) produce the following relationship:

[TABLE]

In Fig. 1, the traveling pulse solution and corresponding phase sensitivity functions $Q_{u}$ and $Q_{v}$ at a phase are shown. The shapes of both the traveling pulse solution $(u_{p},v_{p})$ and the phase sensitivity functions are similar to Fig. 1(a) in Ref.Nakao et al. (2014), which validates our proposed method. In this case, the solution is a traveling one, i.e., $u(x,t)=u(x-ct,0)$ and $v(x,t)=v(x-ct,0)$, where $c$ is the speed of the traveling pulse. Let us assume that the phase sensitivity functions in Fig. 1 is at $\phi=0$. Then, we have the relationships $Q_{u}(x,\phi)=Q_{u}(x-c\phi/\omega,0)$ and $Q_{v}(x,\phi)=Q_{v}(x-c\phi/\omega,0)$. Thus, the phase sensitivity function at another phase can be obtained by the spatial translation of $Q_{u}$ and $Q_{v}$.

Here, the traveling pulse propagates to the right. Thus, the perturbation in the right area of the pulse interacts with the entire wavy tail. Such interaction causes a relatively significant phase shift. In addition, due to the wavy characteristics, both $Q_{u}$ and $Q_{v}$ oscillate spatially. Moreover, the perturbation in the left area of the pulse only interacts with a part of the wavy tail and the phase shift is less significant. In this sense, the right area is the “upstream” of the pulse. The wavy shapes of $Q_{u}$ and $Q_{v}$ on the right of the pulse match these observations.

III.2 Kármán’s vortex street

III.2.1 Methods

We consider the flow past the cylinder in a uniform flow in two-dimensional space. The flow is governed by the incompressible Navier-Stokes equations in the non-dimensional form:

[TABLE]

where $\bm{u}=(u,v)$ is the velocity, $p$ is the pressure, and $Re$ is the Reynolds number.

The diameter of the cylinder is unity, and the uniform flow is represented by $\bm{U}=(1,0)$. The computational domain is a circle of radius $R$ and the center, which is also the center of the cylinder, is the origin of the coordinate.

The boundary conditions in the polar coordinate $(r,\theta)$ are given as follows: the non-slip condition ($\displaystyle\bm{u}=\bm{0}$ and $\displaystyle\frac{\partial p}{\partial r}=0$) was applied to the cylinder ($\displaystyle r=\frac{1}{2}$). On the outer boundary ($\displaystyle r=R$), the inflow condition ($\displaystyle\bm{u}=(1,0)$ and $\displaystyle\frac{\partial p}{\partial r}=0$) was applied in the region $\displaystyle\frac{\pi}{3}<\theta<\frac{5\pi}{3}$, and the outflow condition ($\displaystyle\frac{\partial\bm{u}}{\partial r}=\bm{0}$ and $p=0$) was applied in the region $\displaystyle 0\leq\theta<\frac{\pi}{3},\frac{5\pi}{3}<\theta<2\pi$.

The computational domain was discretized to $N_{r}\times N_{\theta}$ grids, where $N_{r}$ is the division number in the radial direction and $N_{\theta}$ is the division number in the azimuthal direction. The grid was constructed such that the grid spaces were finer near the cylinder and in the downstream area (Fig. 2). For computation, the Navier-Stokes equation was discretized by the finite volume method. The advection term was calculated using the flux splitting methodLiu and Kawachi (1998) with flux of the third order at the boundary of the control volume, and the dissipation term was calculated by the Crank-Nicolson method. The linear equations for the Poisson equation to obtain the pressure and the dissipation term were numerically solved by the BiCGSTAB method using the open software Lis (https://www.ssisc.org/lis/). For time integration, the Adams-Bashforth method was used. The periodic solution was obtained using the same algorithm as that applied for the traveling pulse of the FHN eqs in Sec.III.1.

The computational parameters are $N_{r}=N_{\theta}=60$ and $R=10$. In this condition, the radial grid width ranges from $0.00840$ to $0.308$, and the azimuthal grid width from $0.0525$ to $0.157$. The Reynolds number was $Re=200$. The phase sensitivity functions were calculated with coarser mesh, $N_{r}=N_{\theta}=40$, and found no substantial difference between the results with these meshes.

To calculate the periodic solution representing the Kármán’s vortex street, the period was discretized to 1,000 steps. The time origin $t=0(\phi=0)$ was set at the minimum of the lift coefficient. Eight snapshots of the obtained periodic solution are shown in Fig. 3. The periodic solution gives the mean drag coefficient $\langle C_{D}\rangle=\langle 2F_{x}\rangle=1.346$ and the Strouhal number $St=fD/U=0.201$. These values are close to the values of previous studiesHenderson (1995); Williamson (1996).

Fig. 4 shows the time series of $C_{L}$ for the periodic solution and those started from the perturbed states, to demonstrate the occurrence of phase shift. A perturbation was applied to the single horizontal velocity component of the periodic solution. The position $\bm{x}_{pert}$ was set $(-6.00,0)$, a point in the upstream of the cylinder, which is indicated by the red point in Fig. 2. To perturb the velocity, the discretized velocity component was changed at $\bm{x}=\bm{x}_{pert}$ as $u\mapsto u\pm\epsilon_{0}$ at $t=0$, where $\epsilon_{0}=0.5$, which indicates that the velocity in the corresponding control volume (area) is changed. The time series of $C_{L}$ in the cases with the perturbations $\pm\epsilon_{0}$ and that in the case with no perturbation are shown in Fig. 4, where the cases of $+\epsilon_{0}$ and $-\epsilon_{0}$ cause the advance and delay of the phase, respectively. The changes in the phase and amplitudes of $C_{L}$ become apparent in the time regime $t\gtrsim 4$; the perturbation does not cause significant changes in $C_{L}$ before it reaches the cylinder. Further, amplitude changes rather than phase changes are apparent during $4\lesssim t\lesssim 8$, which implies that a diffused perturbation slightly enhanced (reduced) the flow speed around the cylinder to cause an increase (decrease) in the lift force. After the perturbation was advected downstream, the periodic state of the Kármán’s vortex street recovers whereas the phase shift remained. These behaviors will be discussed in Sec. III.2.2 in detail with the results of the phase reduction theory.

III.2.2 Phase sensitivity functions for Kármán’s vortex street

The proposed method enables us to produce the phase sensitivity function of Karman’s vortex street or the LC of the Navier-Stokes equation. For the following discussion, the phase sensitivity vector for the fluid system, $\bm{Q}(\bm{x};\phi)$, is defined as

[TABLE]

where $\bm{u}_{p}(\bm{x},t)$ is the time periodic solution corresponding to the Karman’s vortex street. The components of the phase sensitivity vector are the phase sensitivity function defined in Ref. Nakao et al. (2014). The phase sensitivity vector causes the phase shift due to a point-wise perturbation to the velocity field. By definition of the functional derivatives, the following relationship holds:

[TABLE]

when $|\bm{u}^{\prime}|$ is small. The normalization of the phase sensitivity vector is given by:

[TABLE]

where $\bm{x}_{i,j}$ is the representative position of the control volume(area) indicated by $(i,j)$ ($i$ and $j$ are the radial index and azimuthal index, respectively), and $\Delta S_{i,j}$ is the area of the control volume. To estimate the phase sensitivity vector $\bm{Q}(\bm{x};\phi)$, $(Z_{i,j}^{u},Z_{i,j}^{v})=(\partial\phi/\partial u_{i,j},\partial\phi/\partial v_{i,j})$ is defined, where $(u_{i,j},v_{i,j})$ is the discretized velocity component at the position indicated by $(i,j)$. Then, the following formula which is similar to Eq. (27) is obtained:

[TABLE]

In Figs. 5 and 6, $\bm{Q}(\bm{x};\phi)$ is shown downstream of the cylinder and the wider region including the upstream region, respectively.

Downstream of the cylinder, a large value of $|\bm{Q}|$ is observed in the area whose size is approximately twice the diameter of the cylinder (green area in Fig. 5). One or two band(s) with particularly large $|\bm{Q}|$ ($|\bm{Q}|>6$) are observed in the downstream region of the cylinder, which change its shape with time (or phase).

The integral curves of $\bm{Q}$ in the large-$|\bm{Q}|$ region constructs several closed curves or spirals, which are referred to as “Q-eddy” henceforth. The direction of the large-$|\bm{Q}|$ band shares the edge of Q-eddy. During a single period, the Q-eddy is generated near the cylinder, and is transferred to the downstream and disappears. At $t=0$, two Q-eddies $Q_{1}^{+}$ and $Q_{2}^{+}$, both of which are counterclockwise, exist near the cylinder; this direction can be defined as positive and vice versa. These Q-eddies are transferred to the downstream to merge with the single Q-eddy $Q_{12}^{+}$ at $t=2/8T$. At the same time, new two Q-eddies $Q_{3}^{-}$ and $Q_{4}^{-}$ were detached from the cylinder. They are also transferred to the downstream to form $Q_{34}^{-}$ ($t=6/8T$).

When the large-$|\bm{Q}|$ band is compared with the flow speed distribution in Fig. 3, the band corresponds to the region where the edge of the low flow speed region at the back of the cylinder. Further, the direction corresponds to the flow vector of the eddy which is attached to the cylinder and will be detached (e.g. lower eddy in Fig.3, $t=0$). This observation matches our expectation that such a perturbation supporting the vortex generation results in earlier separation leading to the phase advance. However, it must be noted that the structure depends significantly on both time and space, and the application of the phase control in this region requires regulating the perturbation distribution to the velocity field.

In the upstream of the cylinder, the region with large $|\bm{Q}|$ spreads wider, although the peak values are not significantly high as those in the downstream (Fig. 6). The region has a triangular shape. A significant feature is that the $x-$component of $\bm{Q}$ in the direct upstream region of cylinder, which may be characterized by $|y|<D/2$ and $x<0$, is always positive. The phase shift due to the perturbation in this region is demonstrated in Sec.III.2.1 (Fig.4).

These characteristics of the $\bm{Q}$ field can be explained as follows: The positive perturbation (e.g. $\bm{u}^{\prime}=(\epsilon,0)$, where $\epsilon$ is a small positive number) in the upstream dissipates to spread during transfer to the downstream, which enhances the local speed near the cylinder. Because the period of the LC for the Karman’s vortex street is related to the uniform flow by a constant Strouhal number, $St=fD/U$, the frequency $f$ is slightly enhanced due to a slight enhancing of the local speed. After the perturbation is transferred downstream, the state converges to the LC, but the existence of a time interval with large frequency shifts the phase in advance. For this scenario, the perturbation position should be located at a certain distance from the cylinder; such distance is needed to spread the perturbation so that it can be regarded as an enhancement of the local speed around the cylinder. These characteristics of the $\bm{Q}$ field suggest that a simple control strategy is possible by using the direct upstream region of the cylinder than other regions.

The Q-eddies in the upstream are located at a constant interval on both the sides of the direct upstream region (e.g. $Q_{6}^{-}$ and $Q_{7}^{-}$ in Fig. 6). They travel downstream as the time (phase) increases (for example, $Q_{5}^{+}$ in Fig. 6(a)-(c)). The spiral structure implies that the sign of the phase shift due to the constant perturbations in this area changes with time. The perturbation added to this region (upper or lower regions of the direct upstream region) is transferred to the side of the cylinder. Therefore, it works to change the separation timing rather than to enhance the flow speed around the cylinder. Thus, the Q-eddy structure should be spatially periodic with the period length roughly estimated by $TU\simeq 5$.

It is important to note that such property of the Q-eddy structure in the upstream is similar to the result in Sec. III.1, where $Q_{u}$ and $Q_{v}$ of the traveling pulse oscillated in the “upstream” of the pulse with oscillatory tail, though the “upstream” for the traveling pulse must be interpreted as the region where the pulse travels.

The phase sensitivity function calculated by the proposed method was compared with the result obtained by the direct method. Several points in the upstream shown in Fig. 7 were selected, which were used for the comparison.

The procedure of the direct method is as follows: the field at $t=0(\phi=0)$ was chosen, and the discretized velocity component at each selected point, $(u,v)$, was perturbed by either $u\mapsto u+\epsilon_{0}$ or $v\mapsto v+\epsilon_{0}$. The calculation of the time integration started from the perturbed field to obtain the time series of $C_{L}$. Let us define $f_{u}(t)$ and $f_{p}(t)$ as the time series of $C_{L}$ for the unperturbed case and that for the perturbed case, respectively. Then, a $L^{2}-$norm of the difference between these two data, $\displaystyle\int_{-T/2}^{T/2}\left(f_{u}(t+t_{0})-f_{p}(t)\right)^{2}dt$, was used to find the minimizer $t_{0}(\in[-T/2,T/2))$ such that the $L^{2}-$norm was minimized. Then, the phase difference between these two data, $\Delta\phi$, is estimated as $\displaystyle\Delta\phi=\left(\frac{2\pi}{T}\right)t_{0}$. For the comparison, $\tilde{Z}_{i,j}^{u}(\phi/\omega)$ and $\tilde{Z}_{i,j}^{v}(\phi/\omega)$ are used, which are defined as:

[TABLE]

Thus, the phase difference by the direct method, $\Delta\phi$, is related to $Z_{i,j}^{u}(\phi)$ by the formula $\tilde{Z}_{i,j}^{u}(\phi/\omega)=t_{0}/\epsilon_{0}$, and similar formula for $Z_{i,j}^{v}(\phi/\omega)$. The shifted data was generated by Fourier transformation and the phase shift of the Fourier coefficients, and the minimizer was obtained by the downhill simplex method. The one period length of data has been during $[9T,10T]$.

Figs. 8(a)-(d) show $\tilde{Z}_{i,j}^{u}$ and $\tilde{Z}_{i,j}^{v}$ obtained by the proposed method and direct method, and the direct method with $\epsilon_{0}=0.01$ and $0.02$. Figs. 8(a) and (c) show the values along the line segment indicated by red triangles in Fig. 7, whereas Figs. 8(b) and (d) show the values along the line segment indicated by blue points in Fig. 7. In all cases, the agreement of the values between the present method and the direct method, as well as the convergence of the direct method with different perturbation amplitudes, is reasonable. It is remarkable that the phase shift represented by the ratio with the period, $t_{0}/T$, is given $\epsilon\tilde{Z}_{i,j}^{u}/T$. Considering the magnitudes of $\tilde{Z}_{i,j}^{u}$ and $\tilde{Z}_{i,j}^{v}(\phi)$ are at most $0.15$, which indicates that $t_{0}/T\approx 0.003$ ($0.3\%$) or less if $\epsilon_{0}=0.01$ and $T\approx 5$. The value of the phase shift can be increased if the wider region is perturbed for a longer time interval according to the information of $\bm{Q}(\bm{x},\phi)$.

IV Summary

In this paper, a method to calculate the phase sensitivity function was developed, which is a fundamental function of the phase reduction theory. This method does not use the explicit form of the linearized matrix around the limit cycle (the Jacobian), which can be applied for the incompressible fluid system, where the Jacobian is determined by the Poisson equation. This method does not need the long-time integration until convergence like the direct method and the adjoint method, which reduces the computation time as well as the memory. Further, two applications were demonstrated: traveling pulse of the FitzHugh Nagumo equation in a periodic domain to validate this method, and the Kármán’s vortex street, to demonstrate the application to incompressible fluid systems.

The phase sensitivity function for the Kármán’s vortex street indicated how the phase shifts due to the external perturbation. Our analysis suggested that the phase shift property can be easily designed in the direct upstream region of the cylinder, where positive perturbations to the $x-$component of the velocity causes the phase advance, regardless of the phase in a wide area. This can be explained by a local speed-up due to the spread of the perturbative flow and constant Strouhal number for the Kármán’s vortex street. Higher values are obtained in the downstream area of the cylinder; however, the effective region is narrow and phase-dependent, which suggests that the control of the phase requires detailed design of the perturbation distribution and direction. This result is fundamental in controlling the phase of the Kármán’s vortex street.

The phase description is a powerful tool to analyze the phenomena with the limit cycle. For a large system, the numerical method to calculate the phase sensitivity function proposed here will be of great use, especially when the synchronization or the entrainment is considered. Further applications will be reported in future studies.

Acknowledgements.

This work was partially supported by the Mazda foundation. The author would like to thank Prof. Hiroya Nakao and Prof. Yoji Kawamura for discussions.