Probing the evolutionary history of epigenetic mechanisms: what can we learn from marine diatoms

1.   Introduction

Waveform relaxation (WR) is an iterative method for systems of ordinary differential equations (ODEs). It is introduced for the first time by Lelarasmee et al. ^[1] for the time domain analysis of large-scale nonlinear dynamical systems. A lot of studies have been done on WR methods and most of them focus on convergence (see ^{[2,3,4,5,6,7,8,9]}). Zero-stability is the basic property of numerical methods of ODEs^[10]. However, to our best knowledge, so far there exists little work on zero-stability of WR methods.

We first of all propose two examples to show the fact that tiny perturbation may lead to huge change of WR approximation solutions, at the same time, to confirm that it is necessary to study zero-stability of WR methods. For this, it is enough to consider a simple method, and suppose only the initial value has tiny perturbations.

Example 1. Consider the WR method

of solving numerically the initial value problem

which is gotten by applying the two step method

that is consistent for any $ \alpha $ and zero-stable for $ -1 < \alpha < 1 $ to the iterative scheme

Suppose that

and

We can get the approximate solution $ \{\{x^{(k)}_n, n = 0, 1, \cdots, N\}, k = 0, 1, \cdots\} $ by using Eqs (1.1)–(1.4). Suppose that tiny perturbations denoted by $ \delta $ lead to the initial values become

and the resulting perturbed solution of Eq (1.1) is denoted by $ \{\tilde{x}^{(k)}_n\} $. Let $ \delta = 1 $ for simplicity and denote $ \tilde{x}^{(k)}_n-x^{(k)}_n $ with $ e^{(k)}_n $. Then we have

By virtue of Eq (1.6) the following numerical results are obtained.

The data in Tables 1 and 2 show that errors of the approximate solution of Eq (1.1) brought by perturbations of initial values are controllable for $ \alpha = 0.5 $, but they are unbounded for $ \alpha = 2 $ when $ h $ tends to zero.

In this numerical examples, taking $ \delta = 1 $ instead of tiny perturbations seem to be unreasonable. The main reason for this is that we are usually interested in the ratio $ \frac{e^{(k)}_N}{\delta} $ instead of $ \delta $ itself. In fact we have also performed some other experiments using smaller $ \delta $, in which the similar ratios are obtained for different $ \delta $.

Example 2. Consider the one-step WR method

where $ 2/h $ is an integer. Suppose that the initial value $ 0 $ becomes $ \delta $ by the influence of perturbations and the solution $ x_n^{(k)} $ becomes $ \tilde{x}_n^{(k)} $ by the influence of $ \delta $, which satisfies

Clearly the numerical solution generated by Eq (1.7) is the constant zero and Eq (1.8) can be regarded as the WR method of solving the following initial value problem

It is easy to derive that $ t_0 = \dfrac{1}{\delta^{0.001}} $ is the blow-up time of the solution of Eq (1.9) and $ t_0\in [0, 2] $ whenever $ \delta > 10^{-301} $. Consequently, $ \sup\limits_{0 < t < 2}x(t) = +\infty $ for any $ \delta > 10^{-301} $. If the "2" in Eqs (1.8) and (1.9) is replaced by $ t_0-\varepsilon $, where $ \varepsilon > 0 $ is small enough, then Eq (1.8) converges to Eq (1.9), that is,

see ^[4]. Thus

that is, the change of the solution of Eq (1.8) brought by the tiny change of initial values is unbounded as $ h $ tends to zero.

The above two examples show that tiny perturbations can lead to huge change of WR approximation solutions, although the methods used in the examples are too simple to be practically useful.

A numerical method is said to be zero-stable if the change of its solution brought by the tiny perturbations is controllable when $ h $ tends to zero (see ^[10]). In this paper we will explore what conditions guarantee the zero-stability of WR methods.

2.   Preliminary

Consider the initial problem

Taking the splitting function $ F(t, x, x) = f(t, x) $ we can construct the iterative scheme

Using one-step methods with variable step-size to discrete Eq (2.1) one arrives at

