A Novel Method for Solving Ordinary Differential Equations with Artificial Neural Networks

Roseline N. Okereke, Olaniyi S. Maliki, Ben I. Oruh

Department of Mathematics, Michael Okpara University of Agriculture, Umudike, Nigeria.

DOI: 10.4236/am.2021.1210059                

Abstract

This research work investigates the use of Artificial Neural Network (ANN) based on models for solving first and second order linear constant coefficient ordinary differential equations with initial conditions. In particular, we employ a feed-forward Multilayer Perceptron Neural Network (MLPNN), but bypass the standard back-propagation algorithm for updating the intrinsic weights. A trial solution of the differential equation is written as a sum of two parts. The first part satisfies the initial or boundary conditions and contains no adjustable parameters. The second part involves a feed-forward neural network to be trained to satisfy the differential equation. Numerous works have appeared in recent times regarding the solution of differential equations using ANN, however majority of these employed a single hidden layer perceptron model, incorporating a back-propagation algorithm for weight updation. For the homogeneous case, we assume a solution in exponential form and compute a polynomial approximation using statistical regression. From here we pick the unknown coefficients as the weights from input layer to hidden layer of the associated neural network trial solution. To get the weights from hidden layer to the output layer, we form algebraic equations incorporating the default sign of the differential equations. We then apply the Gaussian Radial Basis function (GRBF) approximation model to achieve our objective. The weights obtained in this manner need not be adjusted. We proceed to develop a Neural Network algorithm using MathCAD software, which enables us to slightly adjust the intrinsic biases. We compare the convergence and the accuracy of our results with analytic solutions, as well as well-known numerical methods and obtain satisfactory results for our example ODE problems.

Keywords

Ordinary Differential Equations, Multilayer Perceptron Neural Networks, Gaussian Radial Basis Function, Network Training, MathCAD (Computer Aided Design) 14, IBM-SPSS (Statistical Package for Social Science) 23

1. Introduction

The beginning of Neuro-computing is often taken to be the research article of McCulloch and Pitts [1] published in 1943, which showed that even simple types of neural networks could, in principle, compute any arithmetic or logical function, was widely read and had great influence. Other researchers, principally von Neumann, wrote a book [2] in which the suggestion was made that the research into the design of brain-like or brain-inspired computers might be of great interest and benefit to scientific and technological knowledge.

We present a new perspective for obtaining solutions of initial value problems using Artificial Neural Networks (ANN). We discover that neural network based model for the solution of ordinary differential equations (ODE) provides a number of advantages over standard numerical methods. Firstly, the neural network based solution is differentiable and is in closed analytic form. On the other hand most other techniques offer a discretized solution or a solution with limited differentiability. Secondly, the neural network based method for solving a differential equation provides a solution with very good generalization properties. The major advantage here is that our method reduces considerably the computational complexity involved in weight updating, while maintaining satisfactory accuracy.

Neural Network Structure

A neural network is an inter-connection of processing elements, units or nodes, whose functionality resemble that of the human neurons [3]. The processing ability of the network is stored in the connection strengths, simply called weights, which can be obtained by a process of adaptation to, a set of training patterns. Neural network methods can solve both ordinary and partial differential equations. Furthermore, it relies on the function approximation property of feed forward neural networks which results in a solution written in a closed analytic form. This form employs a feed forward neural network as a basic approximation element [4] [5]. Training of the neural network can be done either by any optimization technique which in turn requires the computation of the gradient the error with respect to the network parameters, by regression based model or by basis function approximation. In any of these methods, a trial solution of the differential equation is written as a sum of two parts, proposed by Lagaris [6]. The first part satisfies the initial or boundary conditions and contains no adjustable parameters. The second part contains some adjustable parameters that involves feed forward neural network and is constructed in a way that does not affect the initial or boundary conditions. Through the construction, the trial solution, initial or boundary conditions are satisfied and the network is trained to satisfy the differential equation. The general flowchart for neural network training (or learning) is given below in Figure 1.

2. Neural Networks as Universal Approximators

Artificial neural networks can make a nonlinear mapping from the inputs to the outputs of the corresponding system of neurons, which is suitable for analyzing the problem defined by initial/boundary value problems that have no analytical solutions or which cannot be easily computed. One of the applications of the multilayer feed forward neural network is the global approximation of real valued multivariable function in a closed analytic form. Namely such neural networks are universal approximators. It is discovered in the literature that multilayer feed forward neural networks with one hidden layer using arbitrary squashing functions are capable of approximating any Borel measurable function from one finite dimensional space to another with any desired degree of accuracy. This is made clear in the following theorem.

Universal Approximation Theorem

The universal approximation theorem for MLP was proved by Cybenko [7] and Hornik et al. [8] in 1989. Let $I_n$ represent an n-dimensional unit cube containing all possible input samples $x=\left(x_1,x_2,\cdots ,x_n\right)$ with $x_i\in \left[0,1\right]$, $i=1,2,\cdots ,n$. Let $C\left(I_n\right)$ be the space of continuous functions on $I_n$, given a continuous sigmoid function $\phi (\cdot )$, then the universal approximation theorem states that the finite sums of the form

$y_k=y_k\left(x,w\right)={\displaystyle \underset{i=1}{\overset{N_2}{\sum }}{w}_{ki}^3\phi }\left({\displaystyle \underset{j=0}{\overset{n}{\sum }}{w}_{ki}^2{x}_j}\right)\text{,}k=1,2,\cdots ,m$ (1)

are dense in $C\left(I_n\right)$. This simply means that given any function $f\in C\left(I_n\right)$ and $\epsilon >0$, there is a sum $y\left(x,w\right)$ of the above form that satisfies

$\left|y\left(x,w\right)-f\left(x\right)\right|<\epsilon ,\forall x\in I_n$. (2)

Figure 1. Network training flowchart.

3. Gradient Computation

Minimization of error function can also be considered as a procedure for training the neural network [9], where the error corresponding to each input vector $x$ is the value $f\left(x\right)$ which has to become zero. In computing this error value, we require the network output as well as the derivatives of the output with respect to the input vectors. Therefore, while computing error with respect to the network parameters, we need to compute not only the gradient of the network but also the gradient of the network derivatives with respect to its inputs. This process can be quite tedious computationally, and we briefly outline this in what follows.

3.1. Gradient Computation with Respect to Network Inputs

Next step is to compute the gradient with respect to input vectors, for this purpose let us consider a multilayer perceptron (MLP) neural network with n input units, a hidden layer with m sigmoid units and a linear output unit. For a given input vector $x=\left(x_1,x_2,\cdots ,x_n\right)$ the network output is written:

