数理解析研究所講究録
1129 巻 2000 年 115-128

115

Neuralnet Collocation Method for Solving Differential Equations
Tatsuoki TAKEDA, Makoto FUKUHARA, Ali LIAQAT
Department of Information Mathematics and Computer Science
The University of Electro-Communications
電気通信大学竹田辰興、 福原誠、 リアカト・アリ

Abstract
Collocation method for solving differential equations by using a multilayer neural network is described. Possible applications of the method areinvestigated by considering the distinctive features of the method. Dataassimillation is one of promising application fields of this method.
1.
1.1

Introduction
Neural network
In this article we present an application of a feed-forward multi-layer
neural network (neuralnet) as a solution method of a differential equation.
The feed-forward multi-layer neural network is a system composed of (1)simple processor units (PUs) placed in layers and (2) connections between
the PUs of adjacent layers (Fig.1). Data flow from the input layer to the
output layer through the hidden (intermediate) layers and the connections
to which weights determined by a training process are assigned. At each PUof the hidden layers and the output layer weighted sum of the incoming
data is calculated and the result is nonlinearly transformed by an activ ation
function and transferred to the next layer. Usually some kind of sigmoid
functions is used as the activation function. But at the output layer often
the nonlinear transformation is omitted. The neural network is usually
trained as follows. (1) A lot of datasets composed of input data andcorresponding output data (supervisor data) are prepared. (2) Output data
are calculated by entering the input data to the input PUs according to theabove data flow. (3) Sum of squared errors of the output data from thesupervisor data is calculated. (4) eights assigned to the connections areadjusted so that the above squared sum is minimized. The mostcommonly used algorithm to this process is the gradient method and, inthe case of the feed-forward multi-layer neural network the error back$\mathrm{W}$116
propagation method based on the gradient method is implemented
efficiently.

Fig.1 A schematic diagram of a three-layered neural network
Features of a neural network
The feed-forward multi-layer neural network is characterized by
keywords of training, apping, smoothing, and interpolation. In relation
with the keywords, we make use of the following distinctive features of theneural network for our purpose.
(1)
Training . Optimization of the set of the weights with respect to an
appropriate object function is carried out. The object function of the
1.2$\mathrm{m}$$\cdot$neural network is represented as$E= \sum_{\mathrm{t}l^{nt}em}F(\vec{f}(\vec{x}_{t},\vec{p}n\ell em’)\mathrm{I}$is an input data vector for the trainin and is a vector
representing the weights. For the feed-forward multi-layer neuralnetwork the following type of object functions is usually employed.

where$\vec{p}$$\vec{x}_{\mu \mathrm{f}tem}$$\mathrm{g}$$\mathrm{p}(\vec{f}(_{\vec{X}_{\}nttem}}’\vec{\mathcal{P}}))=[\vec{f}(\vec{X}\vec{p}\mu ttem’)-\vec{f}_{\mu t}tem]\wedge 2$

Usually a sum of squared errors of output data from the preparedsupervisor data is used for the obj ect function as abov e, but it should be117

remarked that there are other possibilities for choice of the object
function. The training by the error back propagation method is
carried out as$\vec{p}\Rightarrow\vec{p}-\alpha^{\frac{\partial E}{\phi^{arrow}}}$(2). By the neural network mapping from the input space to the$\mathrm{M}\mathrm{a}_{\mathrm{P}\mathrm{P}^{\mathrm{i}\mathrm{g}}}\mathrm{n}\cdot$

output space good approximation can be attained in mparison with ausual orthogonal function expansion. It is because in this method not
only expansion coefficients but also the basis functions of the
expansion themselves are optimized.
Smoothing and interpolation: These are carried out in the meaning of
the least square fitting. It is important that these are always associatedwith the above neural network mapping processes.
$\mathrm{c}o$

(3)

1.3