where $ n = 0, 1, \cdots, N-1, k = 0, 1, \cdots $, $ \sum\limits_{n = 0}^{N-1}h_n = T $ and $ x^{(0)}_n = x_0 $ for all $ n $. Applying a class of multi-step methods with fixed step-size to Eq (2.1) we get

where $ n = p, p+1, \cdots, N-1, k = 0, 1, \cdots $, $ h = T/N $, the initial values except $ x^{(0)}_n(n = 0, 1, \cdots, N) $ and $ x^{(k)}_0(k = 0, 1, \cdots) $ are generated by a suitable one-step method. Here we take the initial values as

A classical example of Eq (2.3) is Admas-type linear multi-step methods:

where $ \beta_j, j = 0, 1, \cdots, p+1 $ are constants. The linear multi-step methods use generally the fixed step-size and therefore we only consider the Eq (2.3) with the fixed step-size.

Let $ _ix^{(k)} $ denote the component of $ x^{(k)} $ satisfying

where $ _ix^{(k)}\in \mathbb{R}^{d_i} $, $ x^{(k)}\in \mathbb{R}^{d} $ and $ d_1+d_2+\cdots+d_m = d $. The Eq (2.1) can be rewritten as

with $ i = 1, 2, \cdots, m $, that is, the large Eq (2.1) is divided into $ m $ subsystems. When taking

one arrives at Guass-Jacobi WR method used frequently in actual computation, which consists of $ m $ independent subsystems and is hence parallel in nature(see ^[2,6,9]). In Eq (2.2) the unique mesh $ 0 = t_0 < t_1 < \cdots < t_{N} = T $ is applied to all subsystems. However the subsystems of Gauss-Jocobi WR methods are independent and may have distinct behaviors, and consequently it is better to use different meshes for the subsystems with different behaviors. Applying the mesh

to $ i $th subsystem yields the following multi-rate WR method

where $ i = 1, 2, \cdots, m, k = 0, 1, \cdots $. Here the initial values $ \sideset{_i}{^{(k+1)}_0}{\mathop{x}} = \sideset{_i}{^{(0)}_n}{\mathop{x}} = \sideset{_i}{_0}{\mathop{x}} $ for all $ k $ and $ n $, $ \left(\sideset{_1}{_0}{\mathop{x}}^{\rm{T}}, \sideset{_2}{_0}{\mathop{x}}^{\rm{T}}, \cdots, \sideset{_m}{_0}{\mathop{x}}^{\rm{T}} \right)^{\rm{T}} = x_0 $ and $ \sideset{_j}{^{(k)}}{\mathop{y}}(t) $ is the interpolation function satisfying $ \sideset{_j}{^{(k)}}{\mathop{y}}(\sideset{_j}{_n}{\mathop{t}}) = \sideset{_j}{^{(k)}_n}{\mathop{x}} $ for $ n = 0, 1, \cdots, N_j $.

Because of the presence of errors in actual computing, it is necessary to consider the following perturbed systems of Eqs (2.2), (2.3) and (2.5) generated by the perturbations $ \{\delta^{(k)}_n, n = 0, 1, \cdots, N, k = 0, 1, \cdots\} $

and

where $ i = 1, 2, \cdots, m $ and $ \sideset{_j}{^{(k)}}{\mathop{\tilde{y}}}(t) $ is the interpolation function satisfying $ \sideset{_j}{^{(k)}}{\mathop{\tilde{y}}}(\sideset{_j}{_n}{\mathop{t}}) = \sideset{_j}{^{(k)}_n}{\mathop{\tilde{x}}} $ for $ n = 0, 1, \cdots, N_j, j = 0, 1, \cdots, m $.

Definition 2.1. (Page 32 of ^[10]) Let $ \{\delta^{(k)}_n, n = 0, 1, \cdots, N, k = 0, 1, \cdots\} $ be any perturbation of Eq (2.2) or Eq (2.3) and Eq (2.4), and let $ \{\tilde{x}^{(k)}_n, n = 0, 1, \cdots, N, k = 0, 1, \cdots\} $ be the solution of the resulting perturbed system Eq (2.6) or Eq (2.7). Then if there exist constants $ C $ and $ h_0 $ such that, for all $ h\in(0, h_0] $($ h = \max\limits_n h_n $ for Eq (2.2))

