Deep FBSDE Neural Networks for Solving Incompressible Navier-Stokes Equation and Cahn-Hilliard Equation

The FBSDEs where the randomness in the BSDE driven by a forward stochastic differential equation (SDE), is written in the general form

where $\{W_{s}\}_{0\leq s\leq T}$ is a d-dimensional Brownian motion, $b:[0,T]\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$, $\sigma:[0,T]\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times d}$, $f:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{m}\times\mathbb{R}^{d\times m}\rightarrow\mathbb{R}^{m}$ and $g:[0,T]\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{m}$ are all deterministic mappings of time and space, with the fixed $T>0$. We refer to $Z$ as the $control\ process$ according to the stochastic control terminology. In order to guarantee the existence of a unique solution pair $\{(Y_{s},Z_{s})\}_{0\leq s\leq T}$ adapted to the augmented natural filtration, the standard well-posedness assumptions of [23] are required. Indeed, considering the quasi-linear, parabolic terminal problem

with $u(T,x)=g(x)$ and $\mathcal{L}$ is the second-order differential operator

the nonlinear Feynman-Kac formula indicates that the solution of (1) coincides almost exactly with the solution of (2) (cf., e.g., [23])

As a result, the BSDE formulation offers a stochastic representation to the synchronous solution of a parabolic problem and its gradient, which is a distinct advantage for numerous applications in stochastic control.

Inorder to review the deep BSDE method proposed in [9], we consider the following FBSDEs

which is the integration form of (1). Given a partition of the time interval $[0,T]:0=t_{0}<t_{1}<...<t_{N}=T$, the Euler-Maruyama scheme is used to discretize for $X_{t}$ and $Y_{t}$ and we have

where $\Delta t_{n}=t_{n+1}-t_{n}=\frac{T}{N}$, $\Delta W_{t_{n}}=W_{t_{n+1}}-W_{t_{n}}$. The $Z_{t_{n}}(X_{t_{n}})$ is approximated by a FNN with parameter $\theta_{n}$ for $n=1,\cdots,N-1$. The initial values $Y_{t_{0}}(X_{t_{0}})$ and $Z_{t_{0}}(X_{t_{0}})$ are treated as trainable parameters in the model. To make $Y_{t_{0}}(X_{t_{0}})$ to approximate $u(t_{0},X_{t_{0}})$, the difference in the matching with a given terminal condition is used to define the expected loss function

which represents $M$ different realizations of the underlying Brownian motion, where the subscript $m$ corresponds to the $m$-th realization of the underlying Brownian motion. The process is called the deep BSDE method.

Raissi [6] introduced neural networks called FBSNNs to solve FBSDEs. The unknown solution $u(t,x)$ is approximated by the FNN with the input $(t,x)$ and the required gradient vector $\nabla u(t,x)$ is attained by applying AD. The parameter $\theta$ of FNN can be learned by minimizing the loss function given explicitly in equation (9) obtained from discretizing the FBSDEs (5) using the Euler-Maruyama scheme

where $Y_{t_{n}}(X_{t_{n}})$ represents the estimated value of $u(t_{n},X_{t_{n}})$ given by the FNN and $\widetilde{Y}_{t_{n+1}}(X_{t_{n+1}})$ is the reference value corresponding to $Y_{t_{n+1}}(X_{t_{n+1}})$, which is obtained from the calculation in (8). The loss function is then given by

where the subscripts $m$ is the same meaning as it in (7).

The deep BSDE method only calculates the value of $u(t_{0},X_{t_{0}})$. This means that in order to obtain an approximation to $Y_{t_{n}}(X_{t_{n}})=u(t_{n},X_{t_{n}})$ at a later time $t_{n}>t_{0}$, we have to retrain the algorithm. Furthermore, the number of the FNNs grows with the number of time steps $N$, which makes training difficult. In this article, we use the FBSNNs. The FNN is expected to be able to approximate $u(t,x)$ over the entire computational area instead of only one point. That is, we will use multiple initial points to train the FNN. In order to improve the efficiency, the number of Brownian motion trajectories for each initial point is set as $M=1$.

The Cauchy problem for deterministic backward Navier–Stokes equation for the velocity field of the incompressible and viscous fluid is

which is obtained from the classical Navier–Stokes equation via the time-reversing transformation

