**3. A First Order ODE Example**

Let us consider the first order ODE, describing the initial value problem (IVP) given by

$$
\mathcal{L}y - f = y' + \mathcal{L}xy = 0, \quad x > 0 \tag{8a}
$$

$$\mathcal{B}y - \mathbf{b} = y(0) - 1 = 0,\tag{8b}$$

with the exact solution *y* = *e*−*x*<sup>2</sup> . By means of (6), a cost function suitable for solving (8) may be generated by

$$\begin{split} I = & \left\langle \left[ \sum\_{i=1}^{N} w\_i^{(1)} w\_i^{(0)} \phi \left( w\_i^{(0)} \mathbf{x} + b\_i^{(0)} \right) \left( 1 - \phi \left( w\_i^{(0)} \mathbf{x} + b\_i^{(0)} \right) \right) + \\ & 2 \mathbf{x} \left( \left( \sum\_{i=1}^{N} w\_i^{(1)} \phi \left( w\_i^{(0)} \mathbf{x} + b\_i^{(0)} \right) \right) + b^{(1)} \right) \right]^2 \right\rangle + \left[ y(0) - 1 \right]^2 \end{split} \tag{9}$$

The solution of (8) can be obtained by implementing a training routine which iteratively finds the set of weights *w* and bias *b* that minimises (9) (and similarly for (19) minimising (17)). The most well-known of these is the Gradient Decent method attributed to Cauchy, who first suggested it in 1847 [12]. For an overview, see, e.g., [13].

As mentioned in the previous section, the derivatives of (4) w.r.t. to the weights *w* and bias *b* are required to find them, and *automatic differentiation* is, normally, employed to perform the differentiation. However, here we carry out symbolic differentiation to demonstrate exactly the explicit expressions that constitute the gradient of the cost function. Indeed, by taking the partial derivatives we obtain

$$\frac{\partial y}{\partial w\_i^{(0)}} = \frac{\partial}{\partial w\_i^{(0)}} \left( \left( \sum\_{i=1}^N w\_i^{(1)} \phi \left( w\_i^{(0)} \mathbf{x} + b\_i^{(0)} \right) \right) + b^{(1)} \right) = w\_i^{(1)} \phi' \left( w\_i^{(0)} \mathbf{x} + b\_i^{(0)} \right) \mathbf{x}, \tag{10a}$$

$$\frac{\partial y}{\partial w\_i^{(1)}} = \frac{\partial}{\partial w\_i^{(1)}} \left( \left( \sum\_{i=1}^N w\_i^{(1)} \phi \left( w\_i^{(0)} \mathbf{x} + b\_i^{(0)} \right) \right) + b^{(1)} \right) = \phi \left( w\_i^{(0)} \mathbf{x} + b\_i^{(0)} \right), \tag{10b}$$

$$\frac{\partial y}{\partial b\_j^{(0)}} = \frac{\partial}{\partial b\_j^{(0)}} \left( \left( \sum\_{i=1}^N w\_i^{(1)} \phi \left( w\_i^{(0)} \mathbf{x} + b\_i^{(0)} \right) \right) + b^{(1)} \right) = w\_i^{(1)} \phi' \left( w\_i^{(0)} \mathbf{x} + b\_i^{(0)} \right), \tag{10c}$$

$$\frac{\partial y}{\partial b^{(1)}} = 1.\tag{10d}$$

Moreover, the derivatives of the cost function (5) w.r.t. to the weights and bias are also required. For the derivative w.r.t. *w*(0) *<sup>i</sup>* for the first order ODE (8), this means that

$$\left\langle 2(y' + 2xy) \left( \frac{\partial y'}{\partial w\_i^{(0)}} + 2x \frac{\partial y}{\partial w\_i^{(0)}} \right) \right\rangle + 2(y(0) - 1) \frac{\partial y(0)}{\partial w\_i^{(0)}}.\tag{11}$$

To complete the analysis, we also need expressions for the derivatives of *y* w.r.t. *w*(0) *<sup>i</sup>* , *<sup>w</sup>*(1) *<sup>i</sup>* , *b* (0) *<sup>i</sup>* and *<sup>b</sup>*(1). By the chain rule, the following expressions can be obtained, viz.