whenever

we say that the Eq (2.2) or Eqs (2.3) and Eq (2.4) is zero-stable.

Here and hereafter $ \|\cdot\| $ denotes any norm defined in $ \mathbb{R}^d $.

Definition 2.2. Let $ \{_i\delta^{(k)}_n, i = 1, 2, \cdots, m, n = 0, 1, \cdots, N, k = 0, 1, \cdots\} $ be any perturbation of Eq (2.5), and let $ \{_i\tilde{x}^{(k)}_n, i = 1, 2, \cdots, m, n = 0, 1, \cdots, N, k = 0, 1, \cdots\} $ be the solution of the resulting perturbed Eq (2.8). Then if there exist constants $ h_0 $ and $ C_k $ depended on $ k $ such that

for all $ h = \max\limits_{

$$\begin{subarray}{c}1\leq i\leq m\end{subarray}$$ } \max\limits_{0\leq n\leq N_i}\sideset{_i}{_n}{\mathop{h}}\in(0, h_0] $, whenever

we say that the Eq (2.5) is weakly zero-stable.

The following several lemmas are useful for studying zero-stability of the WR methods mentioned above.

Lemma 2.3. Let $ a, b $ and $ T $ be nonnegative constants, and $ N $ be a positive integer. Suppose that the sets $ \{h_n > 0, n = 0, 1, \cdots, N-1\} $ and $ \{u_n > 0, n = 0, 1, \cdots, N\} $ satisfy

Then

provided that $ \sum\limits_{n = 0}^{N-1}h_n = T $.

Proof. By using Eq (2.9) repeatedly we have

where $ n = 0, 1, \cdots, N-1 $. Let $ t_0 = 0, t_{n+1} = t_n+h_n, n = 0, 1, \cdots, N-1 $, then $ t_N = T $. Noting that $ 1+x < {\mathrm{e}}^x $ for $ x > 0 $ we can derive easily from above inequality

for $ n = 0, 1, \cdots, N-1 $. This together with monotonicity of the function $ {\mathrm{e}}^{-at} $ imply

for all $ n = 0, 1, \cdots, N-1 $.

Lemma 2.4. Let $ a, b, c, d $ and $ T $ be nonnegative constants, $ N $ be a positive integer and $ h_n, n = 0, 1, \cdots, N-1, $ be positive real numbers.

Suppose that the sequence of positive numbers $ \{e_{n}^{(k)} > 0, n = 0, 1, \cdots N, k = 0, 1, \cdots\} $ satisfies

Then for any $ 0 < h = \max\limits_{0\leq n\leq N-1} h_n < \dfrac{1}{2c} $

provided that $ \sum\limits_{n = 0}^{N-1}h_n = T $.

Proof. By Eq (2.10) we have

for $ n = 0, 1, \cdots, N-1 $. Let $ \varepsilon = \max\limits_{0\leq n\leq N}e^{(0)}_n $ and $ u_n = \max\left\{\max\limits_{0\leq k}e^{(k)}_n, \varepsilon\right\}, n = 0, 1, \cdots, N $. Clearly $ \max\left\{\max\limits_{0\leq k}e^{(k+1)}_n, \varepsilon\right\} = u_n $. Hence we can get by Eq (2.12)

Suppose that $ 0 < h < \dfrac{1}{2c} $. Consequently, $ \dfrac{1}{1-ch_n} < 2 $ and the above inequality thus yield

This together with Lemma 2.3 yield

Thus Eq (2.11) holds true. The proof is complete.

Lemma 2.5. Let $ a_j, b_j(j = 0, 1, \cdots, p), c, d $ and $ T $ be nonnegative constants. Let $ N $ be a positive integer and $ h = T/N $.Suppose that the sequence of positive numbers $ \{e_{n}^{(k)} > 0, n = 0, 1, \cdots N, k = 0, 1, \cdots\} $ satisfies