Here $u=u(t,x)$ represents the $d$-dimensional velocity field of a fluid, $p=p(t,x)$ is the pressure, $\nu>0$ is the viscosity coefficient, and $f=f(t,x)$ stands for the external force. We now study the backward Navier–Stokes equation (10) in $\mathbb{R}^{d}$ with different boundary conditions.

Then, the PDE (10) is associated through the nonlinear Feynman-Kac formula to the following FBSDEs

where

Given a partition of $[0,T]:0=t_{0}<t_{1}<...<t_{N}=T$, we consider the Euler-Maruyama scheme with $n=0,...,N-1$ for FBSDEs (12)

where $\Delta t_{n}=t_{n+1}-t_{n}=\frac{T}{N}$ and $\Delta W_{t_{n}}=W_{t_{n+1}}-W_{t_{n}}$. The $(Y_{t_{n}},P_{t_{n}})^{T}$ represents the estimated value of $(u,p)^{T}$ at time $t_{n}$ given by the FNN, respectively. The $\widetilde{Y}_{t_{n+1}}(X_{t_{n+1}})$ is the reference value of $Y_{t_{n+1}}(X_{t_{n+1}})$, which is obtained from the calculation in the second equation in (13). We utilize $K$ different initial sampling points for training the FNN. The algorithm of the proposed scheme is summarized in Algorithm 1. Illustration of the Algorithm 1 for solving the incompressible Navier-Stokes equation is shown in Figure 1.

For the backward Navier–Stokes equation (10) in $\Omega\subset\mathbb{R}^{d}$ with the Dirichlet boundary condition

the corresponding FBSDEs can be rewritten as the following form according to [43] by the nonlinear Feynman-Kac formula

where $a\Lambda b=\min\{a,b\}$, the stopping time $\tau=\inf\{t>0:X_{t}\notin\Omega\}$ be the first time that the process $X_{t}$ exits $\Omega$ and

Through the Euler scheme, the discrete formulation of the FBSDEs (16) can be obtained accordingly

where $\mathbbm{1}_{(0,\tau)}(t_{n})=1,t_{n}\in[0,\tau)$. It should be noted that we will calculate the stop time $\tau$ after the iteration of $X_{t_{n+1}}$ is completed. Supposing $\tau=t_{n+1}$ when $\tau\leq T$, we let $X_{\tau}=X_{t_{n+1}}$ on $\partial\Omega$ and update $\Delta W_{t_{n}}=(X_{t_{n+1}}-X_{t_{n}})/\sqrt{2\nu}$ to satisfy (18). The algorithm of the proposed scheme is similar as Algorithm 1.

We consider the backward Navier–Stokes equation (10) with the Neumann boundary condition

where $\mathbf{n}$ is the unit normal vector at $\partial\Omega$ pointing outward of $\Omega$. Supposing $x\in\partial\Omega,x+\Delta x\notin\Omega$ and $x-\Delta x\in\Omega$, where $x+\Delta x$ and $x-\Delta x$ are symmetric to the boundary $\partial\Omega$. Then we have

If $X_{t_{n+1}}\in\Omega$, let $X^{\prime}_{t_{n+1}}=X_{t_{n+1}}$, and if $X_{t_{n+1}}\notin\Omega$, let $X^{\prime}_{t_{n+1}}\in\Omega$ is the symmetric point of $X_{t_{n+1}}$ to the boundary $\partial\Omega$. The $X^{\prime\prime}_{t_{n+1}}$ is used to denote the intersection of the line segment $\overline{X_{t_{n+1}}X^{\prime}_{t_{n+1}}}$ and $\partial\Omega$. Therefore, the discretization can be rewritten similarly as

where $\Delta X_{t_{n+1}}=X_{t_{n+1}}-X^{\prime}_{t_{n+1}}$. The algorithm of the proposed scheme is similar as Algorithm 1.

We consider the following Cahn-Hilliard equation, which has fourth order derivatives,

where $\phi=\phi(t,x)$ is the unknown, e.g., the concentration of the fluid, $\mu=\mu(t,x)$ is a function of $\phi$, e.g.,the chemical potential, $L_{d}>0$ is the diffusion coefficient and $\gamma>0$ is the model parameter. A first order stabilized scheme [44] for the Cahn-Hilliard equation (22) reads as

where $S$ is a suitable stabilized parameter. It is easy to derive the following equation

The first equation of (23) and (24) can be regarded as the discretization of the following modified Cahn-Hilliard equation in $[0,T]$

