Search code examples
pythonscipyodedifferential-equationspde

Dirichlet boundary conditions using odeint

I am trying to edit the Gray-Scott 1D equation example (last example on the page) in the odeint documentation.

I have the code below and it works for the Neumann boundary conditions but I want a Dirichlet boundary condition on the left at x=0 and retain Neumann on the right at x=L. I tried to reduce the dimension of y0 by 2 as that should remove 2 ODEs. But I was just repeating the same problem with 1 less interior point and getting the same answer, which makes sense since I didn't have any place where I fixed the left boundary values.

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()
Solution

  • You can fix the value of u and v at x = 0 and then keep the Neumann conditions at x = L:

    Code:

    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()

    Prints

    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]]

    enter image description here

    enter image description here

    Comment:

    • @Warren Weckesser: Instead of using the parameters u0_fixed and v0_fixed and setting u[0] and v[0] in each call of grayscott1d(), you can set the corresponding initial conditions to those values (i.e. y0[:2] = [u0_fixed, v0_fixed]. Then, since dudt[0] = 0 and dvdt[0] = 0, the solution will remain constant at those values.

    Code:

    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):
    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
    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
y0[0], y0[1] = 1.0, 0.5

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)

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

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()