Then for any $ 0 < h < \dfrac{1}{2c} $

where $ \varepsilon = \max\left\{ \max\limits_{0\leq n\leq N}e^{(0)}_n, \max\limits_{

$$\begin{subarray}{c} 0\leq k, 0\leq n\leq p \end{subarray}$$ }e^{(k)}_n\right\} $, $ A = \sum\limits_{j = 0}^p (a_j+b_j) $.

Proof. Note that $ e^{(k)}_{n}\leq \varepsilon $ for all $ n\leq p $ and all $ k\geq 0 $. By using Eq (2.13) repeatedly we have

Let $ 0 < h < \dfrac{1}{2c} $. With the notation $ \alpha_h = \dfrac{1+Ah}{1-ch} $ and $ \beta_h = \dfrac{dh}{1-ch} $ we can derive from Eq (2.14)

Noting that $ \varepsilon\leq \alpha_h\varepsilon+\beta_h $ and therefore $ e^{(k)}_{n}\leq \alpha_h\varepsilon+\beta_h $ for all $ n\leq p+1 $ and all $ k\geq 0 $ we thus have by repeating the above process

Consequently, we get

hold true for $ m = 1, 2, \cdots, N-p $. Note that $ \alpha_h\geq 1 $. Thus

for $ k = 0, 1, \cdots $, $ n = 0, 1, \cdots, N $. This together with $ 0 < ch < 1/2 $ yield

for $ k = 0, 1, \cdots, n = 0, 1, \cdots, N $. The proof is complete.

3.   Main results

Theorem 3.1. Suppose that there exist constants $ L_1 $, $ L_2 $, $ L_3 $ and $ L_4 $ such that

for any real numbers $ t, h\in \mathbb{R}, x_i, y_i\in \mathbb{R}^{d}, i = 1, 2, 3, 4 $. Then the WR Eq (2.2) is zero-stable.

Proof. Let $ \varepsilon^{(k)}_n = \|\tilde{x}^{(k)}_n-x^{(k)}_n\| $ and for any $ \varepsilon > 0 $ the perturbations in Eq (2.6) $ \{\delta^{(k)}_n, k = 0, 1, \cdots, n = 0, 1, \cdots, N\} $ satisfies $ \max\limits_{

$$\begin{subarray}{c}0\leq k, 0\leq n\leq N\end{subarray}$$ }\|\delta^{(k)}_n\| < \varepsilon $. Then by virtue of Eqs (2.2) and (2.6) and the condition (3.1) we have

for $ n = 0, 1, \cdots, N-1, k = 0, 1, \cdots $ and

Let $ h_n < \dfrac{1}{2L_1} $. Then $ \dfrac{1}{1-h_nL_1} < 2 $. We therefore derive from Eq (3.2)

for $ n = 0, 1, \cdots, N-1, k = 0, 1, \cdots $. Consequently, by Eqa (3.3), (3.4) and Lemma 2.4 we have

for any $ k = 0, 1, \cdots, n = 0, 1, \cdots, N $ if $ 0 < h < \min\left\{\dfrac{1}{2L_1}, \dfrac{1}{4L_3}\right\} $. The proof is therefore complete.

Theorem 3.2. Suppose that there exist constants $ L_i, i = 1, 2, \cdots, 2p+4 $ such that

for any real numbers $ t, h\in \mathbb{R}, x_i, y_i\in \mathbb{R}^d, i = 1, 2, \cdots, 2p+4 $. Then the WR Eq (2.3) and Eq (2.4) is zero-stable.

Proof. Let $ \varepsilon^{(k)}_n $ denote $ \|\tilde{x}^{(k)}_n-x^{(k)}_n\| $ and the perturbations $ \{\delta^{(k)}_n, k = 0, 1, \cdots, n = 0, 1, \cdots, N\} $ in Eq (2.7) satisfy $ \max\limits_{

$$\begin{subarray}{c}0\leq k, 0\leq n\leq N\end{subarray}$$ }\|\delta^{(k)}_n\| < \varepsilon $, where $ \varepsilon $ is any positive number. Then by virtue of Eqs (2.3), (2.4), (2.7) and the condition (3.5) we have