$N\left(x,p\right)={\displaystyle \underset{i=1}{\overset{m}{\sum }}{v}_j\phi \left(z_j\right)}$, $z_j={\displaystyle \underset{i=1}{\overset{n}{\sum }}{w}_{ji}{x}_i+u_j}$ (3)

$w_{ji}$ denotes the weight from input unit 𝑖 to the hidden unit 𝑗, $v_j$ denotes weight from the hidden unit 𝑗 to the output unit, $u_j$ denotes the biases, and $\phi \left(z_j\right)$ is the sigmoid activation function.

Now the derivative of networks output N with respect to input vector $x_i$ is:

$\frac{\partial }{\partial x_i}N\left(x,p\right)=\frac{\partial }{\partial x_i}\left({\displaystyle \underset{j=1}{\overset{m}{\sum }}{v}_j\phi \left(z_j\right)}\right)={\displaystyle \underset{j=1}{\overset{m}{\sum }}{v}_j{w}_{ji}}{\phi }^{\left(1\right)}$ (4)

where ${\phi }^{\left(1\right)}\equiv \partial \phi \left(x\right)/\partial x$. Similarly, the k^th derivative of N is computed as; ${\partial }^{k}N/\partial x_i^k={\displaystyle \underset{j=1}{\overset{m}{\sum }}{v}_j{w}_{ji}^k{\phi }_j^{(\; k\; )}}$

Where ${\phi }_j\equiv \phi \left(z_j\right)$ and ${\phi }^{\left(k\right)}$ denotes the k^th order derivative of the sigmoid activation function.

3.2. Gradient Computation with Respect to Network Parameters

Network’s derivative with respect to any of its inputs is equivalent to a feed-forward neural network $N_k\left(x\right)$ with one hidden layer, having the same values for the weights $w_{ji}$ and thresholds $u_j$ and with each weight $v_j$ being replaced with $v_j{p}_j$. Moreover, the transfer function of each hidden unit is replaced with the k^th order derivative of the sigmoid function. Therefore, the gradient of $N_k$ with respect to the parameters of the original network can easily obtained as:

3.3. Network Parameter Updation

After computation of derivative of the error with respect to the network parameter has been defined then the network parameters $v_j$, $u_j$ and $w_{ji}$ updation rule is given as,

${\nu }_j\left(t+1\right)={\nu }_j\left(t\right)+\mu \frac{\partial N_k}{\partial {\nu }_j}$, $u_j\left(t+1\right)=u_j\left(t\right)+\eta \frac{\partial N_k}{\partial u_j}$, $$w_{ji}\left(t+1\right)=w_{ji}\left(t\right)+\gamma \frac{\partial N_k}{\partial w_{ji}}$$ (5)

where $\mu$, $\eta$ and $\gamma$ are the learning rates, $i=1,2,\cdots ,n$ and $j=1,2,\cdots ,m$.

Once a derivative of the error with respect to the network parameters has been defined it is then straightforward to employ any optimization technique to minimize error function.

4. General Formulation for Differential Equations

Let us consider the following general differential equations which represent both ordinary and partial differential equations Majidzadeh [10]:

$G\left(x,\psi \left(x\right),\nabla \psi \left(x\right),{\nabla }^2\psi \left(x\right),\cdots \right)=0,\forall x\in D,$ (6)

subject to some initial or boundary conditions, where $x=\left(x_1,x_2,\cdots ,x_n\right)\in {ℝ}^n$, $D\subset {ℝ}^n$ denotes the domain, and $\psi \left(x\right)$ is the unknown scalar-valued solution to be computed. Here, G is the function which defines the structure of the differential equation and $\nabla$ is a differential operator. Let ${\psi }_t\left(x,p\right)$ denote the trail solution with parameters (weights, biases) p. Tian Qiet al. [11], gave the following as the general formulation for the solution of differential equations (5) using ANN. Now, ${\psi }_t\left(x,p\right)$ may be written as the sum of two terms

${\psi }_t\left(x,p\right)=A\left(x\right)+F\left(x,N\left(x,p\right)\right),$ (7)

where $A\left(x\right)$ satisfies initial or boundary condition and contains no adjustable parameters, whereas $N\left(x,p\right)$ is the output of feed forward neural network with the parameters p and input data x. The function $F\left(x,N\left(x,p\right)\right)$ is actually the operational model of the neural network. Feed forward neural network (FFNN) converts differential equation problem to function approximation problem. The neural network $N\left(x,p\right)$ is given by;

$N\left(x,p\right)={\displaystyle \underset{j=1}{\overset{m}{\sum }}{v}_j\sigma \left(z_j\right)}$, $z_j={\displaystyle \underset{i=1}{\overset{n}{\sum }}{w}_{ji}{x}_i+u_j}$ (8)

$w_{ji}$ denotes the weight from input unit i to the hidden unit j, $v_j$ denotes weight from the hidden unit j to the output unit, $u_j$ denotes the biases, and $\sigma \left(z_j\right)$ is the sigmoid activation function.

Neural Network Training

The neural network weights determine the closeness of predicted outcome to the desired outcome. If the neural network weights are not able to make the correct prediction, then only the biases need to be adjusted. The basis function we shall apply in this work in training the neural network is the sigmoid activation function given by

$\sigma \left(z_j\right)={\left(1+{\text{e}}^{-z_j}\right)}^{-1}$. (9)

5. Method for Solving First Order Ordinary Differential Equation

Let us consider the first order ordinary differential equation below

${\psi }^{\prime }\left(x\right)=f\left(x,\psi \right),x\in \left[a,b\right]$ (10)

with initial condition $\psi \left(a\right)=A$. In this case, the ANN trial solution may be written as

${\psi }_t\left(x,p\right)=A+xN\left(x,p\right),$ (11)

where $N\left(x,p\right)$ is the neural output of the feed forward network with one input data x with parameters p. The trial solution ${\psi }_t\left(x,p\right)$ satisfies the initial condition. Now let us consider a first order differential equation:

${\psi }^{\prime }\left(x\right)-\psi =0,x\in \left[0,1\right],\psi \left(0\right)=1$ (12)

with trial solution:

${\psi }_t\left(x,p\right)=A+xN\left(x,p\right),$ (13)

where x is the input to neural network model and p represents the parameters—weights and biases.

$w_t\left(0,p\right)=A+\left(0\right)N\left(0,p\right)=A=1$ $\Rightarrow {\psi }_t\left(x,p\right)=1+xN\left(x,p\right),$ (14)

To solve this problem using neural network, we shall employ a neural network architecture with three layers. One input layer with one neuron; one hidden layer with three neurons and one output layer with one output unit, as depicted below in Figure 2.

