PDE Analysis of Pattern Formation in Biological Systems

• Divide the domain \(x \in [0, L]\) into \(M\) equally spaced points \(x_i\) with spacing \(\Delta x\).
• Field variables \(u_k(x, t)\) are sampled at node points: \(u^n_{k,i} \approx u_k(x_i, t_n)\).
• Derivatives are approximated using centered difference stencils around interior nodes.

• Second derivative (diffusion term) at interior: \[ \frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{\Delta x^2} \]
• Neumann BC at boundaries (zero-flux): \[ u_0' = \frac{u_1 - u_0}{\Delta x}, \quad u_M' = \frac{u_M - u_{M-1}}{\Delta x} \]
• Advection (chemotaxis) term: \[ -\chi \partial_x(u \partial_x v) \approx -\chi \cdot \frac{(u_{i+1} \delta v_{i+1}) - (u_{i} \delta v_i)}{\Delta x} \]

• The discretized equations at each node form a coupled ODE system: \[ \frac{du_{k,i}}{dt} = D_k \cdot \frac{u_{k,i+1} - 2u_{k,i} + u_{k,i-1}}{\Delta x^2} + \text{reaction}_{k,i} + \text{advection}_{k,i} \]
• Reactions are evaluated pointwise: \( f_{k,i} = f(u_i, v_i, w_i) \).

• Enforce zero-flux (Neumann) by symmetric ghost values: \[ u_0 = u_1, \quad u_M = u_{M-1} \]
• Applied in constructing finite difference stencils for boundary nodes.

• Temporal integration is performed using the implicit BDF (backward differentiation formula) method.
• Nonlinear terms solved with Newton iteration at each timestep.
• Recommended stability criterion: \[ \Delta t \leq \frac{\Delta x^2}{2 D_{\text{max}}} \]

# Second derivative with Neumann BCs
                                def second_derivative(u, dx):
                                    d2u = np.zeros_like(u)
                                    d2u[1:-1] = (u[2:] - 2*u[1:-1] + u[:-2]) / dx**2
                                    d2u[0] = (2*u[1] - 2*u[0]) / dx**2
                                    d2u[-1] = (2*u[-2] - 2*u[-1]) / dx**2
                                    return d2u
                            

• Partition the spatial domain \(x\in[0,L]\) into \(N\) control volumes (cells)
• Cell centers at \(x_j\) with width \(\Delta x\), faces at \(x_{j\pm1/2}\)
• Solution represented as cell averages: \[ \bar{u}^n_{k,j} \approx \frac{1}{\Delta x}\int_{x_{j-1/2}}^{x_{j+1/2}} u_k(x,t_n)dx \]

• Integrate PDE over each control volume and apply divergence theorem: \[ \frac{d\bar{u}_{k,j}}{dt} = \frac{1}{\Delta x}\left[F_{k,j+1/2} - F_{k,j-1/2}\right] + \bar{f}_{k,j} \]
• Diffusive flux (central difference): \[ F^{\text{diff}}_{k,j+1/2} = D_k\frac{u_{k,j+1} - u_{k,j}}{\Delta x} \]
• Advective flux (upwind scheme): \[ F^{\text{adv}}_{k,j+1/2} = \chi u_{k,j}\frac{c_{j+1} - c_j}{\Delta x} \]

• Resulting ODE system for each cell: \[ \begin{aligned} \frac{d\bar{u}_{k,j}}{dt} = &\frac{D_k}{\Delta x^2}(\bar{u}_{k,j+1} - 2\bar{u}_{k,j} + \bar{u}_{k,j-1}) \\ &- \frac{\chi}{\Delta x^2}\big[u_{k,j}(c_{j+1}-c_j) - u_{k,j-1}(c_j-c_{j-1})\big] + \bar{f}_{k,j} \end{aligned} \]

• Ghost cell values for no-flux boundaries: \[ u_{k,0} = u_{k,1}, \quad u_{k,N+1} = u_{k,N} \]
• Boundary fluxes: \[ F_{k,1/2} = 0, \quad F_{k,N+1/2} = 0 \]

• Time integration: Adaptive BDF2 method
• Nonlinear handling: Newton iterations for implicit steps
• Stability condition: \(\Delta t \leq \frac{\Delta x^2}{2D_{\text{max}}}\)