Object function and aims of our study
In our study we applied the feed-forward multi-layer neural network
to the solution method of differential equations. Milligen et al., andLagaris et al., have proposed the method and demonstrated successful
results by the method $[1,2]$ . This is a kind of the collocation methods, i.e.,
the whole procedure is described as follows.
(1) From the computational domaina subdomaincomposed ofcollocation points is constructed.
(2) The coordinates of a collocation point are regarded as a pattern for the
input data of the neural network and a solution value corresponding to
the coordinates is regarded as the output data of the neural network.
(3) The residual of the differential equation is calculated from the output
data and the squared residual is used as the object function. In thisprocess it should be noted that the residual is expressed malyticallybecause the output data are the analytical functions of the input
variables with the weight variables as parameters.
Difference of the methods by Milligen et al., and Lagaris et al., is in thetreatment of theconditions. In the method
etal., these conditions are treated as penalty terms in the object functions.
These conditions are, therefore, satisfied only approximately. In themethod by Lagaris et al., on the other hand, these conditions are exactly
satisfied by employing appropriate form factors multiplied to the neural$\Omega$$\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{r}\mathrm{y}/\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}$$\hat{\Omega}$$\mathrm{b}\mathrm{y}\overline{\mathrm{M}}\mathrm{i}\mathrm{l}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{e}\mathrm{n}$118
network output variables. However, it is rather difficult or impossible tofind form factors appropriate for a given problem.
Moreover, the solution method by using the neural network is
generally time-consuming in comparison with the well-studied
and it does notconventional solution methods of differential equati
seem useful to apply the neural network solution method to a usual
problem for solving differential equations.
Our aims of the study are described as follows.
(1) Extension of the method by Lagaris et al., to a method with subdomaindefined form factors: As it is difficult to find a single appropriate form
conditions we divide
factor which describes the whole
the domain into a number of subdom ains and construct a set of a form
factor in each subdomain. We apply the method to a simple model$o\mathrm{n}\mathrm{s}$$\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{r}\mathrm{y}/\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}$(2)problem.
Application of the neural network collocation method for data
assimilation problems: One of the distinctive features of the neuralnetwork differential equation solver is that a smooth solution can benoise is induded in
obtained even if
conditions and so on.
constraining conditions such as
By making use of this feature the method can be applied efficiently tothe data assimilation. We study basic problems concerning thisapplication.
$\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{a}\mathrm{l}/\mathrm{o}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{a}\mathrm{l}$$\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}/\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{r}\mathrm{y}$Neuralnet collocation method (NCM)
Analytical expression of derivatives of solution
In this subsection we consider a three-layered neural network. The
result can be extended easily to a network with more layers. An analyticalexpression of a solution is given as2.
2.1$y_{k}= \sum_{j=1}^{f}w^{(}kj(_{\mathrm{i}=}2)\sum wX\sigma I1j(1i)i+w^{(}j01)1$

are the input variable, the output variable, the weightvariable, and the activation function (usually a sigmoid function:is the offset.
, respectively.
By differentiating the solution with respect to the input variables aderivative of the solution are obtained aswhere$\mathrm{x},$, and$\mathrm{w}$$\mathrm{y},$$\sigma$$\sigma(\mathrm{x})=$$1/(1+\mathrm{e}^{\mathrm{X}})-)$$\mathrm{w}_{\mathrm{i}^{0}}$119$\frac{\phi_{k}}{\partial\kappa_{\ell}}=\sum_{-,j1}^{j}w_{kj}^{(2}w_{j\ell}-)(1)\sigma 1(_{i-}\sum_{- 1}^{I}w_{ji}^{()}1xi^{+w^{(1})}j0)$$\frac{\partial_{\mathcal{Y}_{k}}^{2}}{\partial\kappa_{\ell}\ _{m}}=\sum_{j--1}wwjkjj\ell(2)(1)w^{(1)}jm0\mathrm{J}(_{i}\sum_{=1}w^{(}Xw_{j0}^{(1)}Ijii^{+}1))$

where. By using these expression the residual of theifferential equation is expressed analytically.
$\sigma^{\mathrm{t}}=d\sigma(\mathrm{x})/d\mathrm{X}$

2.2$\mathrm{d}$Collocation method
We consider the following differential equation.
$\vec{x}\in\Omega$$D^{arrow}\mathrm{y}=\vec{g}(_{\vec{X})},$

whereis an operator. We assume the output value of the neuralnetwork is an approximation of the solution of the differential equation, i.e.,$\mathrm{D}$$\vec{y}\simeq\vec{f}(\vec{x},\vec{p}\rangle$

is the weights expressed as a vector.
therefore, expressed as
where$\vec{p}$The object function is,$\mathrm{E}=\int\{D\vec{f}(\vec{X},\vec{\mathcal{P}})-s^{(\vec{X}})\}^{2}d_{\vec{X}}$Then the computational domainof theis discretized to a domain
collocation points and the object function for the collocation method is
obtained as$\hat{\Omega}$$\Omega$$E= \sum_{pat\mathrm{f}e\Gamma n}F(\vec{f}(\vec{\chi}’\vec{p}pat\prime ern))$$F(\vec{f}(\vec{x}_{p}\vec{p}at\prime ern’)\rangle=\{D\vec{f}(\vec{X}\vec{p}paftern’)-\vec{g}(X_{p\mathrm{f}m})ate\}$$\vec{x}_{p\prime}\in atten\hat{\Omega}\subset\Omega$

2.3

conditions
The boundary conditions of the Dirichlet type are given as follows.
$\mathrm{I}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}/\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{r}\mathrm{y}$$\vec{y}(\vec{x}_{b})=\vec{y}_{b}$

where

,

$\vec{x}_{b}\in\Gamma$is the boundary of the computational domain. Essentially thesame discussion holds for other kinds of boundary conditions. Initialconditions are also treated similarly.
According to the Milligen’s method theconditions
are imposed by preparing the following object function.
$\Gamma$$\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}/\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{r}\mathrm{y}$$\hat{E}=E+w_{b}$$E_{b}= \sum_{i}(\vec{y}(_{\vec{X}_{b\mathrm{t}})-_{\hat{\vec{\mathcal{Y}}}}}\wedge bi)^{2}$

It should be noted that in this expression the points where the constraining

120

conditions are given should not necessarily be placed on the computationalboundary.
In the case of the Lagaris’ method the solution of the differential
by using the output data of the neural
equation is approximated as$f(\vec{x},\vec{p})$network$\vec{f_{NN}}(\vec{x},\vec{p})$as$\vec{f}(\vec{x},\vec{p})=\vec{A}(\vec{x})+\vec{H}(\vec{x},\vec{f}_{N}N(\vec{X},\vec{p}))$$\vec{A}(\vec{x}_{b})=\vec{y}_{b}$$\vec{H}(\vec{x}_{b},\vec{f}_{NN}(_{\vec{X}\vec{p})}b’)=0$

For the Dirichlet boundary conditions
product of an$\vec{H}(\vec{x}_{b},\vec{f}_{NN}(\vec{\chi}\vec{p}b’)\mathrm{I}$is expressed as achosen form factor and the output data of the
In this case the form factor vanishes at the$\mathrm{a}\mathrm{p}\mathrm{p}_{arrow}\mathrm{r}\mathrm{o}\mathrm{p}\dot{\mathrm{n}}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{l}\mathrm{y}$neural network
boundary.

$f_{\mathrm{N}N}(\vec{x},\vec{p})$.

Boundary condition assignment by divided form factors
In this section we describe the method to apply the Lagaris’ method toproblems with a complicated boundary shape. This is realized by dividing
the whole computational domain into subdomains with simpler boundary
solve the Poisson equation in a square domain and in a Tshapes.
shaped domain. The Lagaris’ method cannot be applied directly to the
problem in the -shaped domain and we divide the -shaped domain into 7

3.

$\mathrm{W}\mathrm{e}$$\mathrm{T}$$\mathrm{T}$subdomains.

Method of a divided form factor
The conditions which the form factor used for the Dirichlet typeboundary condition should satisfy are summarized as 1) it vanishes on the)in the domain, 3) it should
boundary, 2) it is positive (orbe sufficiently smooth, and 4) derivatives are easy to calculate. As far asthese conditions are satisfied choice of the form factor is arbitrary. In thefollowing we consider two examples, i.e., the case of the square dom ain andthe case of the -shaped domain.
3.1$\mathrm{g}$$- \mathrm{d}\mathrm{e}\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{e}$$\mathrm{n}\mathrm{e}\mathrm{g}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}$$\mathrm{T}$(1) Example 1 (the square domain:$(-1,1)\cross(-\mathrm{L}1)$) (Fig.2)

