\documentclass[12pt]{article}
\usepackage[margin=1in]{geometry}
\usepackage{amsmath,amsthm,amsfonts,amssymb}
\usepackage{graphicx}
\usepackage{listings}
\usepackage{xcolor}
\usepackage{graphicx}
\usepackage{booktabs}


% Define styles for Python code
\lstset{
    language=Python,
    frame=single,
    basicstyle=\ttfamily\footnotesize,
    keywordstyle=\color{blue},
    commentstyle=\color{gray},
    stringstyle=\color{orange},
    numbers=left,
    numberstyle=\tiny,
    stepnumber=1,
    numbersep=5pt,
    showstringspaces=false,
    breaklines=true,
    breakatwhitespace=false,
    tabsize=4,
    backgroundcolor=\color{lightgray!20}
}
\usepackage[english]{babel}

\title{Solving Ordinary Differential Equations using Neural Networks}
\author{Your Name}

\begin{document}

\maketitle

\section{Introduction}

Differential equations are widely used in various disciplines to model real-world phenomena, but finding analytical solutions can be challenging. When analytical solutions are not possible, numerical methods are often employed, though they tend to be computationally expensive. Recently, neural networks, which are known as universal function approximators, have gained interest for solving ordinary differential equations (ODEs) due to their high accuracy and fast convergence.

To solve ODEs using neural networks, a trial solution is constructed to meet the boundary and initial conditions. This solution typically takes the form \( y_{trial}(t) = A(t) + F(t, N(t, \vec{p})) \), where \( N(t, \vec{p}) \) is the neural network output, and \( A(t) \) satisfies the boundary conditions. The network approximates the solution by minimizing the error, often using mean squared error for adjustment.

Neural networks provide several advantages over traditional methods. They can significantly reduce computation time and offer differentiable solutions, unlike methods that rely on domain discretization. While neural network methods are promising, their implementation details, such as network architecture and training parameters, can greatly impact results.

\subsection{Background}
Neural networks are inspired by the biological nervous system, where neurons transmit signals once a certain threshold is reached. Artificial neurons abstract this concept and are designed to learn and adapt based on feedback. The first artificial neuron model was introduced by McCulloch and Pitts in 1943, which demonstrated that logical functions can be represented through networks of neurons using inputs, weights, outputs, and transfer functions. The model used a binary step function to simulate activation, where neurons are either fully activated or not, depending on the input.

\begin{figure}[h]
    \centering
    \includegraphics[width=0.5\textwidth]{neuron.png}
    \caption{A simple illustration of a biological neuron.}
    \label{fig:neuron}
\end{figure}

Mathematically, a neuron’s output can be described by the weighted sum of its inputs, calculated as follows:

\[
S_n = \sum_{i=1}^{n} w_i t_i, \quad n \in \mathbb{N},
\]
where \( w_i \) are the weights and \( t_i \) are the inputs. These weighted sums are then passed through an activation function. In more complex networks, such as the multilayer perceptron (MLP), the network consists of an input layer, one or more hidden layers, and an output layer. This type of feedforward neural network has no cycles and is widely used for its versatility in various applications.

For solving non-trivial problems, nonlinear activation functions are necessary. Two commonly used activation functions in artificial neural networks are the sigmoid and rectified linear unit (ReLU) functions. The sigmoid function is defined as:

\[
\sigma(t) = \frac{1}{1 + e^{-t}}, \quad \sigma: \mathbb{R} \to \mathbb{R},
\]

and the ReLU function is given by:

\[
f(t) = \max(0, t), \quad f: \mathbb{R} \to \mathbb{R},
\]

where \( t \in \mathbb{R} \).

In MLPs, backpropagation is the standard algorithm for learning. Backpropagation adjusts the weights based on the gradient of the loss function, typically computed using gradient descent or stochastic gradient descent. These methods aim to find a local minimum of the loss function, thus minimizing the error in the network’s output.