Each neuron is connected to other neurons of the previous layer through adaptable synaptic weights $w_j$ and biases $u_j$. Now, ${\psi }_t\left(x_i,p\right)=1+x_{i}N\left(x_i,p\right)$, with $N\left(x_i,p\right)={\displaystyle \underset{j=1}{\overset{3}{\sum }}{v}_j\sigma \left(z_j\right)}$, and

${\displaystyle \underset{j=i}{\overset{3}{\sum }}{v}_j\sigma \left(z_j\right)}=v_1\sigma \left(z_1\right)+v_2\sigma \left(z_2\right)+v_3\sigma \left(z_3\right)$, (15)

where $z_1=x_1{w}_{11}+u_1,z_2=x_1{w}_{12}+u_2,z_3=x_1{w}_{13}+u_3$.

Figure 2. Schematic for $N\left(x,p\right)$.

Now, in solving ordinary differential equations, we assume a solution to the homogeneous part and approximate the function using SPSS model, which estimates the regression coefficients in the multiple regression mode. These coefficients are what we use as the weights from the input layer to the hidden layer. The condition placed on $y_a\left(x\right)=f\left(x\right)$, say, where $y_a\left(x\right)$ is the assumed solution is that $f\left(x\right)\ne 0$.

Any exponential function, $y\left(x\right)={\text{e}}^{\alpha x}$, where $\alpha \in ℝ$, is a part of solution to any first order ordinary differential equation. We regress that function using excel spreadsheet and SPSS model as follows:

Assuming we let $y_a\left(x\right)=f\left(x\right)={\text{e}}^{2x},x\in \left[0,1\right]$, be a solution to a given first order ordinary differential equation, $y^{\prime }-2y=f\left(x\right)$, defined on the given interval. Then dividing the interval into 10 equidistant points, we use excel spreadsheet to find values of $y_a\left(x\right)$ at all the points shown in Table 1, then use SPSS to regress and get the weights which we designate weights from input layer to hidden layer.

The above is followed by the display of the SPSS 20 output of the data in Table 2, from which we pick the weights.

Looking at Table 2, we see that the cubic curve fits perfectly the assumed solution. Therefore, we pick the coefficients: 2.413, 0.115 and 3.860 as the weights from input layer to the hidden layer.

The next task is to obtain the weights from hidden layer to the output layer. We shall find $f\left(x\right)$, a real function of a real valued vector $x={\left(x_1,x_2,\cdots ,x_d\right)}^{\text{T}}$ and a set of functions, $\left\{{\phi }_i\left(x\right)\right\}$ called the elementary functions such that

$\hat{f}\left(x,v\right)={\displaystyle \underset{i=1}{\overset{N}{\sum }}{v}_i{\phi }_i\left(x\right)}$ (16)

is satisfied, where $v_i$ are real valued constants such that $\left|f\left(x\right)-\hat{f}\left(x,v\right)\right|<\epsilon$. When one can find coefficients $v_i$ that make $\epsilon$ arbitrarily small for any function $f(.)$ over the domain of interest, we say that the elementary function set $\left\{{\phi }_i(.)\right\}$ has the property of universal approximation over the class of functions $f(.)$. There are many possible elementary functions we can choose from which we will show later. Now if the number of input vectors $x_i$ is made equal to the

Table 1. Table of values for the relationship between y_a andx.

Table 2. Model summary and parameter estimates.

number of elementary functions ${\phi }_i(.)$, then the normal equations can be given as;

$\left[\begin{array}{ccc}{\phi }_1\left(x_1\right)& \cdots & {\phi }_N\left(x_1\right)\\ ⋮& \ddots & ⋮\\ {\phi }_1\left(x_N\right)& \cdots & {\phi }_N\left(x_N\right)\end{array}\right]\left[\begin{array}{c}{v}_1\\ ⋮\\ v_N\end{array}\right]=\left[\begin{array}{c}f\left(x_1\right)\\ ⋮\\ f\left(x_N\right)\end{array}\right]$ (17)

and the solution becomes $v={\varphi }^{-1}f$, where v becomes a vector with the coefficients, f is a vector composed of the values of the function at the N points, and $\varphi$ the matrix with entries given by values of the elementary functions at each of the N points in the domain. An important condition that must be placed in the elementary functions is that the inverse of $\varphi$ must exist. In general, there are many sets $\left\{{\phi }_i(.)\right\}$ with the property of universal approximation for a class of functions. We would prefer a set $\left\{{\phi }_i(.)\right\}$ over another $\left\{{\gamma }_i(.)\right\}$ if $\left\{{\phi }_i(.)\right\}$ provides a smaller error ε for a pre-set value of N. This means that the speed of convergence of the approximation is also an important factor in the selection of the basis. As we mentioned earlier, there are many possible elementary basis functions we can choose from. A typical example of basis function is the Gaussian basis specified by:

$G\left(z\right)=G\left(\Vert x-c_i\Vert \right)=\exp \left(-z^2\right);i=1,2,\cdots ,N$

However in neuro-computing, the most popular choice for elementary functions is the radial basis function (RBFs) where ${\phi }_i\left(x\right)$ is given;

${\phi }_i\left(x\right)=\exp \left(-\frac{{\left|x-x_i\right|}^2}{2{\sigma }^2}\right)$ (18)

where ${\sigma }^2$ is the variance of input vectors x. It is this last basis function that we shall adopt in generating our weights from hidden layer to output layer. At this point we shall divide a given interval into a certain number of points equidistant from each other and choose input vectors $x_i$ to conform with the number of elementary functions ${\phi }_i(.)$ to make the basis function implementable. Here N = 3 and the solution $v={\varphi }^{-1}f$ is given by;

$\left[\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right]={\left[\begin{array}{ccc}{\phi }_1\left(x_1\right)& {\phi }_2\left(x_1\right)& {\phi }_3\left(x_1\right)\\ {\phi }_1\left(x_2\right)& {\phi }_2\left(x_2\right)& {\phi }_3\left(x_2\right)\\ {\phi }_1\left(x_3\right)& {\phi }_2\left(x_3\right)& {\phi }_3\left(x_3\right)\end{array}\right]}^{-1}\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\end{array}\right]$ (19)

where $v,\varphi ,f$ are as defined before, and

${\phi }_i\left(x\right)=\exp \left(-\frac{{\left|x-x_i\right|}^2}{2{\sigma }^2}\right),{\sigma }^2=\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left(x_i-\bar{x}\right)}^2},\bar{x}=\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{x}_i}$ (20)