In this case the simplest form factor is expressed as$g(X,\mathrm{y})=(1+X)(1-\chi\rangle(1+\nu)(1-\nu)$

121

$x$

Fig.2 The square domain.
(2) Example 2 (the$\mathrm{F}\mathrm{l}\mathrm{g}.3$

lhe 1-shaped domain.

-shaped domain: $(-1,1)\cross(0,1)\cup(-1/2,1/2)\cross(-1,1)$ ) (Fig.3)
For this case an example of the divided form factor is given as$\mathrm{T}$$x(1+x)\mathrm{y}(1-y)$$(\overline{\Omega}_{1})$$x(1-x)y(1-y)$$(\overline{\Omega}_{2})$$y(1-y)/4$$g(X,\nu^{)}=\{$$(1/4-x)2(1-\mathrm{y})2/4$$(\overline{\Omega}_{4})$$r_{1}(1-r_{1})/4$$(\overline{\Omega}_{5})$$r_{1}(1-\gamma)1/4$$(\overline{\Omega}_{\text{\’{o}}})$

1/16$(_{r_{2}\sqrt{(_{X}-1/2)^{2}+\nu^{2}}}^{r_{1}}=\sqrt{(_{X}+1/2)^{2}+y2})=$

3.2$(\overline{\Omega}_{3}\rangle$$(\overline{\Omega}_{7})$( $g:C^{1}$class)

Model problems
As a model equation we consider the following Poisson equation,$\Delta u=-(20X-3\frac{15}{2}\chi\lambda y-3\nu)-6(X^{5}-\frac{5}{4}X^{3}+\frac{1}{4}\chi\psi$in$\Omega$on
XVe solve the above equation for the two types of computational domains$u=0$(Fig.4) described in the previous section.

$\Gamma$

The exact solution is given for theabove problems as$u(x,y)=(x-1)(X- \frac{1}{2})\chi(X+1)(x+\frac{1}{2}\rangle(\mathrm{y}-1)y(\mathrm{y}+1)$

122

Fig.4 The computational domains for the model problems.
For simplicity of programing we employed a training algorithm composed ofa random search and the linear least square method instead of theconventionaly used error backpropagation method. The structure of theneural network is composed of two PUs in the input layer, 40 PUs in thelayer, and 1 PU in the output layer (we denote this structure as 2-401). The results of the calculations are shown in Figs.5-8 and the magnitudeof the errors is summarized in Table 1, where the maximum random value
denotes the maximum range of the initially assigned weights.
$\mathrm{l}\dot{\mathrm{u}}\mathrm{d}\mathrm{d}\mathrm{e}\mathrm{n}$Tahle 1$\mathrm{F}_{\lrcorner}1^{\cdot}t0\mathrm{r}_{\mathrm{l}\mathrm{s}}$of solutions of the two model problems

123

(a)

(b)

(c)

(a)

$(\mathrm{b}\rangle$(c)

Fig.5 Birdeye view of the solution
in the square domain

Fig.6 Birdeye view of the solutionin -shaped domain(a) NN collocation, (b) exact,(a) NN collocation, (b) exact.

(c) absolute error(c) absolute error$\mathrm{T}$124

Numerical solution
(b)

(a)

Fig.7 Contour plot of the solution (a) in the square domain and thecorresponding error.

Absolute error

Numerical solution

(b)

(a)

Fig.8 Contour plot of the solution (a) in the -shaped domain and thecorresponding error.
$\mathrm{T}$Application of NCM to a data assimilation problem
In this chapter we study basic issues of the neuralnet collocation
method relating to the data assimilation problems.
conditions
4.1
Irregular
In order to investigate the possibility to apply the neuralnet collocation
method to the data assimilation problem we solve differential equations
with constraining conditions given irregularly in space or time.

4.

$\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{r}\mathrm{y}/\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}$125Solution of Lorenz equation
In this section we solve the Lorenz equation by assigning initialconditions at different temporal points,
toY’ md for three variables, X,, and, respectively.
For training the network we employed a
combination of the error back-propagation method and the quasi-N wton
method. The Lorenz equation is
4.2$\mathrm{t}_{0\cross},$$\mathrm{Y}$$\mathrm{t}_{0\mathrm{Z}}$$\mathrm{Z}$$e$$\frac{dX}{dt}=\sigma(Y-\mathrm{x})$$\frac{dY}{dt}=r\mathrm{X}-\gamma-\mathrm{X}\mathrm{Z}$$\frac{dZ}{dt}=\mathrm{X}Y-b\mathrm{z}$$X(t_{0\mathrm{X}})=\mathrm{x}_{0},$$Y(t_{0Y})=Y_{0},$ $Z(t_{0\mathrm{z}^{)}}=Z_{0}$The object function of the neuralnet collocation method is given as follows.

$E= \sum_{t\in\hat{D}}\wedge\lfloor_{- t}\wedge$$\hat{D}\subset Dand|\hat{D}|<\infty$$X_{t}=A+(t-t_{0\mathrm{x}})f_{\mathrm{N}}N(\mathrm{f},\vec{p}_{\mathrm{X}})$$Y_{t}=B+(t-t_{0)}’)f_{N}N(t,\vec{\mathcal{P}}\gamma)$$Z_{t}=C+(t-\mathrm{f}_{0\mathrm{Z}})f_{NN}(t,\vec{\mathcal{P}}\mathrm{z})$$t\in D,$

$D=[0,0.4]$ ,$A=x(t_{\mathrm{o}\mathrm{x}}),$$B=\mathrm{Y}(to\mathrm{Y}),$$C=Z(t_{0})\mathrm{z}$At first, above equation is solved by the Runge-Kutta method for thefollowing conditions to obtain a base solution for comparison of theneuralnet collocation method.
$t\in D=[0,0.4]$$X_{0}=0.0,$ $Y_{0}=1.0,$ $Z_{0}=0.0$$\sigma=1.0,$,$t_{0\mathrm{X}}=t_{0\mathrm{Y}}=\mathrm{f}_{0Z}=t_{0}=0,0$$r=28.0,b=8/3$

Next, we solve the same problem by the neuralnet collocation method

imposing the initial condition at. The results by the Runge-Kutta
method and the neuralnet collocation method are shown in Fig.9, where
they are indistinguishable each other in this scale. The error of the results
is 2.22E-2 after the error back-propagation, and 2.16E-3 after the quasi$\mathrm{t}=0$

126

Newton procedure.
Then we solved the Lorenz equation for initial conditions imposed
separately to different variables, X, , and Z. The initial conditions wereobtained from the base solution by the Runge-Kutta method as$\mathrm{Y}$$(X_{0}(t_{0\mathrm{X}}),Y_{0}(t_{0\mathrm{Y}}),Z0(t_{0\mathrm{z}}))$$=(0.1961(0.2),1.0348^{0}\cdot 1),0.0011(0.05))$

The$\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{D}\mathrm{u}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{i}_{\mathrm{o}\mathrm{n}}\mathrm{a}1$results are shown in$\mathrm{F}\mathrm{i}\epsilon.10$

la)$\iota \mathrm{a}|$

1D1

(Dl

(C)

(c)

Fig. 9 Solutions of the Lorenz equation
for the same initial conditions for, and.
$\mathrm{X}(\mathrm{a}),$$\mathrm{Y}(\mathrm{b})$$\mathrm{Z}(\mathrm{c})$with the base solutions,Fig.10 Solutions of the Lorenz
equation for separately
given initial conditions.

127

Solution of a heat equation
In order to investigate a possibility to apply the neuralnet collocation
method to the data assimilation problem we consider to solve the heat
equation for irregularly imposed initial conditions.
4.3$\alpha^{2}\frac{\partial^{2}u}{h^{2}}+\sin 3m-\frac{\partial u}{\partial t}=0$,$(x,t)\in D=[0,1]\cross[0.1]$$\alpha=1$$u(\mathrm{O},t)=u(1,t)=0.\mathrm{O}$$u(x, 0)=\sin\pi x$

The analytical solution to the above problem is given as$u(x,t)= \exp(-(m)2t)\sin\pi X+\frac{1-\mathrm{e}\mathrm{x}_{\mathrm{P}}(-(3m)2t)}{(3\pi\alpha)^{2}}\sin 3\pi x$

The same problem was solved by the neuralnet
method. The
result is shown with the analytical solution in Fig. 11. Results with theaverage errors of 3.49E-3 and 1.33E-3 are obtained after 20000 iteration of
the ba&-propagation method, and 5 iterations of the quasi-Newton method,respectively. By giving initial conditions forplaced spatial
points almost similar results were obtained.
$\mathrm{c}\mathrm{o}\mathbb{I}_{0}\mathrm{c}\mathrm{a}\mathrm{t}\mathrm{i}_{\mathrm{o}\mathrm{n}}$$\mathrm{i}\mathrm{r}\mathrm{r}\mathrm{e}\omega \mathrm{a}\mathrm{r}\mathrm{l}\mathrm{y}$晶億億$\mathrm{m}\iota \mathrm{K}$h’–輸 d.

–$\blacksquare$

5.
(1)

Discussion and conclusion
Extension of Lagaris’ method to a case of subdomain-defined form factor

128

conditions is successfully applied to solution of
for thePoisson equation in a square domain and in a -shaped domain. To
attain higher accuracy some improvements are necessary.
In order to study a poSsibilit‘y to apply the neuralnet collocation method
for the data assimilation problem simple model problems on theLorentz equation and the heat equation were solved. Promising
results were obtained.
$\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}/\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{r}\mathrm{y}$$\mathrm{T}$(2)

References[1]
B.Ph. van Milligen, et al., Phys. Rev. Letters 75 (1995) 3594-3597.
[2]
I.E. Lagaris, et al., Comput. Phys. Commun. 104 (1997) 1-14