The choice of activation function and network architecture plays a critical role in solving different types of problems. The multilayer perceptron’s structure and learning capabilities, combined with effective activation functions, allow it to solve complex differential equations with high accuracy.


\section{Methodology}

In this thesis, the focus is on solving second-order ordinary differential equations (ODEs) of the form:

\begin{equation}
G(t, y(t), y'(t), y''(t)) = 0, \quad t \in D \subset \mathbb{R},
\end{equation}
where \( t \) is the independent variable, and \( y(t) \) is the solution. While the method is generalizable to higher-order differential equations, the scope is limited to second-order ODEs due to their numerous applications in science and engineering.

\subsection{Trial Solution Method}

The core method involves constructing a trial solution that satisfies both the differential equation and its boundary conditions. The trial solution is represented as:

\begin{equation}
y_{\text{trial}}(t) = A(t) + F(t, N(t, \vec{p})),
\end{equation}
where \( A(t) \) satisfies the boundary conditions, and \( N(t, \vec{p}) \) is the output of a neural network with parameters \( \vec{p} \), designed such that \( F(t, N(t, \vec{p})) \) does not affect the boundary conditions.
For first-order ODEs of the form:

\begin{equation}
\frac{dy}{dt} = f(t, y(t)),
\end{equation}
a trial solution might be:

\begin{equation}
y_{\text{trial}}(t) = A + tN(t, \vec{p}),
\end{equation}
ensuring the boundary conditions are satisfied. The neural network adjusts its weights to minimize the error between the trial solution and the actual solution of the ODE. The error function used is the Mean Squared Error (MSE), defined as:

\begin{equation}
E[\vec{p}] = \frac{1}{n} \sum_{i=1}^{n} \left( \frac{dy_{\text{trial}}(t_i)}{dt} - f(t_i, y_{\text{trial}}(t_i)) \right)^2,
\end{equation}
where \( t_i \in D \) is a discretization of the domain.
\vspace{12pt}

For second-order ODEs, the trial solution and error function are modified to account for higher derivatives. For example, with the boundary conditions \( y(0) = A \) and \( y'(0) = A' \), the trial solution becomes:

\begin{equation}
y_{\text{trial}}(t) = A + A't + t^2 N(t, \vec{p}),
\end{equation}
and the corresponding error function is:

\begin{equation}
E[\vec{p}] = \frac{1}{n} \sum_{i=1}^{n} \left( \frac{d^2 y_{\text{trial}}(t_i)}{dt^2} - f\left( t_i, y_{\text{trial}}(t_i), \frac{dy_{\text{trial}}(t_i)}{dt} \right) \right)^2.
\end{equation}

\subsection{Application to Systems of ODEs}

The methodology is extended to systems of ODEs, where a trial solution is constructed for each equation in the system. For a system of \( K \) first-order ODEs:

\begin{equation}
\frac{dy_i}{dt} = f_i(t, y_1, y_2, \dots, y_K),
\end{equation}with \( y_i(0) = A_i \), the trial solution for each equation is given by:
\begin{equation}
y_{\text{trial},i}(t) = A_i + t N_i(t, \vec{p_i}).
\end{equation}
The error function is generalized for the system as:
\begin{equation}
E[\vec{p}] = \frac{1}{n} \sum_{k=1}^{K} \sum_{i=1}^{n} \left( \frac{dy_{\text{trial},k}(t_i)}{dt} - f_k\left( t_i, y_{\text{trial},1}, \dots, y_{\text{trial},K} \right) \right)^2,
\end{equation}
where \( K \) is the number of equations in the system. The neural network is trained to minimize this error by adjusting its parameters. 
This method leverages neural networks to approximate solutions to ODEs and systems of ODEs, ensuring accuracy and consistency with boundary and initial conditions.
\subsection{Illustrative Example}

To demonstrate the methodology, we consider a second-order ordinary differential equation given by:

\[
\frac{d^2 y}{dt^2} + y(t) = \tan(t), \quad t > 0,
\]
with initial conditions:

\[
y(0) = 1, \quad \frac{dy}{dt}(0) = 1.
\]
The exact analytical solution for this equation is:

\[
y(t) = \cos(t) + 2 \sin(t) - \cos(t) \ln \left( \sec(t) + \tan(t) \right).
\]
In order to construct a neural network-based solution that satisfies these boundary conditions, we define a trial solution as:
\[
y_{\text{trial}}(t) = 1 + t + t^2 N(t, \vec{p}),
\]
where \( N(t, \vec{p}) \) represents the neural network’s output and \( \vec{p} \) refers to the parameters (weights and biases) of the neural network. The term \( 1 + t \) is selected to satisfy the initial conditions, while \( t^2 N(t, \vec{p}) \) ensures that the boundary conditions remain unaffected as \( t \rightarrow 0 \).

\subsubsection{Code Implementation}

To solve the differential equation using neural networks, we define and train a neural network using the following code snippets. The code is written in Python using the \texttt{autograd} library, which allows automatic differentiation. Below, we break down the core functions:

\begin{lstlisting}[caption={Sigmoid Activation Function}, label={lst:sigmoid}]
def sigmoid(z):
    """
    Sigmoid activation function.
    :param z: The input value
    :return: The output of the activation function
    """
    return 1 / (1 + np.exp(-z))
\end{lstlisting}

The \texttt{sigmoid} function in Listing~\ref{lst:sigmoid} is used as the activation function for the neurons. This non-linear transformation helps the network approximate complex functions. The next step is setting up the neural network:

\begin{lstlisting}[caption={Neural Network Structure}, label={lst:neural_network}]
def neural_network(parameters, x_in):
    """
    Sets up the neural network and calculates the output.
    :param parameters: The weights and biases of the network
    :param x_in: The input values
    :return: The output of the neural network
    """
    # Number of layers
    n_hidden = len(parameters) - 1
    x_prev = np.atleast_2d(x_in)
    
    for i in range(n_hidden):
        # Finds the weights for each hidden layer.
        w_hidden = parameters[i]
        x_prev = np.vstack((np.ones((1, x_prev.shape[1])), x_prev))
        # Hidden layer output
        z_hidden = np.dot(w_hidden, x_prev)
        x_prev = sigmoid(z_hidden)

    # Output layer
    w_output = parameters[-1]
    x_prev = np.vstack((np.ones((1, x_prev.shape[1])), x_prev))
    z_output = np.dot(w_output, x_prev)
    return z_output
\end{lstlisting}

The function in Listing~\ref{lst:neural_network} defines the structure and forward pass of the neural network. It takes the input values and a set of network parameters, passing them through each layer until reaching the output.

The training process is handled by the \texttt{neural\_network\_solve\_ode} function. This function initializes the network weights, defines a cost function, and iteratively updates the weights using gradient descent.

\begin{lstlisting}[caption={Neural Network Training}, label={lst:nn_training}]
def neural_network_solve_ode(x_in, num_neurons, num_iter, learning_rate):
    """
    Sets up the parameters for the neural network and solves the given ODE.
    :param x_in: The input values
    :param num_neurons: The number of neurons in the network
    :param num_iter: The number of iterations used for training
    :param learning_rate: The rate at which the network learns
    :return: The network parameters
    """
    n_hidden = len(num_neurons)
    parameters = [npr.randn(num_neurons[0], 1 + 1)]  # Initial parameters
    
    for j in range(1, n_hidden):
        parameters.append(npr.randn(num_neurons[j], num_neurons[j - 1] + 1))
    
    parameters.append(npr.randn(1, num_neurons[-1] + 1))

    def cost_function(parameters, x_in):
        """
        Returns the mean square error of the approximation.
        :param parameters: The weights and biases of the network
        :param x_in: The input values
        :return: The mean square error
        """
        g_t = g_trial_solution(x_in, parameters)
        d2_g_t = elementwise_grad(elementwise_grad(g_trial_solution, 0))(x_in, parameters)
        right_side = np.tan(x_in) - g_t

        err_sqr = (d2_g_t - right_side) ** 2
        cost_sum = np.sum(err_sqr)
        return cost_sum / np.size(err_sqr)

    cost_function_grad = grad(cost_function, 0)
    
    for i in range(num_iter):
        cost_grad = cost_function_grad(parameters, x_in)
        for j in range(n_hidden + 1):
            parameters[j] -= learning_rate * cost_grad[j]

    return parameters
\end{lstlisting}

The code iteratively minimizes the error by updating network parameters. This function is then used to solve the ODE.

\begin{lstlisting}[caption={Main Script for Solving the ODE}, label={lst:main_script}]
if __name__ == '__main__':
    npr.seed(4155)
    nx = 10
    x = np.linspace(0, 1, nx)
    num_hidden_neurons = [200, 100]
    num_iter = 1000
    learning_rate = 1e-3

    t1_start = process_time()
    P = neural_network_solve_ode(x, num_hidden_neurons, num_iter, learning_rate)
    t1_stop = process_time()

    g_dnn_ag = g_trial_solution(x, P)
    g_analytical = g_exact(x)

    plt.figure(figsize=(10, 10))
    plt.plot(x, g_analytical, marker="o")
    plt.plot(x, g_dnn_ag[0, :])
    plt.legend(['Exact Solution', 'Neural Network Approximation'])
    plt.xlabel('t')
    plt.ylabel('y(t)')
    plt.title('Approximate solution to ODE')
    plt.show()
\end{lstlisting}

The complete implementation in Listings~\ref{lst:sigmoid} to \ref{lst:main_script} demonstrates solving a second-order ODE using a neural network approach.


\section{Results}



\subsection{Performance Metrics}

To quantitatively evaluate the performance of the neural network in solving the ODE, the table below summarizes key metrics, including the Mean Square Error (MSE), time elapsed during training, and the maximum absolute difference between the exact and predicted solutions.

\begin{table}[ht!] % Use [h!] to enforce the table position
    \centering
    \begin{tabular}{|c|c|}
        \hline
        \textbf{Metric} & \textbf{Value} \\ \hline
        Mean Square Error (MSE) & 5.67587e-05 \\ % Update this with your computed MSE
        Time Elapsed (s) & 2.52915 seconds \\ % Update this with the time your code took to run
        Maximum Absolute Difference & 1.21569e-02 \\ % Update this with the maximum absolute error
        \hline
    \end{tabular}
    \caption{Performance Metrics of the Neural Network-Based ODE Solution}
    \label{table:performance_metrics}
\end{table}



\subsection{Graphical Representation of the Predicted Solution}

The figure below illustrates the comparison between the neural network-predicted solution and the exact analytical solution for the given ODE. The predicted solution closely follows the true solution, demonstrating the neural network's ability to approximate differential equations accurately.

\begin{figure}[ht!] % Use [ht!] to control the placement more strictly
    \centering
    \includegraphics[width=0.7\linewidth]{ode_results.png} % Replace 'ode_results.png' with the path to your image file
    \caption{Comparison between Exact and Predicted Solution}
    \label{fig:predicted_solution}
\end{figure}
\section{Conclusion}

This study explored the application of neural networks to solve ordinary differential equations (ODEs) and systems using a trial solution approach. The neural network method ensures compliance with initial and boundary conditions while leveraging the universal approximation capabilities of neural networks, producing a fully differentiable solution.

Although the use of neural networks for solving differential equations is not entirely new, it remains a growing research area with promising potential. Compared to traditional methods that often rely on discretization and can be time-consuming, neural networks offer an efficient alternative for modeling complex real-world systems where exact solutions are difficult to obtain. As the field of machine learning continues to evolve, these approaches are expected to become more effective and widely applicable for solving intricate differential equations in various domains.


\end{document}

Why are images not showing and format incorrectg