$$\begin{split} \frac{\partial \boldsymbol{y}^{\prime}}{\partial \boldsymbol{w}\_{i}^{(0)}} &= \frac{\partial}{\partial \boldsymbol{w}\_{i}^{(0)}} \sum\_{i=1}^{N} \boldsymbol{w}\_{i}^{(1)} \boldsymbol{w}\_{i}^{(0)} \boldsymbol{\phi}^{\prime} \left( \boldsymbol{w}\_{i}^{(0)} \boldsymbol{x} + \boldsymbol{b}\_{i}^{(0)} \right) = \\ \boldsymbol{w}\_{i}^{(1)} \boldsymbol{\phi}^{\prime} \left( \boldsymbol{w}\_{i}^{(0)} \boldsymbol{x} + \boldsymbol{b}\_{i}^{(0)} \right) &+ \boldsymbol{x} \boldsymbol{w}\_{i}^{(1)} \left( \boldsymbol{w}\_{i}^{(0)} \right)^{2} \boldsymbol{\phi}^{\prime\prime} \left( \boldsymbol{w}\_{i}^{(0)} \boldsymbol{x} + \boldsymbol{b}\_{i}^{(0)} \right), \end{split} \tag{12a}$$

$$\frac{\partial y'}{\partial w\_i^{(1)}} = \frac{\partial}{\partial w\_i^{(1)}} w\_i^{(1)} w\_i^{(0)} \phi' \left( w\_i^{(0)} \mathbf{x} + b\_i^{(0)} \right) = w\_i^{(0)} \phi' \left( w\_i^{(0)} \mathbf{x} + b\_i^{(0)} \right), \tag{12b}$$

$$\frac{\partial y'}{\partial b\_i^{(0)}} = \frac{\partial}{\partial b\_i^{(0)}} \sum\_{i=1}^N w\_i^{(1)} w\_i^{(0)} \phi' \left( w\_i^{(0)} \mathbf{x} + b\_i^{(0)} \right) = w\_i^{(1)} w\_i^{(0)} \phi'' \left( w\_i^{(0)} \mathbf{x} + b\_i^{(0)} \right), \tag{12c}$$

$$\frac{\partial y'}{\partial b^{(1)}} = 0.\tag{12d}$$

What remains now is to obtain expressions for *y*(0) and the partial derivatives of *y*(0), w.r.t. to the weights and bias. Let us start with *y*(0). With *y*(*x*) given by (4) we directly have

$$\log(0) = \left(\sum\_{i=1}^{N} w\_i^{(1)} \phi \left(b\_i^{(0)}\right)\right) + b^{(1)}\text{.}\tag{13}$$

which, in turn, means that

$$\frac{\partial \hat{y}(0)}{\partial w\_i^{(0)}} = 0,\tag{14a}$$

$$\frac{\partial \mathfrak{y}(0)}{\partial w\_i^{(1)}} = \frac{\partial}{\partial w\_i^{(1)}} \left( \left( \sum\_{i=1}^N w\_i^{(1)} \mathfrak{o} \left( b\_i^{(0)} \right) \right) + b^{(1)} \right) = \mathfrak{o} \left( b\_i^{(0)} \right), \tag{14b}$$

$$\frac{\partial y(0)}{\partial b\_i^{(0)}} = \left( \left( \sum\_{i=1}^N w\_i^{(1)} \phi \left( b\_i^{(0)} \right) \right) + b^{(1)} \right) = w\_i^{(1)} \phi' \left( b\_i^{(0)} \right) \tag{14c}$$

$$\frac{\partial y(0)}{\partial b^{(1)}} = 1.\tag{14d}$$

The PINN (following the architecture presented above) was implemented as computer program in MATLAB. The program was employed to obtain a numerical solution to the IVP in (8), using the parameters in Table 1.

**Table 1.** Parameters used to defined the PINN to for the IVP in (8).


The weights *w*(0) *<sup>i</sup>* and bias *b* (0) *<sup>i</sup>* were initialised using randomly generated and uniformly distributed numbers in the interval [−2, 2], while the weights *<sup>w</sup>*(1) *<sup>i</sup>* were initially set to zero and the bias *b*(1) to one, to ensure fulfilment of the initial condition (*y*(0) = 1).

Table 2 lists the weights an bias corresponding to the solution presented in Figure 2. We note that, with the weights and bias given by Table 2, the trained network's prediction exhibits the overall error

$$\frac{1}{N\_i} \sqrt{\sum\_{k=1}^{N\_i} \left( \varepsilon^{-x\_k^2} - y(x\_k) \right)^2} = 5.8 \times 10^{-4} \text{ } \tag{15}$$

and 1 − *<sup>y</sup>*(0) = 2.2 × <sup>10</sup><sup>−</sup>4, when comparing against the initial condition.


**Table 2.** Parameters used to defined the PINN for the IVP (8).

**Figure 2.** The solution to the IVP (8), predicted by the PINN (red line with circle markers) and the exact solution obtained by integration (blue continuous line).

#### **4. A PINN for the Classical Reynolds Equation**

The Reynolds equation for a one-dimensional flow situation, where the lubricant is assumed to be incompressible and iso-viscous, is a second order Boundary Value Problem (BVP), which in dimensionless form can be formulated as

$$\frac{d}{d\mathbf{x}}\left(c(\mathbf{x})\frac{d\underline{y}}{d\mathbf{x}}\right) = f(\mathbf{x}), \quad 0 < \mathbf{x} < 1,\tag{16a}$$

$$y(0) = 0, \quad y(1) = 0,\tag{16b}$$

where *c*(*x*) = *H*3, *f*(*x*) = *dH dX* and *H* is the dimensionless film thickness, if it is assumed that the pressure *y* at the boundaries is zero. For the subsequent analysis it is, however, more suitable work with a condensed form which can be obtained by defining the operators L and B as

$$
\mathcal{L}y = c(x)y'' + c'(x)y',\tag{17a}
$$

$$\mathcal{B}y = \begin{bmatrix} y(0) \\ y(1) \end{bmatrix}. \tag{17b}$$

The Reynolds BVP given by (16) can then be presented as

$$
\mathcal{L}y - f = 0, \quad 0 < x < 1,\tag{18a}
$$

$$
\partial \mathbf{y} - \mathbf{b} = \mathbf{0},
\tag{18b}
$$

where **b** = **0**.

For the Reynolds BVP, the cost function (5) becomes

$$l = \left\langle \left( \varepsilon(\mathbf{x}) y'' + \varepsilon'(\mathbf{x}) y' - f \right)^2 \right\rangle + y^2(0) + y^2(1),\tag{19}$$

and from the analysis presented for the IVP in Section 3 above, we have all the "ingredients" except for the partial derivatives of *y*- and *y*(1) w.r.t. to the weights and bias. For *y*--, based on (7) and (12), we obtain

$$\begin{split} \frac{\partial \boldsymbol{y}^{\prime\prime}}{\partial w^{(0)}\_{i}} &= \frac{\partial}{\partial w^{(0)}\_{i}} \sum\_{i=1}^{N} w^{(1)}\_{i} \left( w^{(0)}\_{i} \right)^{2} \boldsymbol{\phi}^{\prime\prime} \left( w^{(0)}\_{i} \, \mathbf{x} + b^{(0)}\_{i} \right) = \\ \boldsymbol{\phi} &= 2w^{(1)}\_{i} w^{(0)}\_{i} \boldsymbol{\phi}^{\prime\prime} \left( w^{(0)}\_{i} \, \mathbf{x} + b^{(0)}\_{i} \right) + x w^{(1)}\_{i} \left( w^{(0)}\_{i} \right)^{2} \boldsymbol{\phi}^{\prime\prime\prime} \left( w^{(0)}\_{i} \, \mathbf{x} + b^{(0)}\_{i} \right), \end{split} \tag{20a}$$

$$\frac{\partial y^{\prime\prime}}{\partial w\_{i}^{(1)}} = \frac{\partial}{\partial w\_{i}^{(1)}} \sum\_{i=1}^{N} w\_{i}^{(1)} \left(w\_{i}^{(0)}\right)^{2} \phi^{\prime\prime} \left(w\_{i}^{(0)} \mathbf{x} + b\_{i}^{(0)}\right) = \left(w\_{i}^{(0)}\right)^{2} \phi^{\prime\prime} \left(w\_{i}^{(0)} \mathbf{x} + b\_{i}^{(0)}\right),\tag{20b}$$

$$\frac{\partial y^{\prime\prime}}{\partial b\_{i}^{(0)}} = \frac{\partial}{\partial b\_{i}^{(0)}} \sum\_{i=1}^{N} w\_{i}^{(1)} \left(w\_{i}^{(0)}\right)^{2} \phi^{\prime\prime} \left(w\_{i}^{(0)} \mathbf{x} + b\_{i}^{(0)}\right) = w\_{i}^{(1)} \left(w\_{i}^{(0)}\right)^{2} \phi^{\prime\prime\prime} \left(w\_{i}^{(0)} \mathbf{x} + b\_{i}^{(0)}\right), \tag{20c}$$

$$\frac{\partial y''}{\partial b^{(1)}} = 0,\tag{20d}$$

where the third derivative of the Sigmoid function (1) is required. It yields

$$\begin{split} &\frac{d}{d\tilde{\xi}}\left(\phi^{\prime\prime}(\tilde{\xi})\right) = \frac{d}{d\tilde{\xi}}\left(\phi^{\prime}(\tilde{\xi})(1-2\phi(\tilde{\xi}))\right) = \phi^{\prime\prime}(\tilde{\xi})(1-2\phi(\tilde{\xi})) - 2\left(\phi^{\prime}(\tilde{\xi})\right)^{2} = \\ &= \phi(\tilde{\xi})(1-\phi(\tilde{\xi}))(1-2\phi(\tilde{\xi}))^{2} - 2\left(\phi(\tilde{\xi})(1-\phi(\tilde{\xi}))\right)^{2} = \\ &= \phi(\tilde{\xi})(1-\phi(\tilde{\xi}))^{2}(1-3\phi(\tilde{\xi})) .\end{split}$$

For *y*(1) we obtain

$$\frac{\partial \mathfrak{y}(1)}{\partial w\_{i}^{(0)}} = \frac{\partial}{\partial w\_{i}^{(0)}} \left( \left( \sum\_{i=1}^{N} w\_{i}^{(1)} \phi \left( w\_{i}^{(0)} + b\_{i}^{(0)} \right) \right) + b^{(1)} \right) = w\_{i}^{(1)} \phi' \left( w\_{i}^{(0)} + b\_{i}^{(0)} \right), \tag{22a}$$

$$\frac{\partial y(1)}{\partial w\_i^{(1)}} = \frac{\partial}{\partial w\_i^{(1)}} \left( \left( \sum\_{i=1}^N w\_i^{(1)} \phi \left( w\_i^{(0)} + b\_i^{(0)} \right) \right) + b^{(1)} \right) = \phi \left( w\_i^{(0)} + b\_i^{(0)} \right), \tag{22b}$$

$$\frac{\partial y(1)}{\partial b\_i^{(0)}} = \frac{\partial}{\partial b\_i^{(0)}} \left( \left( \sum\_{i=1}^N w\_i^{(1)} \phi \left( w\_i^{(0)} + b\_i^{(0)} \right) \right) + b^{(1)} \right) = w\_i^{(1)} \phi' \left( w\_i^{(0)} + b\_i^{(0)} \right), \tag{22c}$$

$$\frac{\partial y(1)}{\partial b^{(1)}} = \frac{\partial}{\partial b^{(1)}} \left( \left( \sum\_{i=1}^{N} w\_i^{(1)} \phi \left( w\_i^{(0)} + b\_i^{(0)} \right) \right) + b^{(1)} \right) = 1,\tag{22d}$$

and we now have all the "ingredients" required to fully specify (19). To test the performance of the PINN, a Reynolds BVP was specified for a linear slider with dimensionless film thickness defined by

$$H(\mathbf{x}) = \mathbf{1} + K - \mathbf{K}\mathbf{x}.\tag{23}$$

This means that *<sup>c</sup>*(*x*)=(<sup>1</sup> + *<sup>K</sup>* − *Kx*)<sup>3</sup> and *<sup>f</sup>*(*x*) = *dH dx* = −*K* and that the exact solution is

$$y\_{\text{exact}}(\mathbf{x}) = \left[ \frac{1}{K} \left( \frac{1}{1+K-\mathbf{K}\mathbf{x}} - \frac{1+K}{2+K} \frac{1}{(1+\mathbf{K}-\mathbf{K}\mathbf{x})^2} - \frac{1}{2+K} \right) \right],\tag{24}$$

see, e.g., [14].

The PINN (following the architecture suggested herein) was implemented in MATLAB and a numerical solution to (16) was obtained using the parameters in Table 3. As for the IVP, addressed in the previous section, the weights *w*(0) *<sup>i</sup>* and bias *b* (0) *<sup>i</sup>* were, again, initialised using randomly generated numbers, uniformly distributed in [−2, 2], while the weights *w*(1) *<sup>i</sup>* and the bias *<sup>b</sup>*(1) was initially set to zero, to ensure fulfilment of the boundary conditions.

Figure 3 depicts solution predicted by the PINN (red line with circle markers) and the exact solution obtained by integration (blue continuous line).


**Table 3.** Parameters used to defined the ANN to for the Reynolds equation.

**Figure 3.** The solution achieved by the ANN (red line with circle markers) and the exact solution obtained by integration (blue continuous line).

Table 4 lists the weights an bias corresponding to the solution presented in Figure 3.


 −9.3674 11.4571 0.3843 −4.5473 3.3266 0.0305 −2.4464 −1.9884 0.1188 −0.1365 −0.1674 0.4155 7 0.8581 0.5253 0.5089 8 1.0901 2.0858 0.3348 9 0.2085 0.2523 −0.2024 −3.2168 5.9722 −0.9899

**Table 4.** Parameters used to defined the ANN.


$$\frac{1}{N\_i}\sqrt{\sum\_{k=1}^{N\_i} \left(y\_{\text{exact}l}(\mathbf{x}\_k) - y(\mathbf{x}\_k)\right)^2} = 6.2 \times 10^{-5} \text{,}\tag{25}$$

while *<sup>y</sup>*(0) = 4.1 × <sup>10</sup>−<sup>4</sup> and *<sup>y</sup>*(1) = −4.0 × <sup>10</sup><sup>−</sup>4.

**Remark 1.** *The formulation of the PINN presented here is applicable as a numerical solution procedure for the Reynolds BVP* (16) *and it does not consider the effect of cavitation. Exactly how the effect of cavitation can be included is, however, out of the scope of this paper.*