Now to compute $f\left(x\right)$ depends on the nature of a given differential equation. For first, second or higher order homogeneous ordinary differential equations, we form linear, quadratic or higher order polynomial equations incorporating the default signs of the terms in the differential equations. For non-homogeneous ordinary differential equations, we use the forcing functions. When the weights of the neural network are obtained by these systematic ways, there is no need to adjust all the parameters in the network, as postulated by previous researchers, in order to achieve convergence. All that is required is a little adjustment of the biases, and these are fixed to lie in a given interval and convergence to a solution with an acceptable minimum error is achieved. When the problem of obtaining the parameters is settled, it becomes easy to solve any first, second or higher order ordinary differential equation using the neural network model appropriate to it. We shall restrict this study to first and second order linear ordinary homogeneous differential equations. To compute the prediction error we use the squared error function defined as follows:

$E=0.5{\left({\psi }_d-{\psi }_p\right)}^2$ (21)

where ${\psi }_d$ represents the desired output and ${\psi }_p$ the predicted output. The same procedure can be applied to second and third order ODE. For the second order initial value problem (IVP):

${\psi }^{″}\left(x\right)=f\left(x,\psi ,{\psi }^{\prime }\left(x\right)\right),\psi \left(0\right)=A,{\psi }^{\prime }\left(0\right)=B$ (22)

The trial solution is written

${\psi }_t\left(x\right)=A+Bx+x^{2}N\left(x,p\right)$ (23)

where $A,B\in ℝ$ and $N\left(x_i,p\right)={\displaystyle \underset{j=1}{\overset{3}{\sum }}{v}_j\sigma \left(z_j\right)}$, and for two point BC: $\psi \left(0\right)=A,\psi \left(1\right)=B$ the trial solution in this case is written: ${\psi }_t\left(x\right)=A\left(1-x\right)+Bx+x\left(1-x\right)N\left(x,p\right)$.

Now, we first demonstrate the computational complexity involved in adjusting all parameters in order to update the weights and getting a close approximation to the desired result. Subsequently, we proceed to our main results and analysis that displays the ease of computation achieved by our novel method of adjusting only the biases. The former adjustment or weight updation is done using the backpropagation algorithm. Therefore, we need to train the network so we can apply the backpropagation algorithm. The basis function we shall apply in this work in training the neural network is the sigmoid activation function given by Equation (9). In a simple neural model as depicted in Figure 3, where there is no hidden layer and $x=\left(x_1,x_2\right),w=\left(w_1,w_2\right)$, we shall assume some figures for illustration. Let the training data be given as $x=\left(x_1,x_2\right)=\left(0.1,0.2\right)$ ; desired output $y_d=0.02$ ; initial weights $w=\left(w_1,w_2\right)=\left(0.4,0.1\right)$ ; the bias $b=1.78$ and predicted output $y_p$. The diagram below shows the neural network training model for the sample data we are considering. Now we proceed to train the network to get the predicted output.

5.1. Network Training

First we compute $z=x_1\cdot w_1+x_2\cdot w_2+b$, which is sum of products and bias, i.e.;

$z=x_1\left(w_1\right)+x_2\left(w_2\right)+b=0.1\left(0.4\right)+0.2\left(0.1\right)+1.78=1.84$ (24)

Figure 3. Schematic for $N\left(x,p\right)$.

Next we apply z, as input to the activation function, which in this case is the sigmoid activation function:

$\sigma \left(z\right)={\left(1+{\text{e}}^{-z}\right)}^{-1}={\left(1+{\text{e}}^{-1.84}\right)}^{-1}=0.863$. Hence the predicted output $y_p=0.863$.

We have seen that the predicted output does not correspond to the desired output, therefore we have to train the network to reduce the prediction error. To compute the prediction error, we use the squared error (21) function defined as follows. Considering the predicted output we calculated above, the prediction error is:

$E=0.5{\left(0.02-0.863\right)}^2=0.355$. (25)

We observe that the prediction error is huge, so we must attempt to minimize it. We noted previously that the weights determine how close a predicted output is to the desired output. Therefore to minimize the error we have to adjust the weights. This can be achieved using the formula;

$w_n=w_o+\eta \left(y_d-y_p\right)x$ (26)

where $w_n$ and $w_o$ represent new and old weights respectively. We update the weights using the following:

$w_o$ : current weight (1.78, 0.4, 0.1);

$\eta$ : network learning rate = 0.01;

$y_d$ : desired output = 0.02;

x: current input vector = (+1, 0.1, 0.2).

$\therefore w_n=\left[1.78,0.4,0.1\right]+0.01\left[0.02-0.863\right]\left[+1,0.1,0.2\right]=\left[1.772,0.399,0.098\right]$.

With this information we adjust the model and retrain the neural network to get;

$z=0.1\left(0.399\right)+0.2\left(0.098\right)+1.772=1.79$ $\Rightarrow \sigma \left(z\right)={\left(1+{\text{e}}^{-1.79}\right)}^{-1}=0.857$

$\therefore E=0.5{\left(0.02-0.857\right)}^2=0.350$.

5.2. Computation of the Gradient

The error computation not only involves the outputs but also the derivatives of the network output with respect to its inputs. So, it requires computing the gradient of the network derivatives with respect to its inputs. Let us now consider a multilayered perceptron with one input node, a hidden layer with m nodes, and one output unit. For the given inputs $x=\left(x_1,x_2,\cdots ,x_n\right)$, the output is given by

$N\left(x,p\right)={\displaystyle \underset{j=1}{\overset{m}{\sum }}{v}_j\sigma \left(z_j\right)}$ (27)

where $z_j={\displaystyle {\sum }_{i=1}^n{w}_{ji}{x}_i+u_j}$, $w_{ji}$ denotes the weight from input unit i to the hidden unit j, $v_j$ denotes weight from the hidden unit j to the output unit, $u_j$ denotes the biases, and $\sigma \left(z_j\right)$ is the sigmoid activation function. The derivatives of $N\left(x,p\right)$ with respect to input $x_i$ is

$\frac{{\partial }^{k}N}{\partial x_i^k}={\displaystyle \underset{j=1}{\overset{m}{\sum }}{v}_j{w}_{ji}^k{\sigma }_j^{\left(k\right)},}$ (28)

where $\sigma =\sigma \left(z_j\right)$ and ${\sigma }^{\left(k\right)}$ denotes the k^th order derivative of sigmoid function.

Let $N_{\theta }$ denote the derivative of the network with respect to its inputs and then we have the following relation

$N_{\theta }=D^{n}N={\displaystyle \sum v_i{p}_i{\sigma }_i^{\left(n\right)}}$ ; $p_j={\displaystyle \underset{k=1}{\overset{n}{\prod }}{w}_{jk}^{{\lambda }^k}}$, $k={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\lambda }_i}$ (29)

The derivative of $N_{\theta }$ with respect to other parameters may be obtained as

$\frac{\partial N_{\theta }}{\partial v_j}=p_j{\sigma }_j^{\left(k\right)}$, $\frac{\partial N_{\theta }}{\partial v_j}=v_j{p}_j{\sigma }_j^{\left(k+1\right)}$,