where $\delta=O(\Delta t)$. By reversing the time and defining

the $(\phi,\mu)$ satisfies the following backward Cahn-Hilliard equation in $[0,T]$

In order to satisfy the nonlinear Feynman-Kac formula and utilize the FBSNNs, we treat the backward Cahn-Hilliard equation (26) as a semilinear parabolic differential equation

with $\psi=(\phi,\mu)^{T}$, $A=\left(\begin{array}[]{ll}0&L_{d}\\ -\frac{\gamma^{2}}{\delta}&S\end{array}\right)$ and $F=\left(f,\frac{S}{L_{d}}f-\frac{1}{\delta}(\mu+\phi-\phi^{3})\right)^{T}$. Supposing $\lambda_{1}$ and $\lambda_{2}$ are two different eigenvalues of the coefficient matrix $A$, then the coefficient matrix $A$ can be diagonalized by $R$ and $R^{-1}$ so that $D=\mbox{diag}(\lambda_{1},\lambda_{2})=R^{-1}AR$, where $R$ is a matrix of eigenvectors. The system (27) becomes

where $\hat{\psi}=R^{-1}\psi=(\hat{\phi},\hat{\mu})^{T}$, $\hat{F}=R^{-1}F=(\hat{F}^{\hat{\phi}},\hat{F}^{\hat{\mu}})^{T}$ and $S$ is chosen so that $S>2\gamma\sqrt{\frac{L_{d}}{\delta}}$ for $\lambda_{1,2}>0$, where $\lambda_{1,2}=\frac{S\pm\sqrt{S^{2}-\frac{4\gamma^{2}L_{d}}{\delta}}}{2}$.

Therefore, the system (28) is decomposed into two independent PDEs and the corresponding FBSDEs can be obtained as follows

where

and

The $\{X_{s}^{\hat{\phi}}\}_{0\leq s\leq T}$ and $\{X_{s}^{\hat{\mu}}\}_{0\leq s\leq T}$ are the forward stochastic processes corresponding to $\hat{\phi}$ and $\hat{\mu}$ respectively, which are constrained by $x_{0}$ at the initial time.

Given a partition of $[0,T]:0=t_{0}<t_{1}<...<t_{N}=T$, we consider the simple Euler scheme for the FBSDEs (29) with $n=0,...,N-1$

where $\Delta t_{n}=t_{n+1}-t_{n}=\frac{T}{N}=\delta,\ \Delta W_{t_{n}}^{\hat{\phi}}=W_{t_{n+1}}^{\hat{\phi}}-W_{t_{n}}^{\hat{\phi}}$ and $\Delta W_{t_{n}}^{\hat{\mu}}=W_{t_{n+1}}^{\hat{\mu}}-W_{t_{n}}^{\hat{\mu}}$. The $(Y_{t_{n}}^{\hat{\phi}}(X_{t_{n}}^{\hat{\phi}}),Y_{t_{n}}^{\hat{\mu}}(X_{t_{n}}^{\hat{\mu}}))^{T}$ represents the estimated value of $(\hat{\phi}(t_{n},X_{t_{n}}^{\hat{\phi}}),\hat{\mu}(t_{n},X_{t_{n}}^{\hat{\mu}}))^{T}$. The $(\widetilde{Y}_{t_{n+1}}^{\hat{\phi}}(X_{t_{n+1}}^{\hat{\phi}}),\widetilde{Y}_{t_{n+1}}^{\hat{\mu}}(X_{t_{n+1}}^{\hat{\mu}}))^{T}$ is the reference value of $(Y_{t_{n+1}}^{\hat{\phi}}(X_{t_{n+1}}^{\hat{\phi}}),$ $Y_{t_{n+1}}^{\hat{\mu}}(X_{t_{n+1}}^{\hat{\mu}}))^{T}$, which is obtained from the last two equations in (30).

The $(Y_{t_{n}}^{\phi},Y_{t_{n}}^{\mu})^{T}$ represents the estimated value of $(\phi,\mu)^{T}$ at time $t_{n}$ given by the FNN. Due to the diagonalization, we have

where subscript $1$ or $2$ is used to represent the $1$-th or $2$-th component of the vector. We utilize $K$ different initial sampling points for training the FNN. The algorithm of the proposed scheme is summarized as Algorithm 2. Illustration of the Algorithm 2 for solving the Cahn-Hilliard equation is shown in Figure 2.