计算科学 - 随机微分方程系统 (SDE)：在双标量中遇到溢出 - 吾爱随笔录

我正在尝试集成 Dekker 等人的以下 SDE 系统。^[1]

{\begin{cases} \frac{d x}{d t} = a_{1} x^{3} + a_{2} x + ϕ + ζ_{x} \\ \frac{d y}{d t} = b_{1} z + b_{2} (κ (x) - (y^{2} + z^{2})) y + ζ_{y} \\ \frac{d z}{d t} = c_{1} y + b_{2} (κ (x) - (y^{2} + z^{2})) z + ζ_{z} \end{cases}

$\begin{cases} \frac{dx}{dt}=a_1x^3+a_2x+\phi+\zeta_x\\ \frac{dy}{dt}=b_1z+b_2(\kappa(x)-(y^2+z^2))y+\zeta_y\\ \frac{dz}{dt}=c_1y+b_2(\kappa(x)-(y^2+z^2))z+\zeta_z \end{cases}$

为此，我使用Euler-Maruyama 方法，该方法基本上是前向欧拉方案加上高斯噪声项，其中均值和标准差参数是从 Dekker 等人的表 2 中选取的。文章。

我们有

噪音平均值：0
噪声方差：0.1
积分时间：500
时间步长：0.5

参数 $a_1, a_2,...$ 以及耦合项 $\kappa(x)$ ，积分域、步长时间和初始条件均取自表 2。

但是，运行该模拟我得到一个RuntimeWarning: overflow encountered in double-scalars. 似乎步时间太大，但这是文章中推荐的。

我想知道为什么会出现此错误以及如何解决？

如您所见，经过几次迭代后，x、y 和 z 的小数点会呈指数级增长。它不是来自随机项，因为下面捕获中的打印（即数组 [x, y, z]）是仅来自前向欧拉的打印。

代码

# 1D EDO for Fold bifurcation paramaters
a1 = -1
a2 = 1

# 2D EDO for Hopf bifurcation parameters
b1 = b2 = 1
c1 = -1
c2 = 1

# coupling parameters
gamma1 = -0.1
gamma2 = 0.12

# initial conditions
x0 = -0.5
y0 = 1.0
z0 = -1.0

r0 = 5.0

# time range
t_init = 0
t_fin = 500
time_step = 0.01

# gaussian noise parameters
mean = 0.0
variance = 0.1

# parameters for stochastic plot
a1_stoch = -1
a2_stoch = 1
b1_stoch = 0.1
b2_stoch = 1
c1_stoch = -0.5
c2_stoch = 1
gamma1_stoch = -0.2
gamma2_stoch = 0.3
time_step_stoch = 0.5

import numpy as np
from consts import (a1_stoch, a2_stoch, b1_stoch, b2_stoch, c1_stoch, c2_stoch, gamma1_stoch, gamma2_stoch)


# linear coupling parameter
# proposed by Dekker et al. article
def gamma(x):
  return gamma1_stoch + (gamma2_stoch * x)

def fold_hopf_stoch(v, phi):
  print(v, phi)
  return np.array([
    a1_stoch * (v[0] ** 3) + a2_stoch * v[0] + phi,
    b1_stoch * v[2] + b2_stoch * (gamma(v[0]) - (v[1]**2 + v[2]**2))*v[1],
    c1_stoch * v[1] + c2_stoch * (gamma(v[0]) - (v[1]**2 + v[2]**2))*v[2]
  ])

import numpy as np
from consts import (mean, variance)

"""
forward euler method with a stochastical term
x_{i+1} = x_i * dt + zeta + sqrt(dt)
(forward-euler + gaussian noise * sqrt(dt))
"""
def forward_euler_maruyama(edo, v, dt, *args):
  zeta = np.random.normal(loc=mean, scale=np.sqrt(variance))
  print(zeta, edo(v, *args) * dt)
  return edo(v, *args) * dt + zeta * np.sqrt(dt)

import numpy as np
import matplotlib.pyplot as plt

from tqdm import tqdm

from .edo import (fold_hopf)
from .edo_stoch import (fold_hopf_stoch)
from .rk4 import rk4
from .euler import forward_euler_maruyama

from consts import (x0, y0, z0, t_init, t_fin, time_step, time_step_stoch)

plt.style.use("seaborn-whitegrid")
np.set_printoptions(precision=3, suppress=True)

class time_series():
  def __init__(self):
    # initial conditions
    self.initial_conditions = np.array([[x0, y0, z0]])
    # forcing parameter
    self.phi = -2
    # time
    self.t0 = t_init
    self.tN = t_fin
    # stochastic number of iterations
    self.niter = 10
    # legend
    self.legends = ["$x$ (leading)", "$y$ (following)", "$z$ (following)"]

  def solve(self, solver, edo, dt, nt):
    # [x0, y0, z0]
    #v = self.initial_conditions[0]

    # vector mesh -- will receive the result
    v_mesh = np.ones((nt, 3)) #
    # set inital conditions
    v_mesh[0] = self.initial_conditions[0]
    print(edo, dt)

    for t in tqdm(range(0, nt - 1)):
      v_mesh[t + 1] = v_mesh[t] + solver(edo, v_mesh[t], dt, self.phi)

      # increase forcing parameter
      self.phi += 0.001

    return v_mesh

  def basic(self):
    dt = time_step
    nt = int((self.tN - self.t0) / dt)
    time_mesh_basic = np.linspace(start=self.t0, stop=self.tN, num=nt)

    # [ [x0, y0, z0], [x1, y1, z1], ..., [xN, yN, zN] ]
    results = self.solve(rk4, fold_hopf, dt, nt)
    return time_mesh_basic, results

  def stochastic(self):
    dt = time_step_stoch
    nt = int((self.tN - self.t0) / dt)
    time_mesh_stoch = np.linspace(start=self.t0, stop=self.tN, num=nt)


    stochastic_results = np.ones((self.niter, nt, 3))
    # compute a lot of simulations
    for i in tqdm(range(0, self.niter)):
      results = self.solve(forward_euler_maruyama, fold_hopf_stoch, dt, nt)
      stochastic_results[i] = results
      return
    return time_mesh_stoch, stochastic_results

  def plot(self):
    #time_mesh_basic, basic_results = self.basic()
    time_mesh_stoch, stochastic_results = self.stochastic()
    return

    # 3 lines, 1 column
    fig, ((ax1), (ax2)) = plt.subplots(2, 1, figsize=(15, 7))
    fig.suptitle("Série temporelle")

    ax1.plot(time_mesh_basic, basic_results[:,0])
    ax1.plot(time_mesh_basic, basic_results[:,1])
    ax1.plot(time_mesh_basic, basic_results[:,2])

    ax1.set_xlabel("$t$")
    ax1.set_ylabel("$x$, $y$, $z$")
    ax1.set_xlim(0,500)
    ax1.set_ylim(-1.5,1.5)
    ax1.legend(self.legends, loc="upper right")
    ax1.set_title("Basique")

    print(stochastic_results)

    ax2.stackplot(time_mesh_stoch, stochastic_results[:,:,0])
    ax2.stackplot(time_mesh_stoch, stochastic_results[:,:,1])
    ax2.stackplot(time_mesh_stoch, stochastic_results[:,:,2])

    ax2.set_xlabel("$t$")
    ax2.set_ylabel("$x$, $y$, $z$")
    ax2.set_xlim(0,500)
    ax2.set_ylim(-1.5,1.5)
    ax2.legend(self.legends, loc="center left", bbox_to_anchor=(1,0.5))
    ax2.set_title("Stochastique")

    plt.tight_layout()
    plt.show()

参考：

德克尔，MM；冯德海特，AS；Dijkstra, HA 气候系统的级联转变。地球系统。达因。 2018, 9 (4), 1243–1260。DOI：10.5194/esd-9-1243-2018。