$\frac{\partial N_{\theta }}{\partial w_{ji}}=x_j{v}_j{p}_j{\sigma }_j^{\left(k+1\right)}+v_j{\lambda }_i{w}_{ji}^{{\lambda }_i-1}\left({\displaystyle \underset{k=1,k\ne 1}{\prod }{w}_{ji}^{{\lambda }_k}}\right){\sigma }_j^{\left(k\right)}$ (30)

Now after getting all the derivatives we can find out the gradient of error. Using general learning method for supervised training we can minimize the error to the desired accuracy. We illustrate the above using the first order ordinary differential equation below

${\psi }^{\prime }\left(x\right)=f\left(x,\psi \right),x\in \left[a,b\right]$ (31)

with initial condition $\psi \left(a\right)=A$. In this case, the ANN trial solution may be written as

${\psi }_t\left(x,p\right)=A+xN\left(x,p\right),$ (32)

where $N\left(x,p\right)$ is the neural output of the feed forward network with one input data x with parameters p. The trial solution ${\psi }_t\left(x,p\right)$ satisfies the initial condition. We differentiate the trial solution ${\psi }_t\left(x,p\right)$ to get

$\frac{\text{d}{\psi }_t\left(x,p\right)}{\text{d}x}=N\left(x,p\right)+x\frac{\text{d}N\left(x,p\right)}{\text{d}x},$ (33)

For evaluating the derivative term in the right hand side of (32), we use equations (6) and (26)–(31). The error function for this case may be formulated as

$E\left(p\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(\frac{\text{d}{w}_t\left(x_i,p\right)}{\text{d}{x}_i}-f\left(x_i,w_t\left(x_i,p\right)\right)\right)}^2}$ (34)

The weights from input to hidden are modified according to the following rule

$w_{ji}^{r+1}=w_{ji}^r-\eta \left(\frac{\partial E}{\partial w_{ji}^r}\right)$ (35)

where

$\frac{\partial E}{\partial w_{ji}^r}=\frac{\partial }{\partial w_{ji}^r}\left({\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(\frac{\text{d}{w}_t\left(x_i,p\right)}{\text{d}{x}_i}-f\left(x_i,w_t\left(x_i,p\right)\right)\right)}^2}\right)$ (36)

Here, 𝜂 is the learning rate and r is the iteration step. The weights from hidden to output layer may be updated in a similar formulation as done for input to hidden. Now going back to Equation (27), we recall that $z_j={\displaystyle {\sum }_{i=1}^n{w}_{ji}{x}_i+u_j}$

and $N\left(x,p\right)={\displaystyle \underset{j=1}{\overset{m}{\sum }}{v}_j\sigma \left(z_j\right)}$.This implies that;

$\begin{array}{c}\frac{\text{d}N\left(x,p\right)}{\text{d}x}=\frac{\text{d}}{\text{d}x}{\displaystyle \underset{j}{\overset{m}{\sum }}{v}_j\sigma \left(w_{ji}{x}_i+u_j\right)}={\displaystyle \underset{j}{\overset{m}{\sum }}{v}_j\frac{\text{d}}{\text{d}x}\sigma \left(w_{ji}{x}_i+u_j\right)}\\ ={\displaystyle \underset{j}{\overset{m}{\sum }}{v}_j\frac{\text{d}}{\text{d}x}\sigma \left(z_j\right)}={\displaystyle \underset{j}{\overset{m}{\sum }}{v}_j\frac{\text{d}}{\text{d}{z}_j}\sigma \left(z_j\right)\cdot }\frac{\text{d}{z}_j}{\text{d}x}\end{array}$ (37)

If the neural network model is a simple one as we saw in Section 3.3, then,

$N\left(x,p\right)=\sigma \left(z\right)=\sigma \left(x_1{w}_1+x_2{w}_2+u\right)$ $\Rightarrow \frac{\text{d}N}{\text{d}{x}_1}=\frac{\partial \sigma }{\partial z}\cdot \frac{\text{d}z}{\text{d}{x}_1};\frac{\text{d}N}{\text{d}{x}_2}=\frac{\partial \sigma }{\partial z}\cdot \frac{\text{d}z}{\text{d}{x}_2}$

Now let us consider a first order differential equation:

${\psi }^{\prime }\left(x\right)-\psi =0,x\in \left[0,1\right],\psi \left(0\right)=1$ (38)

with trial solution:

${\psi }_t\left(x,p\right)=A+xN\left(x,p\right)$ (39)

where x is the input to neural network model and p represents the parameters—weights and biases.

$w_t\left(0,p\right)=A+\left(0\right)N\left(0,p\right)=A=1$ $\Rightarrow {\psi }_t\left(x,p\right)=1+xN\left(x,p\right)$ (40)

To solve this problem using neural network (NN), we shall employ a NN architecture given in Figure 3. Now, ${\psi }_t\left(x_i,p\right)=1+x_{i}N\left(x_i,p\right)$, with $N\left(x_i,p\right)={\displaystyle {\sum }_{j=1}^3{v}_j\sigma \left(z_j\right)}$, and

${\displaystyle \underset{j=1}{\overset{3}{\sum }}{v}_j\sigma \left(z_j\right)}=v_1\sigma \left(z_1\right)+v_2\sigma \left(z_2\right)+v_3\sigma (\; z\; 3\; )$

where $z_1=x_1{w}_{11}+u_1$, $z_2=x_1{w}_{12}+u_2$ and $z_3=x_1{w}_{13}+u_3$.

If the neural network model is not able to predict correctly the solution of the differential equation with the given initial parameters—weights and biases, we need to find the prediction error given by

$E\left(p\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(\frac{\text{d}{w}_t\left(x_i,p\right)}{\text{d}{x}_i}-f\left(x_i,w_t\left(x_i,p\right)\right)\right)}^2}$. (41)

If the prediction error does not satisfy an acceptable threshold, then the parameters need to be adjusted using the equation,

$w_{ji}^{r+1}=w_{ji}^r-\eta \left(\frac{\partial E}{\partial w_{ji}^r}\right)$, where $\frac{\partial E}{\partial w_{ji}^r}=\frac{\partial }{\partial w_{ji}^r}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(\frac{\text{d}{w}_t\left(x_i,p\right)}{\text{d}{x}_i}-f\left(x_i,w_t\left(x_i,p\right)\right)\right)}^2}$ (42)

Recall that: ${\psi }_t\left(x_i,p\right)=1+x_{i}N\left(x_i,p\right)$, and