where $ k = 0, 1, \cdots $. Let $ h < \dfrac{1}{2L_1} $. Consequently, $ \dfrac{1}{1-hL_1} < 2 $. We therefore derive from Eq (3.6)

for $ n = 0, 1, \cdots, N-1, k = 0, 1, \cdots $. By Eq (3.7) and Lemma 2.5 we have

for any $ k = 0, 1, \cdots, n = 0, 1, \cdots, N $ if $ 0 < h < \min\left\{\dfrac{1}{2L_1}, \dfrac{1}{4L_{p+3}}\right\} $. The proof is complete.

Theorem 3.3. For Eqs (2.5) and (2.8) let $ \Phi = (\Phi_1^{\mathrm{T}}, \Phi_2^{\mathrm{T}}, \cdots, \Phi_m^{\mathrm{T}})^{\mathrm{T}} $ satisfy the Lipschitz condition:there exist constants $ L_1 $ and $ L_2 $ such that

and the interpolation functions $ \sideset{_j}{^{(k)}}{\mathop{\tilde{y}}}(t) $ and $ \sideset{_j}{^{(k)}}{\mathop{y}}(t) $ satisfy that

where $ C $ is the positive constant. Then the multi-rate Eq (2.5) is weakly zero-stable.

Proof. Define $ \|x\|_a = \max\limits_{1\leq i\leq m}\|_ix\| $, $ \|x\|_b = \sum\limits_{i = 1}^m\|_ix\| $ for $ x = (\sideset{_1}{^{\rm{T}}}{\mathop{x}}, \sideset{_2}{^{\rm{T}}}{\mathop{x}}, \cdots, \sideset{_m}{^{\rm{T}}}{\mathop{x}})^{\rm{T}}\in \mathbb{R}^d $, where $ _ix\in \mathbb{R}^{d_i}, i = 1, 2, \cdots, m, d_1+d_2+\cdots+d_m = d $. It is easy to prove that $ \|\cdot\|_a $ and $ \|\cdot\|_b $ are the norm. Consequently, there exist constants $ C_1 $ and $ C_2 $ such that

Eqs (3.8) and (3.10) thus yield

for any $ x_1, x_2, y_1, y_2\in \mathbb{R}^d $ and therefore

With the notation $ _ie^{(k)}_n = \|_i\tilde{x}^{(k)}_n-\sideset{_i}{^{(k)}_n}{\mathop{x}}\| $ we can derive from Eqs (2.5), (2.8), (3.9) and (3.11)

Let $ h = \max\limits_{

$$\begin{subarray}{c}1\leq i\leq m \end{subarray}$$ }\max\limits_{0\leq n\leq N_i} \sideset{_i}{_n}{\mathop{h}} < \dfrac{C_1}{2L_1C_2C} $. Consequently, $ _ih_n\dfrac{L_1C_2C}{C_1} < \dfrac12 < 1-\sideset{_i}{_n}{\mathop{h}}\dfrac{L_1C_2C}{C_1} $ which together with Eq (3.12) imply that

with $ a = \dfrac{2(L_1+L_2)C_2C}{C_1} $. By using Eq (3.13) with $ k = 0 $ and Lemma 2.3 we get

Thus

with $ a_1 = {\mathrm{e}}^{aT}+\dfrac{1}{a}({\mathrm{e}}^{aT}-1)(a+2) $. Similarly, using the above inequality, Eq (3.13) with $ k = 1 $ and Lemma 2.3 we have

with $ a_2 = {\mathrm{e}}^{aT}+({\mathrm{e}}^{aT}-1)(ma_1+\dfrac{2}{a}) $. Repeating the above process yields that for any $ k $ there exists the constant $ a_k $ such that

for any mesh $ \{_it_n, n = 0, 1, \cdots, N, i = 1, 2, \cdots, m\} $. The proof is therefore complete.

