使用 odeint 的狄利克雷边界条件

Question

我正在尝试编辑 odeint 文档中的 Gray-Scott 一维方程示例（页面上的最后一个示例）。

我有下面的代码，它适用于诺依曼边界条件，但我想要左侧 x=0 处的狄利克雷边界条件，并在右侧 x=L 处保留诺依曼。我尝试将 y0 的维度减少 2，因为这样应该删除 2 个 ODE。但我只是用少 1 个内点重复同样的问题并得到相同的答案，这是有道理的，因为我没有任何地方可以修复左边界值。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import time
from scipy.integrate import odeint

def G(u, v, f, k):
    return f * (1 - u) - u*v**2

def H(u, v, f, k):
    return -(f + k) * v + u*v**2

def grayscott1d(y, t, f, k, Du, Dv, dx):
    """
    Differential equations for the 1-D Gray-Scott equations.

    The ODEs are derived using the method of lines.
    """
    # The vectors u and v are interleaved in y.  We define
    # views of u and v by slicing y.
    u = y[::2]
    v = y[1::2]

    # dydt is the return value of this function.
    dydt = np.empty_like(y)

    # Just like u and v are views of the interleaved vectors
    # in y, dudt and dvdt are views of the interleaved output
    # vectors in dydt.
    dudt = dydt[::2]
    dvdt = dydt[1::2]

    # Compute du/dt and dv/dt.  The end points and the interior points
    # are handled separately.
    dudt[0]    = G(u[0],    v[0],    f, k) + Du * (-2.0*u[0] + 2.0*u[1]) / dx**2
    dudt[1:-1] = G(u[1:-1], v[1:-1], f, k) + Du * np.diff(u,2) / dx**2
    dudt[-1]   = G(u[-1],   v[-1],   f, k) + Du * (- 2.0*u[-1] + 2.0*u[-2]) / dx**2
    dvdt[0]    = H(u[0],    v[0],    f, k) + Dv * (-2.0*v[0] + 2.0*v[1]) / dx**2
    dvdt[1:-1] = H(u[1:-1], v[1:-1], f, k) + Dv * np.diff(v,2) / dx**2
    dvdt[-1]   = H(u[-1],   v[-1],   f, k) + Dv * (-2.0*v[-1] + 2.0*v[-2]) / dx**2

    return dydt

rng = np.random.default_rng(0)
num_discretize_x = 100 #2500 gives their solution
y0 = rng.standard_normal(2*num_discretize_x)
y0 = np.ones(2*num_discretize_x)*3
t = np.linspace(0, 4, 100)
f = 0.024
k = 0.055
Du = 0.01
Dv = 0.005
dx = 0.01
x = np.linspace(0, 1, num_discretize_x)
t0 = time.time()
#sola = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx))
#t1 = time.time()
solb = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx), ml=2, mu=2)
print('solb shape: ', solb.shape)
t2 = time.time()
#print(f'No banding takes {t1-t0} s, while banding takes {t2-t1} s.')

u = solb[:,::2]
v = solb[:,1::2]

print(u.T)

t_grid, x_grid = np.meshgrid(t, x)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(t_grid, x_grid, u.T)
plt.xlabel('$t$')
plt.ylabel('$x$')
plt.show()
plt.close('all')


plt.plot(t, u[:,-1], label='$u(x=L)$')
plt.plot(t, v[:,-1], label='$v(x=L)$')
plt.xlabel('$t$')
plt.ylabel('$u$ or $v$')
plt.legend(loc=0)
plt.show()

Answer 1

您可以将 u 和 v 的值固定在 x = 0 处，然后将诺依曼条件保持在 x = L 处：

代码：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import time
from scipy.integrate import odeint


def G(u, v, f, k):
    return f * (1 - u) - u * v**2


def H(u, v, f, k):
    return -(f + k) * v + u * v**2


def grayscott1d(y, t, f, k, Du, Dv, dx, u0_fixed=1.0, v0_fixed=0.5):
    u, v = y[::2], y[1::2]
    dydt = np.empty_like(y)
    dudt, dvdt = dydt[::2], dydt[1::2]
    dudt[0], dvdt[0] = 0, 0
    u[0], v[0] = u0_fixed, v0_fixed

    dudt[1:-1] = G(u[1:-1], v[1:-1], f, k) + Du * np.diff(u, 2) / dx**2
    dudt[-1] = G(u[-1], v[-1], f, k) + Du * (-2.0 * u[-1] + 2.0 * u[-2]) / dx**2
    dvdt[1:-1] = H(u[1:-1], v[1:-1], f, k) + Dv * np.diff(v, 2) / dx**2
    dvdt[-1] = H(u[-1], v[-1], f, k) + Dv * (-2.0 * v[-1] + 2.0 * v[-2]) / dx**2

    return dydt


rng = np.random.default_rng(0)
num_discretize_x = 100
y0 = rng.standard_normal(2 * num_discretize_x)
y0 = np.ones(2 * num_discretize_x) * 3
t = np.linspace(0, 4, 100)
f = 0.024
k = 0.055
Du = 0.01
Dv = 0.0005
dx = 0.01
x = np.linspace(0, 1, num_discretize_x)

t0 = time.time()
solb = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx), ml=2, mu=2)
print('solb shape: ', solb.shape)
t2 = time.time()

u = solb[:, ::2]
v = solb[:, 1::2]

print(u.T)

t_grid, x_grid = np.meshgrid(t, x)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(t_grid, x_grid, u.T)
plt.xlabel('$t$')
plt.ylabel('$x$')
plt.show()
plt.close('all')

plt.plot(t, u[:, -1], label='$u(x=L)$')
plt.plot(t, v[:, -1], label='$v(x=L)$')
plt.xlabel('$t$')
plt.ylabel('$u$ or $v$')
plt.legend(loc=0)
plt.show()

打印

solb shape:  (100, 200)
[[3.00000000e+00 1.00000000e+00 1.00000000e+00 ... 1.00000000e+00
  1.00000000e+00 1.00000000e+00]
 [3.00000000e+00 1.17379393e+00 8.52324657e-01 ... 3.41682696e-01
  3.41131018e-01 3.40594796e-01]
 [3.00000000e+00 1.34513346e+00 7.43775028e-01 ... 9.71922350e-02
  9.68656847e-02 9.65485574e-02]
 ...
 [3.00000000e+00 1.77119093e+00 6.99074218e-01 ... 1.18100847e-03
  1.18800490e-03 1.19504096e-03]
 [3.00000000e+00 1.77119093e+00 6.99074218e-01 ... 1.18100847e-03
  1.18800490e-03 1.19504096e-03]
 [3.00000000e+00 1.77119093e+00 6.99074218e-01 ... 1.18100847e-03
  1.18800490e-03 1.19504096e-03]]