$\frac{\text{d}{w}_t}{\text{d}{x}_i}\left(x_i,p\right)=\frac{\text{d}}{\text{d}{x}_i}\left(1+x_{i}N\left(x_i,p\right)\right)=N\left(x_i,p\right)+x_i\frac{\text{d}}{\text{d}{x}_i}N\left(x_i,p\right)$. (43)

$\begin{array}{l}\therefore \frac{\text{d}}{\text{d}{x}_1}N\left(x_1,p\right)=\frac{\text{d}}{\text{d}{x}_1}{\displaystyle \underset{j=1}{\overset{3}{\sum }}{v}_j\sigma \left(z_j\right)}=v_j\frac{\text{d}}{\text{d}{x}_1}{\displaystyle \underset{j=1}{\overset{3}{\sum }}\sigma \left(z_j\right)}\\ =v_1\frac{\text{d}\sigma }{\text{d}{z}_1}\frac{\text{d}{z}_1}{\text{d}{x}_1}+v_2\frac{\text{d}\sigma }{\text{d}{z}_2}\frac{\text{d}{z}_2}{\text{d}{x}_1}+v_3\frac{\text{d}\sigma }{\text{d}{z}_3}\frac{\text{d}{z}_3}{\text{d}{x}_1}\end{array}$

$\frac{\text{d}}{\text{d}{x}_1}N\left(x_1,p\right)=v_1{\sigma }^{\prime }\left(z_1\right)w_{11}+v_2{\sigma }^{\prime }\left(z_2\right)w_{12}+v_3{\sigma }^{\prime }\left(z_3\right)w_{13}$ (44)

Putting Equations (42) and (43) into Equation (40) gives

$E\left(p\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(\left(N\left(x_i,p\right)+x_i\frac{\text{d}}{\text{d}{x}_i}N\left(x_i,p\right)\right)-\left(1+xN\left(x,p\right)\right)\right)}^2}$

$\therefore E\left(p\right)={\left({\displaystyle \underset{j=1}{\overset{3}{\sum }}{v}_j\sigma \left(z_j\right)+x_1\left(\frac{\text{d}}{\text{d}{x}_1}N\left(x_i,p\right)\right)}-1-x_1\left({\displaystyle \underset{j=1}{\overset{3}{\sum }}{v}_j\sigma \left(z_j\right)}\right)\right)}^2$ (45)

$\begin{array}{l}\Rightarrow E\left(p\right)=(v_1\sigma \left(z_1\right)+v_2\sigma \left(z_2\right)+v_3\sigma \left(z_3\right)\\ +x_1\left(v_1{\sigma }^{\prime }\left(z_1\right)w_{11}+v+v_2{\sigma }^{\prime }\left(z_2\right)w_{12}+v_3{\sigma }^{\prime }\left(z_3\right)w_{13}\right)\\ -1-{x_1\left(v_1\sigma \left(z_1\right)+v_2\sigma \left(z_2\right)+v_3\sigma \left(z_3\right)\right))}^2\end{array}$ (46)

We noted earlier that when the neural network is not able to predict acceptable solution, that is, solution with minimum error, the weights and biases need to be adjusted. This involves complex derivatives with multivariable chain rule [11] [12]. From the ongoing, we need to compute the derivative of Equation (46) with respect to the weights and biases. That is:

$\begin{array}{c}\frac{\partial E}{\partial w_{ji}^r}=\frac{\partial }{\partial w_{ji}^r}(v_1\sigma \left(z_1\right)+v_2\sigma \left(z_2\right)+v_3\sigma \left(z_3\right)\\ +x_1\left(v_1{\sigma }^{\prime }\left(z_1\right)w_{11}+v+v_2{\sigma }^{\prime }\left(z_2\right)w_{12}+v_3{\sigma }^{\prime }\left(z_3\right)w_{13}\right)\\ -1-{x_1\left(v_1\sigma \left(z_1\right)+v_2\sigma \left(z_2\right)+v_3\sigma \left(z_3\right)\right))}^2\end{array}$ (47)

Similarly, to update the weights from the hidden layer to output layer, we compute $\partial E/\partial v_j,j=1,2,3$. Finally, we update the biases by computing $\partial E/\partial u_j,j=1,2,3$. Equation (47) together with $\partial E/\partial v_j$ and $\partial E/\partial u_j$ are used to update the weights and biases. The superscript (r) denotes the r^th iteration.

It is important to note that the foregoing is necessary in order to achieve the number of iterations required. Sometimes, it may be necessary to do up to 30 or more iterations for the solution to converge within an acceptable threshold. It is on the basis of this complex derivatives involved in solving ODE with neural network, especially of higher order, and the many iterations required for convergence that motivated our search for a more efficient and accurate way of solving the given problem, but avoiding the inherent computational complexity.

5.3. Results

We begin with a couple of examples on first and second order differential equations.

5.4. Example

Consider the initial value problem;

$y^{\prime }-\alpha y=0;y\left(0\right)=1,x\in \left[0,1\right],\alpha \equiv 1.$ (48)

This equation has been solved by Mall and Chakraverty [13]. As discussed in the previous section, the trial solution is given by;

$y_t\left(x\right)=A+xN\left(x,p\right)$

Applying the initial condition gives $A=1$, therefore $y_t\left(x\right)=1+xN\left(x,p\right)$.

To obtain the weights from input to hidden layer, it is natural to assume $y_a\left(x\right)={\text{e}}^x$. We approximate this function using regression in SPSS. We shall train the network for 10 equidistant points in [0, 1], and then employ excel spreadsheet to find values of $y_a\left(x\right)={\text{e}}^x$ at all the points. This leads us to the following data in Table 3.

We then perform a curve-fit polynomial regression built into IBM SPSS 23 [14], for the function $y_a\left(x\right)={\text{e}}^x$. The output is displayed below in Table 4.

Table 4 displays the linear, quadratic and cubic regression of $y_a\left(x\right)={\text{e}}^x$. The quadratic and cubic curves show perfect goodness of fit, $R^2=1$. Using the cubic curve, we pick our weights from input layer to hidden layer as:

Table 3. Values of $\left(x,y_a\left(x\right)={\text{e}}^x\right)$ for problem (47).

Table 4. Model summary and parameter estimates.

$w_1{}_1=1.016,w_{12}=0.423,w_{13}=0.279$. Now to compute the weights from hidden layer to the output layer, we find a function $\vartheta \left(x\right)$ such that $v={\varphi }^{-1}f$ ; $v,f$ and $\varphi$ are as defined in Section 3. In particular, $f\left(x\right)={\left(\vartheta \left(x_1\right)\text{,}\vartheta \left(x_2\right)\text{,}\vartheta \left(x_3\right)\right)}^{\text{T}}$. We now form a linear function based on the default sign of the differential equation, i.e. $\vartheta \left(x\right)=ax-b$, where a is the coefficient of the derivative of y and b is the coefficient of y Thus;