Remark 1. In the above discussions we assume that the perturbation at one step is $ \delta_n^{(k)} h $, which depends on step-size $ h $. This assumption is necessary for the proof of the theorems mentioned above. However it seems to be more reasonable to regard $ \delta_n^{(k)} $, a small number, as the perturbation, if the step-size is very small. In this case the tiny perturbation at each step may bring out the huge change of the numerical solution. We will explore these in next section by some computer experiments.

Remark 2. Note that the convergence property of the WR methods is not be used in the proof of Theorem 3.1-3.3. Hence, a divergent WR method may be zero-stable. In other words, zero-stability does not imply convergence.

4.   Numerical experiments

In this section some numerical experiments are performed to verify the theorems obtained in Section 3.

Consider the equation

where $ x\in \mathbb{R}^{15} $, $ T_i = \left(

$$\begin{array}{cccc} -2^i & 1 & & \\ 1 & -2^i & \ddots & \\ & \ddots & \ddots & 1\\ & & 1 & -2^i\\ \end{array}$$ \right)\in \mathbb{R}^{5\times 5} $ and $ I $ is the identity matrix. The special case of such system has been examined by Burrage ^[11], which is obtained by discretizing the heat equation in two variables.

Let $ A_1 = \left(

$$\begin{array}{ccc} T_1 & & \\ & T_2 & \\ & & T_3 \\ \end{array}$$ \right) $ and $ A_2 = A-A_1 $. First, choose the following WR methods of solving the Eq (4.1):

(i) Euler method with variable step-size

(ii) Two step methods with fixed step-size

where $ x^{(k+1)}_1 $ is obtained by a suitable one-step method. Here let

(iii) Multi-rate Jacobi methods based on Euler method

Second, assume the initial values and the approximation solution at each calculation step have tiny perturbations as in Eqs (2.6)–(2.8) and let $ \delta^{(k)}_n\equiv \Delta = \delta\cdot (1, 1, \cdots, 1)^{\rm{T}}\in \mathbb{R}^{15} $, $ _i\delta^{(k)}_n\equiv \Delta_p = \delta\cdot (1, 1, 1, 1, 1)^{\rm{T}}\in \mathbb{R}^5 $, $ \tilde{x}^{(k)}_n $ denote the resulting perturbed solution and $ u^{(k)}_n $ denote $ \tilde{x}^{(k)}_n-x^{(k)}_n $, which yields

and

Lastly, reach a reasonable conclusion by analysing data generated by Eqs (4.2)–(4.4).

We choose $ A $, $ A_1 $, $ A_2 $ and $ \delta^{(k)}_n $ as above to guarantee that $ u^{(k)}_n $ increase with $ k $ and $ n $ that is the worst case of error propagation.

In our experiments the step-sizes $ h_n $ in Eq (4.2) and $ _ih_n $ in Eq (4.4) are taken as

and

where $ h > 0 $, $ \{\xi_0, \xi_1, \cdots\} $ is a sequence of independent and identically distributed random variables satisfying $ \xi_i\sim N(0, 1) $.

The data in Tables 3, 5 and 7 show clearly the errors resulted by the perturbations $ h_n\Delta(h\Delta/2) $ in Eq (4.2) (Eq (4.3) are controllable as $ k $ and $ n $ tend to infinity which support the theorems developed in Section 3). However the data in Table 4 and 6 show that the errors increase with $ k $ and $ 1/h $ when the perturbations $ h_n\Delta (h\Delta/2) $ are replaced by the constant perturbation $ \Delta(\Delta/2) $.

5.   Further work

It is well known that a linear multi-step method with fixed step-size is convergent if, and only if, it is both consistent and zero-stable. In this paper we only present the sufficient conditions of zero-stability for some special methods. We will explore the relationship between convergence and zero-stability for WR methods based on general linear multi-step methods in the future. Moreover, we shall apply our methods to the numerical study of some inverse problems ^[12,13,14].

Acknowledgments

Research is supported by the Natural Science Foundation of Fujian Province, China (grant number 2021J011031).

Conflict of interest

The author declares that he has no conflicts of interest regarding the publication of this paper.