$\vartheta \left(x\right)=x+1,f\left(x\right)={\left(\vartheta \left(x_1\right)\text{,}\vartheta \left(x_2\right)\text{,}\vartheta \left(x_3\right)\right)}^{\text{T}}={\left(1.1,\text{1}\text{.2,1}\text{.3}\right)}^{\text{T}}$

The neural architecture for the neural network is shown in Figure 2, so we let $N=3$. We take $x={\left(0.1,0.2,0.3\right)}^{\text{T}}$ and $f\left(x\right)={\left(1.1,1.2,1.3\right)}^{\text{T}}$. It then follows that

$v={\varphi }^{-1}f$, $\Rightarrow \left[\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right]={\left[\begin{array}{ccc}{\phi }_1\left(x_1\right)& {\phi }_2\left(x_1\right)& {\phi }_3\left(x_1\right)\\ {\phi }_1\left(x_2\right)& {\phi }_2\left(x_2\right)& {\phi }_3\left(x_2\right)\\ {\phi }_1\left(x_3\right)& {\phi }_2\left(x_3\right)& {\phi }_3\left(x_3\right)\end{array}\right]}^{-1}\left[\begin{array}{c}{\vartheta }_1\\ {\vartheta }_2\\ {\vartheta }_3\end{array}\right]$ (49)

where

${\phi }_i\left(x_j\right)=\exp \left(-\frac{{\left(\left|x_i-x_j\right|\right)}^2}{2{\sigma }^2}\right),i=1,2,3;j=1,2,3.$ (50)

Substituting the given values of the vectors x and f, we obtain the weights from the hidden layer to the output layer,

$\left[\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right]={\left[\begin{array}{ccc}\text{1}& \text{0}\text{.94}& \text{0}\text{.78}\\ \text{0}\text{.94}& \text{1}& \text{0}\text{.94}\\ \text{0}\text{.78}& \text{0}\text{.94}& \text{1}\end{array}\right]}^{-1}\left[\begin{array}{c}\text{1}\text{.1}\\ \text{1}\text{.2}\\ \text{1}\text{.3}\end{array}\right]=\left[\begin{array}{c}\text{5}\text{.17}\\ -\text{9}\text{.375}\\ \text{6}\text{.08}\end{array}\right]$, (51)

Therefore the weights from the hidden layer to the output layer are; $v_1=5.17,v_2=-9.375,v_3=6.08$.

The biases are fixed between –1 and 1. We now train the network with the available parameters using MathCAD 14 software [15] as follows:

$\begin{array}{l}{w}_1:=1.016w_2:=0.423w_3:=0.279x:=1\\ v_1:=5.17v_2:=6.08u_1:=1u_2:=0.2251u_3:=-0.1\\ z_1:=w_1\cdot x+u_1=2.016z_2:=w_2\cdot x+u_2=0.6481z_3:=w_3\cdot x+u_3=0.179\\ \sigma \left(z_1\right):={\left[1+\exp \left(z_1\right)\right]}^{-1}=0.882467,\sigma \left(z_2\right):={\left[1+\exp \left(z_2\right)\right]}^{-1}=0.656582,\\ \sigma \left(z_3\right):={\left[1+\exp \left(z_3\right)\right]}^{-1}=0.544631\\ N:=v_1\cdot \sigma \left(z_1\right)+v_2\cdot \sigma \left(z_2\right)+v_3\cdot \sigma \left(z_3\right)=1.718251\\ y_p\left(x\right):=1+x\cdot N=2.718251,y_d\left(x\right):={\text{e}}^x=2.718282\\ E:=0.5\cdot {\left(y_d\left(x\right)-y_p\left(x\right)\right)}^2=4.707964\times {10}^{-10}\end{array}$.

Here y_d and y_p are respectively the desired output (exact solution) and the predicted output (trial solution). From the indicated error value, this is an acceptable accuracy. We compare our results with the neural network results obtained by Otadi and Mosleh [16] and find them to be in reasonable agreement. This is depicted in Table 5, as well as the graphical profile in Figure 4 below.

The perfect accuracy is evident in the graphical profile depicted in Figure 4.

5.5. Remark

In what follows, we consider a non-homogeneous second order linear differential equation. It is important to recall that for any second order non-homogeneous differential equation of the form $y^{″}\left(x\right)+a\left(x\right)y^{\prime }+b\left(x\right)y=f\left(x\right)$, the non-homogeneous term $f\left(x\right)$ is termed the forcing function. In this section, we shall employ the forcing function to compute the weights from hidden layer to the output layer. This is made clear in the following example.

5.6. Example

Consider the initial value problem;

$y^{″}-4y=24\cos \left(2x\right);y\left(0\right)=3,y^{\prime }\left(0\right)=4,x\in \left[0,1\right]$.

The trial solution is $y_t\left(x\right)=A+Bx+x^{2}N\left(x,p\right)$. Applying the initial condition gives $A=3,B=4$.

Therefore, $y_t\left(x\right)=3+4x+x^{2}N\left(x,p\right)$, $y_a\left(x\right)={\text{e}}^{-2x}+{\text{e}}^{-3x}$.

We use excel spreadsheet to find values of $y_a\left(x\right)$ at all the x points, as displayed in Table 6.

Using regression and SPSS model we find weights from input layer to hidden layer. From Table 7 and using the cubic curve fit with R² = 1, we pick our weights from input layer to hidden layer as:

$w_1{}_1=0.488,w_{12}=1.697,w_{13}=3.338$.

Now to compute the weights from hidden layer to the output layer, we use the function:

Table 5. Comparison of the results.

Figure 4. Plot of Y exact and Y predicted for Example 1.

$\begin{array}{l}\vartheta \left(x\right)=24\cos \left(2x\right),\\ f\left(x\right)={\left(\vartheta \left(x_1\right)\text{,}\vartheta \left(x_2\right)\text{,}\vartheta \left(x_3\right)\right)}^{\text{T}}={\left(23.999854,23.999415\text{,23}\text{.998684}\right)}^{\text{T}}\end{array}$

With $x={\left(0.1,0.2,0.3\right)}^{\text{T}}$. Hence, the weights from the hidden layer to the output layer given by $v={\varphi }^{-1}f$ are;

$\left[\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right]={\left[\begin{array}{ccc}\text{1}& \text{0}\text{.94}& \text{0}\text{.78}\\ \text{0}\text{.94}& \text{1}& \text{0}\text{.94}\\ \text{0}\text{.78}& \text{0}\text{.94}& \text{1}\end{array}\right]}^{-1}\left[\begin{array}{c}23.999854\\ 23.999415\\ 23.998684\end{array}\right]$ $\Rightarrow \left[\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right]=\left[\begin{array}{c}112.489\\ -187.474\\ 112.483\end{array}\right]$

The weights from the hidden layer to the output layer are; $v_1=112.489,v_2=-187.474,v_3=112.483$.

The biases are fixed between −1 and 1. We now train the network with the available parameters using our MathCAD 14 algorithm as follows:

$\begin{array}{l}{\text{w}}_{\text{1}}:=0.488{\text{w}}_{\text{2}}:=1.697{\text{w}}_3:=3.338\text{x}:=1\\ {\text{v}}_1:=112.489{\text{v}}_{\text{2}}:=-187.474{\text{v}}_3:=112.483{\text{u}}_1:=1{\text{u}}_2:=1{\text{u}}_3:=-0.1691\\ {\text{z}}_1:={\text{w}}_{\text{1}}\cdot \text{x}+{\text{u}}_1=1.488{\text{z}}_2:={\text{w}}_2\cdot \text{x}+{\text{u}}_2=2.697{\text{z}}_3:={\text{w}}_3\cdot \text{x}+{\text{u}}_3=3.1689\\ \text{σ}\left({\text{z}}_1\right):={\left[1+\exp \left(-{\text{z}}_1\right)\right]}^{-1}=0.815778,\text{σ}\left({\text{z}}_2\right):={\left[1+\exp \left(-{\text{z}}_2\right)\right]}^{-1}=0.936849,\\ \text{σ}\left({\text{z}}_3\right):={\left[1+\exp \left(-{\text{z}}_3\right)\right]}^{-1}=0.959647\\ N:={\text{v}}_1\cdot \text{σ}\left({\text{z}}_1\right)+{\text{v}}_2\cdot \text{σ}\left({\text{z}}_2\right)+{\text{v}}_3\cdot \text{σ}\left({\text{z}}_3\right)=24.075112\\ {\text{y}}_{\text{p}}\left(\text{x}\right):=3+4\cdot \text{x}+{\text{x}}^2\cdot \text{N}=31.075112,\\ {\text{y}}_{\text{d}}\left(\text{x}\right):=4\cdot {\text{e}}^{2\cdot \text{x}}+2\cdot {\text{e}}^{-2\cdot \text{x}}-3\cdot \cos \left(2\cdot \text{x}\right)=31.075335\\ \text{E}:=0.5\cdot {\left({\text{y}}_{\text{d}}\left(\text{x}\right)-{\text{y}}_{\text{p}}\left(\text{x}\right)\right)}^2=2.502862\times {10}^{-8}\end{array}$.

We compare the exact and approximate solution in Table 8. The accuracy is clearly depicted graphically in Figure 5.

Table 6. Values of $\left(x,y_a\left(x\right)={\text{e}}^{2x}+{\text{e}}^{-2x}\right)$.

Table 7. Model summary and parameter estimates.

Table 8. Comparison of the results.

Figure 5. Plot of Y exact and Y predicted for Example 2.

6. Conclusion

In this paper, we have presented a novel approach for solving first and second order linear ordinary differential equations with constant coefficients. Specifically, we employ a feed-forward Multilayer Perceptron Neural Network (MLPNN), but avoid the standard back-propagation algorithm for updating the intrinsic weights. This greatly reduces the computational complexity of the given problem. Our results are validated by the near perfect approximations achieved in comparison with the exact solutions, as well as demonstrating the function approximation capabilities of ANN. This then proves the efficiency of our neural network procedure. We employed Excel spreadsheet, IBM SPSS 23, and MathCAD 14 algorithm to achieve this task.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	McCulloch, W.S. and Pitts W. (1943) A Logical Calculus of the Ideas Immanent in Nervous Activity. The Bulletin of Mathematical Biophysics, 5, 115-133. https://doi.org/10.1007/BF02478259
[2]	Neumann, J.V. (1951) The General and Logical Theory of Automata. Wiley, New York.
[3]	Graupe, D. (2007) Principles of Artificial Neural Networks. Vol. 6, 2nd Edition, World Scientific Publishing Co. Pte. Ltd., Singapore.
[4]	Rumelhart, D.E. and McClelland, J.L. (1986) Parallel Distributed Processing, Explorations in the Microstructure of Cognition I and II. MIT Press, Cambridge. https://doi.org/10.7551/mitpress/5236.001.0001
[5]	Werbos, P.J. (1974) Beyond Recognition, New Tools for Prediction and Analysis in the Behavioural Sciences. Ph.D. Thesis, Harvard University, Cambridge.
[6]	Lagaris, I.E., Likas, A.C. and Fotiadis, D.I. (1997) Artificial Neural Network for Solving Ordinary and Partial Differential Equations. arXiv: physics/9705023v1.
[7]	Cybenco, G. (1989) Approximation by Superposition of a Sigmoidal Function. Mathematics of Control, Signals and Systems, 2, 303-314. https://doi.org/10.1007/BF02551274
[8]	Hornic, K., Stinchcombe, M. and White, H. (1989) Multilayer Feedforward Networks Are Universal Approximators. Neural Networks, 2, 359-366. https://doi.org/10.1016/0893-6080(89)90020-8
[9]	Lee, H. and Kang, I.S. (1990) Neural Algorithms for Solving Differential Equations. Journal of Computational Physics, 91, 110-131. https://doi.org/10.1016/0021-9991(90)90007-N
[10]	Majidzadeh, K. (2011) Inverse Problem with Respect to Domain and Artificial Neural Network Algorithm for the Solution. Mathematical Problems in Engineering, 2011, Article ID: 145608, 16 p. https://doi.org/10.1155/2011/145608
[11]	Chen, R.T.Q., Rubanova, Y., Bettencourt, J. and Duvenaud, D. (2018) Neural Ordinary Differential Equations. arXiv: 1806.07366v1.
[12]	Okereke, R.N. (2019) A New Perspective to the Solution of Ordinary Differential Equations using Artificial Neural Networks. Ph.D Dissertation, Mathematics Department, Michael Okpara University of Agriculture, Umudike.
[13]	Mall, S. and Chakraverty, S. (2013) Comparison of Artificial Neural Network Architecture in Solving Ordinary Differential Equations. Advances in Artificial Neural Systems, 2013, Article ID: 181895. https://doi.org/10.1155/2013/181895
[14]	IBM (2015) IBM SPSS Statistics 23 http://www.ibm.com
[15]	PTC (Parametric Technology Corporation) (2007) Mathcad Version 14. http://[e-mail]
[16]	Otadi, M. and Mosleh, M. (2011) Numerical Solution of Quadratic Riccati Differential Equations by Neural Network. Mathematical Sciences, 5, 249-257.