# GENERALIZED KURAMOTO-SIVASHINSKY EQUATION SOLVING BY B-SPLINES BASED NEURAL NETWORK

Adetola Jamal, Charbel Zeus Jupiter Mamlankou, Brice Adjalla

► To cite this version:

Adetola Jamal, Charbel Zeus Jupiter Mamlankou, Brice Adjalla. GENERALIZED KURAMOTO-SIVASHINSKY EQUATION SOLVING BY B-SPLINES BASED NEURAL NETWORK. 2025. hal-04943043

HAL Id: hal-04943043

https://hal.science/hal-04943043v1

Preprint submitted on 12 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# GENERALIZED KURAMOTO-SIVASHINSKY EQUATION SOLVING BY B-SPLINES BASED NEURAL NETWORK

JAMAL ADETOLA<fnref n="1" />, CHARBEL MAMLANKOU<fnref n="2" /> AND BRICE ADJALLA<fnref n="3" />

**Abstract.** This study explores the resolution of the generalized Kuramoto-Sivashinsky equation using Variational Physics-Informed Neural Networks (VPINNs) combined with B-Spline test functions. The Kuramoto-Sivashinsky equation, a nonlinear fourth-order PDE, is widely used to model instabilities in fluid dynamics, combustion, and pattern formation. Our approach leverages the variational formulation to improve solution accuracy while ensuring convergence. Numerical experiments demonstrate that VPINNs effectively capture the complex dynamics of the system, achieving an error order below $10^{-2}$. The robustness of the model is analyzed through different network architectures and optimization strategies, including Adam and L-BFGS. The results highlight VPINNs' potential for solving high-order nonlinear PDEs while maintaining computational efficiency. Future work will focus on optimizing computational cost and exploring adaptive collocation techniques, as well as exploring higher orders?

**Keywords:** Kuramoto-Sivashinsky equation, VPINNs, B-Splines, numerical analysis, nonlinear PDEs, variational methods

**2020 Mathematics Subject Classification.** 65E02, 65E02.

## CONTENTS

<table><tr><td>Introduction</td><td>2</td></tr><tr><td>1. Preliminaries</td><td>2</td></tr><tr><td>2. Background</td><td>8</td></tr><tr><td>3. Numerical experiments</td><td>10</td></tr><tr><td>4. Discussion</td><td>12</td></tr><tr><td>Conclusion</td><td>13</td></tr><tr><td>References</td><td>14</td></tr></table>

<footnote>*Keywords and phrases:* Kuramoto-Sivashinsky equation, VPINNs, B-Splines, numerical analysis, nonlinear PDEs, variational methods</footnote>

<footnote n="1">¹ adetolajamal65@gmail.com</footnote>

<footnote n="2">² charbelzeusmamlankou@gmail.com</footnote>

<footnote n="3">³ bricebonaparteadjalla2@gmail.com</footnote>
# INTRODUCTION

The Kuramoto-Sivashinsky equation (KSE) is a nonlinear fourth-order partial differential equation that was introduced to model instabilities in various physical systems. It was independently derived by Kuramoto in the context of phase turbulence in chemical reactions [16] and by Sivashinsky to describe thermo-diffusive instabilities in laminar flame fronts [21]. Since its inception, the KSE has been widely studied due to its ability to exhibit complex dynamical behaviors, including chaotic solutions, bifurcations, and pattern formation [12, 14].

This equation, which involves fourth-order spatial derivatives coupled with nonlinear terms, describes the evolution of a physical variable such as velocity or interface height. The KSE serves as a canonical model for studying spatiotemporal chaos and nonlinear dynamics in dissipative systems [5]. Its solutions exhibit a rich variety of behaviors depending on the domain size and boundary conditions, including steady states, traveling waves, and fully chaotic regimes [4, 24].

The applications of the KSE are numerous and span multiple fields. In combustion theory, it is used to model instabilities in laminar flame fronts [22]. In fluid dynamics, it describes thin film flows and interfacial instabilities in viscous fluids [18]. In applied mathematics, it has become a paradigmatic example for studying deterministic chaos and attractors in infinite-dimensional dynamical systems [23]. Extensions of the KSE to higher dimensions or with additional terms (e.g., Benney-type equations) have been employed to study multi-scale interactions and stabilization mechanisms in complex systems [3, 20].

The Kuramoto-Sivashinsky equation remains a cornerstone in the study of nonlinear dynamics due to its mathematical simplicity and its ability to capture key features of chaotic systems. Recent research has focused on advanced numerical methods and machine learning approaches, such as physics-informed neural networks (PINNs), to further explore its solutions and their physical implications [13, 19].

# 1. PRELIMINARIES

In this section, we establish the mathematical foundations necessary for the study of the generalized Kuramoto-Sivashinsky. We begin by defining the functional spaces in which the problem is well-posed, followed by a discussion on existence and uniqueness results. Key theorems, lemmas, and fundamental properties related to the equation, including energy estimates, Galilean invariance, and stability considerations, are presented to provide a rigorous theoretical framework for our numerical approach.

Consider the generalized Kuramoto-Sivashinsky equation, written as follows:

$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \alpha \frac{\partial^2 u}{\partial x^2} + \beta \frac{\partial^3 u}{\partial x^3} + \gamma \frac{\partial^4 u}{\partial x^4} = 0, \quad x \in \Omega \subset \mathbb{R}, t > 0 \quad (1.1) $$

where $\alpha, \beta$ and $\gamma$ are arbitrary real constants. to ensure a stabilizing diffusion term, we impose $\alpha, \beta, \gamma > 0$, $(\gamma < \beta)$.

In the context of the Kuramoto-Sivashinsky equation, the variable $u$ typically represents a physical quantity that varies with both space and time. Specifically, $u(x, t)$ can be interpreted a height of an interface, velocity field and temperature or concentration.

For the rest of our work, we introduce functional spaces, in particular the Sobolev space for the $u$ solution, if it exists (we'll prove this below).
**Proposition 1.1.** To study the existence and uniqueness of solutions, it is natural to consider the Sobolev spaces $H^s(\Omega)$, defined as:

$$H^s(\Omega) = \left\{ u \in L^2(\Omega) \mid \frac{\partial^k u}{\partial x^k} \in L^2(\Omega), \forall k \le s \right\}.$$

The analysis is often conducted in a bounded domain $\Omega = [-L, L]$ with periodic boundary conditions $u(-L, t) = u(L, t)$ or Dirichlet-type conditions $u = u_0(x) \in H_{per}^2(\Omega)$, $u_x(-L, t) = u_x(L, t)$ on the boundary. The best choice of space to control the regularity of fourth-order of space derivative is $H_{per}^s$, $s \ge 2$, the Sobolev space with periodic boundary conditions. Let introduce some theorem.

**Theorem 1.2** (Lumer-Phillips [17]). Let A be a linear operator defined on a linear subspace $D(A)$ of the Banach space X. Then A generates a contraction semigroup if and only if:

(1) $D(A)$ is dense in X

(2) A is dissipative, and

(3) $A - \lambda_0 I$ is surjective for some $\lambda_0$, where I denotes the identity operator

An operator satisfying the last two conditions is called maximally dissipative.

**Theorem 1.3** (Uniformly Continuous Semigroup [7]). A strongly continuous semigroup $T(t)$ is said to be **uniformly continuous** if

$$\lim_{t \to 0^+} \|T(t) - I\| = 0.$$

In this case, the infinitesimal generator A of T is bounded, and we have $D(A) = X$, as well as the exponential representation:

$$T(t) = e^{At} := \sum_{k=0}^{\infty} \frac{A^k}{k!} t^k.$$

Conversely, any bounded operator $A: X \to X$ is the infinitesimal generator of a uniformly continuous semigroup given by:

$$T(t) := e^{At}.$$

Thus, a linear operator A is the infinitesimal generator of a uniformly continuous semigroup if and only if A is a bounded linear operator. Moreover, if X is a finite-dimensional Banach space, then any strongly continuous semigroup is uniformly continuous. Finally, for a strongly continuous semigroup that is not uniformly continuous, the infinitesimal generator A is unbounded, and in this case, the exponential series $e^{At}$ may not converge.

**Theorem 1.4** (Wirtinger's Inequality [9]). Let a, b be two real numbers such that $a < b$, and let f be a continuous function of class $C^1$ piecewise on $[a, b]$, with complex values. If the integral of f over $[a, b]$ is zero and if $f(a) = f(b)$, then the following upper bound, known as Wirtinger's inequality, holds:

$$\int_a^b |f(t)|^2 dt \le \left(\frac{b-a}{2\pi}\right)^2 \int_a^b |f'(t)|^2 dt.$$

**Theorem 1.5** (Existence and uniqueness). The generalize Kuramoto-Sivashinsky equation (1.1) has a unique solution in $H^s(\Omega)$
**Proof:** Let $A(u) = -\alpha\partial_x^2 u - \beta\partial_x^3 u - \gamma\partial_x^4 u$ and $N(u) = -u\partial_x u$. The equation can therefore be written as a Cauchy problem in the following form:

$$\partial_t u = A(u) - N(u)$$

where $A$ is the linear operator and $N$, the nonlinear operator. The operator $A(u)$ is defined on $D(A) = H^s(\Omega)$ with periodic edge conditions. We use the theorem (1.2) to show that this operator generates a contracting semigroup. First, let's show that $A$ is dissipative. This is equivalent to showing that for all $u \in D(A)$, we have $(A(u), u) \le 0$, where $(\cdot, \cdot)$ is scalar product. For each term of $A$, we have after after several integrations by parts:

$$(-\alpha \frac{\partial^2 u}{\partial x^2}, u) = \alpha \int_{\Omega} \left(\frac{\partial u}{\partial x}\right)^2 dx \ge 0$$

because the integral function is positive [11] for the positive function $(\frac{\partial u}{\partial x})^2$. As same for the rest of term of $A$, we have:

$$(-\beta \frac{\partial^3 u}{\partial x^3}, u) = -\beta \int_{\Omega} \left(\frac{\partial^2 u}{\partial x^2}\right)^2 dx \le 0, \quad (-\gamma \frac{\partial^4 u}{\partial x^4}, u) = \gamma \int_{\Omega} \left(\frac{\partial^2 u}{\partial x^2}\right)^2 dx \ge 0$$

Remember that all the terms in square brackets in the integral are due to the periodic conditions we impose on our model.

$$(A(u), u) = \alpha \int_{\Omega} \left(\frac{\partial u}{\partial x}\right)^2 dx + (\gamma - \beta) \int_{\Omega} \left(\frac{\partial^2 u}{\partial x^2}\right)^2 dx, \quad (\gamma - \beta) < 0$$

It is therefore necessary to show that the term with the second derivative is dominant in this equality. Using the Wirtinger inequality (1.4), we have:

$$\int_{\Omega} \left(\frac{\partial u}{\partial x}\right)^2 dx \le C \int_{\Omega} \left(\frac{\partial^2 u}{\partial x^2}\right)^2 dx, \quad C = \frac{L^2}{\pi^2} \\ \Rightarrow \alpha \int_{\Omega} \left(\frac{\partial u}{\partial x}\right)^2 dx \le \alpha C \int_{\Omega} \left(\frac{\partial^2 u}{\partial x^2}\right)^2 dx$$

So for any $(\beta - \gamma) \ge \alpha$, we can be sure that the term with the second derivative is dominant and therefore for any $u \in D(A)$, $(A(u), u) \le 0$. This new condition will be applied when choosing model parameters for numerical experiments. A second necessary condition is to use the notion of image to satisfy the second point of the Lumer-Phillips theorem (1.2). The proof of the surjectivity of $(\lambda_0 I - A)$. For an operator $(\lambda_0 I - A)$, $\lambda_0 \in \mathbb{R}$, we define by $R(\lambda_0 I - A)$, the set of all values that the operator $(\lambda_0 I - A)$ can take.

$$R(\lambda_0 I - A) = \{f \in H^s(\Omega) \mid \exists u \in D(A) \text{ such that } (\lambda_0 I - A)u = f\}.$$

In our context, showing that $R(\lambda_0 I - A) = H^s(\Omega)$ means that for all $f \in H^s(\Omega)$, there exists a solution $u \in D(A)$ to the equation $(\lambda_0 I - A)u = f$, that is, the equation has a solution for any right-hand side $f$. This equation can be rewritten as:

$$\lambda_0 u + \alpha \frac{\partial^2 u}{\partial x^2} + \beta \frac{\partial^3 u}{\partial x^3} + \gamma \frac{\partial^4 u}{\partial x^4} = f.$$
Taking the Fourier transform [1] on both sides and denoting $\hat{u}$ as the Fourier transform of $u$, we obtain:

$$ \lambda_0 \hat{u}(\xi) + \alpha \xi^2 \hat{u}(\xi) + i\beta \xi^3 \hat{u}(\xi) + \gamma \xi^4 \hat{u}(\xi) = \hat{f}(\xi). $$

Solving for $\hat{u}(\xi)$, we get $\hat{u}(\xi) = \frac{\hat{f}(\xi)}{\lambda_0 + \alpha\xi^2 + i\beta\xi^3 + \gamma\xi^4}$. The $H^s$ norm can be expressed as:

$$ \|u\|_{H^s}^2 = \|u\|_{L^2}^2 + \left\| \frac{\partial u}{\partial x} \right\|_{L^2}^2 + \left\| \frac{\partial^2 u}{\partial x^2} \right\|_{L^2}^2 + \cdots + \left\| \frac{\partial^s u}{\partial x^s} \right\|_{L^2}^2 $$

In the Fourier domain, this becomes:

$$ \|u\|_{H^s}^2 = \int |\hat{u}(\xi)|^2 d\xi + \int |\xi|^2 |\hat{u}(\xi)|^2 d\xi + \int |\xi|^4 |\hat{u}(\xi)|^2 d\xi + \cdots + \int |\xi|^{2s} |\hat{u}(\xi)|^2 d\xi $$

This sum can be rewritten in a more compact form as $\|u\|_{H^s}^2 = \int (1 + |\xi|^2)^s |\hat{u}(\xi)|^2 d\xi$. To ensure that this solution belongs to $H^s(\Omega)$, we must verify that: $\int (1 + |\xi|^2)^s |\hat{u}(\xi)|^2 d\xi < \infty$.

Substituting the expression for $\hat{u}(\xi)$, this reduces to proving the convergence of:

$$ \int \frac{(1 + |\xi|^2)^s |\hat{f}(\xi)|^2}{|\lambda_0 + \alpha\xi^2 + i\beta\xi^3 + \gamma\xi^4|^2} d\xi. $$

To analyze the denominator $|\lambda_0 + \alpha\xi^2 + i\beta\xi^3 + \gamma\xi^4|^2$, we observe that for sufficiently large $|\xi|$, the dominant term is $\gamma\xi^4$. Thus, for sufficiently large $\lambda_0$, we obtain a lower bound:

$$ |\lambda_0 + \alpha\xi^2 + i\beta\xi^3 + \gamma\xi^4|^2 \geq C(\lambda_0 + |\xi|^4), C > 0. $$

Therefore, we can estimate $\int (1 + |\xi|^2)^s |\hat{u}(\xi)|^2 d\xi \leq C \int \frac{(1 + |\xi|^2)^s |\hat{f}(\xi)|^2}{(\lambda_0 + |\xi|^4)^2} d\xi$.

For sufficiently large $\lambda_0$, this integral remains finite for $f \in H^s(\Omega)$, because the numerator grows as $|\xi|^{2s}$ and the denominator grows as $|\xi|^8$. We conclude that $u \in H^s(\Omega)$ and establish the bound:

$$ |u|_{H^s} \leq C|f|_{H^s}. $$

Thus, we conclude that $R(\lambda I - A) = H^s(\Omega)$ for sufficiently large $\lambda$. This estimates proves that the solution $u$ exists, it is unique, and depends continuously on $f$.

**Theorem 1.6 (Galilean Invariance [8]).** *The generalized Kuramoto-Sivashinsky equation (1.1) is invariant under the Galilean transformation:*

$$ x' = x - ct, \quad t' = t, \quad u' = u + c, $$

where $c$ is a constant representing the velocity of the reference frame.

**Proof:** We first compute the transformed derivatives under the Galilean transformation. Using the chain rule, we obtain:

$$ \frac{\partial u'}{\partial t'} = \frac{\partial u}{\partial t} - c \frac{\partial u}{\partial x}, \quad \frac{\partial u'}{\partial x'} = \frac{\partial u}{\partial x}, \quad \frac{\partial^2 u'}{\partial x'^2} = \frac{\partial^2 u}{\partial x^2}, \quad \frac{\partial^3 u'}{\partial x'^3} = \frac{\partial^3 u}{\partial x^3}, \quad \frac{\partial^4 u'}{\partial x'^4} = \frac{\partial^4 u}{\partial x^4} $$

Substituting these expressions into the transformed equation:

$$ \frac{\partial u'}{\partial t'} + u' \frac{\partial u'}{\partial x'} + \alpha \frac{\partial^2 u'}{\partial x'^2} + \beta \frac{\partial^3 u'}{\partial x'^3} + \gamma \frac{\partial^4 u'}{\partial x'^4} = 0, $$
we obtain:

$$ \left( \frac{\partial u}{\partial t} - c \frac{\partial u}{\partial x} \right) + (u+c) \frac{\partial u}{\partial x} + \alpha \frac{\partial^2 u}{\partial x^2} + \beta \frac{\partial^3 u}{\partial x^3} + \gamma \frac{\partial^4 u}{\partial x^4} = 0. $$

Expanding the terms:

$$ \frac{\partial u}{\partial t} - c \frac{\partial u}{\partial x} + u \frac{\partial u}{\partial x} + c \frac{\partial u}{\partial x} + \alpha \frac{\partial^2 u}{\partial x^2} + \beta \frac{\partial^3 u}{\partial x^3} + \gamma \frac{\partial^4 u}{\partial x^4} = 0. $$

The terms $-c\frac{\partial u}{\partial x}$ and $c\frac{\partial u}{\partial x}$ cancel out, leading back to the original equation:

$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \alpha \frac{\partial^2 u}{\partial x^2} + \beta \frac{\partial^3 u}{\partial x^3} + \gamma \frac{\partial^4 u}{\partial x^4} = 0. $$

Since the transformed equation retains the same form as the original, the generalized Kuramoto-Sivashinsky equation is invariant under Galilean transformations. This invariance is significant because it ensures that the equation maintains its structure in any Galilean reference frame, which aligns with the principle of Galilean relativity in physics.

**Proposition 1.7 (Reflection Symmetry [10]).** If $u(x, t)$ is a solution of the generalized Kuramoto-Sivashinsky equation (1.1), then the function $u^*(x, t) = -u(-x, t)$ is also a solution.

**Proof:** We define the transformed function $u^*(x, t) = -u(-x, t)$. We substitute $u^*(x, t)$ into the equation and verify that it satisfies the same equation. Using the chain rule, we compute the derivatives:

$$ \frac{\partial u^*}{\partial t} = -\frac{\partial u}{\partial t}(-x, t), \frac{\partial u^*}{\partial x} = -\frac{d}{dx}u(-x, t) = -\left(-\frac{\partial u}{\partial x}(-x, t)\right) = \frac{\partial u}{\partial x}(-x, t), \frac{\partial^2 u^*}{\partial x^2} = \frac{\partial^2 u}{\partial x^2}(-x, t) $$

$$ \frac{\partial^3 u^*}{\partial x^3} = -\frac{\partial^3 u}{\partial x^3}(-x, t), \frac{\partial^4 u^*}{\partial x^4} = \frac{\partial^4 u}{\partial x^4}(-x, t) $$

We now substitute these expressions into (1.1):

$$ \frac{\partial u^*}{\partial t} + u^* \frac{\partial u^*}{\partial x} + \alpha \frac{\partial^2 u^*}{\partial x^2} + \beta \frac{\partial^3 u^*}{\partial x^3} + \gamma \frac{\partial^4 u^*}{\partial x^4} = 0. $$

Replacing the transformed derivatives:

$$ -\frac{\partial u}{\partial t}(-x, t) + (-u(-x, t)) \frac{\partial u}{\partial x}(-x, t) + \alpha \frac{\partial^2 u}{\partial x^2}(-x, t) - \beta \frac{\partial^3 u}{\partial x^3}(-x, t) + \gamma \frac{\partial^4 u}{\partial x^4}(-x, t) = 0. $$

Expanding $-\left[\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \alpha \frac{\partial^2 u}{\partial x^2} + \beta \frac{\partial^3 u}{\partial x^3} + \gamma \frac{\partial^4 u}{\partial x^4}\right](-x, t) = 0$. Since $u(x, t)$ is a solution, we know that:

$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \alpha \frac{\partial^2 u}{\partial x^2} + \beta \frac{\partial^3 u}{\partial x^3} + \gamma \frac{\partial^4 u}{\partial x^4} = 0. $$

Thus, $-[0] = 0$. Since $u^*(x, t)$ satisfies the same equation, it is also a solution.
**Proposition 1.8** (Energy Functional [6]). Let $u(x, t)$ be a sufficiently smooth solution of the generalized Kuramoto-Sivashinsky equation. The associated energy functional is defined as:

$$E(t) = \frac{1}{2} \int_{\Omega} u^2(x, t) dx.$$

**Remark 1.9.** The energy functional time evolution satisfies the dissipation law:

$$\frac{dE}{dt} = -\alpha \int_{\Omega} (u_x)^2 dx - \gamma \int_{\Omega} (u_{xx})^2 dx.$$

Therefore, the system exhibits an energy dissipation mechanism governed by the parameters $\alpha$ and $\gamma$.

In fact, we multiply (1.1) by $u(x, t)$ and integrate over $\Omega$:

$$\int_{\Omega} u \frac{\partial u}{\partial t} dx + \int_{\Omega} uu_x dx + \alpha \int_{\Omega} uu_{xx} dx + \beta \int_{\Omega} uu_{xxx} dx + \gamma \int_{\Omega} uu_{xxxx} dx = 0.$$

By analyzing each term, we have:

$$\int_{\Omega} u \frac{\partial u}{\partial t} dx = \frac{1}{2} \frac{d}{dt} \int_{\Omega} u^2 dx = \frac{dE}{dt} \quad \text{and} \quad \int_{\Omega} uu_x dx = \frac{1}{2} \int_{\Omega} \frac{d}{dx} (u^2) dx.$$

If the boundary conditions ensure that this term vanishes (periodic conditions), then $\int_{\Omega} uu_x dx = 0$. By integrating by parts the diffusion term, we get $\int_{\Omega} uu_{xx} dx = -\int_{\Omega} (u_x)^2 dx$. This term is always negative and favors energy dissipation. After integration by parts of $\int_{\Omega} uu_{xxx} dx$ and assuming adequate boundary conditions, the third derivative term give us $\int_{\Omega} uu_{xxx} dx = 0$. It does not directly contribute to energy evolution. After two integrations by parts for the stabilization term (fourth derivative) and using boundary conditions, we have $\int_{\Omega} uu_{xxxx} dx = -\int_{\Omega} (u_{xx})^2 dx$. This term is also negative and stabilizes the system.

By combining these results, we obtain:

$$\frac{dE}{dt} = -\alpha \int_{\Omega} (u_x)^2 dx - \gamma \int_{\Omega} (u_{xx})^2 dx.$$

This result shows that energy decreases over time, indicating a dissipation effect. More precisely, $\alpha > 0$ ensures energy dissipation through viscous diffusion and $\gamma > 0$ stabilizes high frequencies by damping rapid oscillations. We also note that the term $\beta u_{xxx}$ does not contribute directly to dissipation, as it is purely dispersive. In view of all our above assumptions, we consider the following model [25] for the rest of the study over the domain $\Omega = [-10, 10]$:

$$\left\{
\begin{array}{ll}
\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \alpha \frac{\partial^2 u}{\partial x^2} + \beta \frac{\partial^3 u}{\partial x^3} + \gamma \frac{\partial^4 u}{\partial x^4} = 0, & x \in \Omega, \quad t \in [0, T], \\
u(x, 0) = 15 \tanh\left(\frac{x}{2}\right) - 15 \tanh\left(\frac{x}{2}\right) - 15 \tanh\left(\frac{x}{2}\right) + 11, & x \in \Omega, \\
u(10, t) = u(-10, t) = u_L(t), & t > 0, \\
\frac{\partial u}{\partial x}(10, t) = \frac{\partial u}{\partial x}(-10, t) = v_L(x, t), & t > 0,
\end{array}
\right.
\quad (1.2)$$
where

$$u_L(x) = -15 \tanh(t-5)^3 - 15 \tanh(t-5)^2 + 15 \tanh(t-5) + 11$$

and

$$v_L(x,t) = -15(\tanh^2(t-5) - 1)\tanh(t-5) - 15\left(\frac{3}{2}\tanh^2(t-5) - \frac{3}{2}\right)\tanh^2(t-5) + \\ + \frac{15}{2}\tanh^2(t-5) - \frac{15}{2}$$

The rest of our study is subdivided into 3 main sections, in which we detail the neural network method we are applying to our system (section 2), then move on to a numerical experiment (section 3) to find out the predictive power of our model (1.2), and finally a discussion (section 4), which will serve as a point of appreciation of the work and the prospects that lie ahead.

## 2. BACKGROUND

Several neural network-based techniques exist for solving partial differential equations, but the variational [15] approach has proved robust for nonlinear problems. Called Variational Physics-Informed Neural Networks (VPINNs), we combine a finite element approach with determining the weak formulation and minimizing this formulation as residuals.

Let's define the approximation form $\tilde{u}(x, t) = u_{NN}(x, t; w, b)$ where w and b are the weights and biases of the neural network respectively.

To determine the loss function, we write the weak formulation and multiply by a test function $v(x, t)$ before proceeding to an integration over the $\Omega$ spatial domain:

$$\int_{\Omega} \left( \frac{\partial \tilde{u}}{\partial t} + \tilde{u} \frac{\partial \tilde{u}}{\partial x} + \alpha \frac{\partial^2 \tilde{u}}{\partial x^2} + \beta \frac{\partial^3 \tilde{u}}{\partial x^3} + \gamma \frac{\partial^4 \tilde{u}}{\partial x^4} \right) v \, dx = 0.$$

Applying integration by parts where necessary and incorporating boundary conditions, we obtain the final weak formulation:

$$\int_{\Omega} \frac{\partial \tilde{u}}{\partial t} v \, dx + \int_{\Omega} \tilde{u} \frac{\partial \tilde{u}}{\partial x} v \, dx - \alpha \int_{\Omega} \frac{\partial \tilde{u}}{\partial x} \frac{\partial v}{\partial x} \, dx + \beta \int_{\Omega} \frac{\partial^2 \tilde{u}}{\partial x^2} \frac{\partial v}{\partial x} \, dx + \gamma \int_{\Omega} \frac{\partial^2 \tilde{u}}{\partial x^2} \frac{\partial^2 v}{\partial x^2} \, dx = 0.$$

This weak form is used in the VPINNs framework, where the loss function is constructed as follows:

$$\mathcal{L} = \mathcal{L}_{res} + \lambda_{IC}\mathcal{L}_{IC} + \lambda_{BC}\mathcal{L}_{BC}, \qquad (2.1)$$

where:

$$\mathcal{L}_{res} = \sum_{i=1}^{N_r} \left| \int_{\Omega} (\tilde{u}_t v + \tilde{u} \tilde{u}_x v - \alpha \tilde{u}_x v_x + \beta \tilde{u}_{xx} v_x + \gamma \tilde{u}_{xx} v_{xx}) dx \right|^2, \quad (2.2)$$

$$\mathcal{L}_{IC} = \sum_{i=1}^{N_{IC}} |\tilde{u}(x_i, 0) - g_{IC}(x_i)|^2, \qquad (2.3)$$

$$\mathcal{L}_{BC} = \sum_{i=1}^{N_{BC}} (|\tilde{u}(x_i, t) - g_{BC}(x_i, t)|^2 + |\tilde{u}_x(x_i, t) - h_{BC}(x_i, t)|^2), \quad (2.4)$$
where $N_r$, $N_{IC}$, and $N_{BC}$ are the number of collocation points for the residual, initial, and boundary conditions respectively. The functions $g_{IC}$ and $g_{BC}$ are the prescribed initial and boundary conditions, and $h_{BC}$ represents the boundary derivative constraints.

We introduce then some theoretical foundations of **B-Splines**: B-splines of degree $p$ are recursively defined on a knot sequence $T = \{t_i\}$ as follows:

$$ \text{For } p = 0, N_{i,0}(x) = \begin{cases} 1 & \text{if } t_i \le x < t_{i+1} \\ 0 & \text{otherwise} \end{cases} $$

$$ \text{For } p > 0, N_{i,p}(x) = \omega_1(x)N_{i,p-1}(x) + \omega_2(x)N_{i+1,p-1}(x) $$

$$ \text{where } \omega_1(x) = \frac{x-t_i}{t_{i+p}-t_i} \quad \omega_2(x) = \frac{t_{i+p+1}-x}{t_{i+p+1}-t_{i+1}} $$

**Proposition 2.1** (Compact Support of B-Splines [2]). Let $N_{i,p}$ be a B-spline basis function of degree $p$ defined on a knot sequence $T = \{t_j\}_{j=1}^{n+p+1}$. Then:

(1) The support of $N_{i,p}$ is compact and given by $\text{supp}(N_{i,p}) = [t_i, t_{i+p+1}]$

(2) For any $x \notin [t_i, t_{i+p+1}]$, $N_{i,p}(x) = 0$

(3) If $x \in (t_j, t_{j+1})$ where $i \le j < i + p$, then $N_{i,p}(x) > 0$

**Proof:** For $p = 0$, by definition:

$$ N_{i,0}(x) = \begin{cases} 1 & \text{if } x \in [t_i, t_{i+1}) \\ 0 & \text{otherwise} \end{cases} $$

Thus, $\text{supp}(N_{i,0}) = [t_i, t_{i+1}]$, verifying the property for $p = 0$.

For degree $p$, consider the recursive formula:

$$ N_{i,p}(x) = \frac{x-t_i}{t_{i+p}-t_i} N_{i,p-1}(x) + \frac{t_{i+p+1}-x}{t_{i+p+1}-t_{i+1}} N_{i+1,p-1}(x) $$

By the induction hypothesis:

$$ \text{supp}(N_{i,p-1}) = [t_i, t_{i+p}] $$

$$ \text{supp}(N_{i+1,p-1}) = [t_{i+1}, t_{i+p+1}] $$

Therefore, for $x < t_i$: both terms are zero since $x$ is outside both supports; for $x > t_{i+p+1}$: both terms are zero for the same reason; for $x \in [t_i, t_{i+p+1}]$: at least one term is non-zero when the denominator is non-zero.

For the positivity claim in (3), the coefficients in the recursive formula are positive in $(t_j, t_{j+1})$. By induction, at least one of $N_{i,p-1}$ or $N_{i+1,p-1}$ is positive, which implies that $N_{i,p}(x) > 0$ in this interval

**Corollary 2.2** (Local Support Property). At any point $x \in [t_p, t_{n-p}]$, at most $p+1$ B-spline basis functions of degree $p$ are non-zero.

**Proposition 2.3** (Partition of Unity for B-Splines). Let $\{N_{i,p}\}_{i=1}^n$ be the set of B-spline basis functions of degree $p$ defined on a knot sequence $T = \{t_i\}_{i=1}^{n+p+1}$ where knots satisfy $t_1 \le t_2 \le \dots \le t_{n+p+1}$. Then, the B-spline basis functions form a partition of unity on the interval $[t_p, t_{n-p}]$, i.e.,

$$ \sum_{i=1}^{n} N_{i,p}(x) = 1 \quad \forall x \in [t_p, t_{n-p}] $$
Moreover, this property holds independently of the knot multiplicities, provided they do not exceed $p + 1$.

**Proof:** The proof proceeds by induction on the degree $p$.

For $p = 0$, the property is trivial since the B-splines are characteristic functions of intervals $[t_i, t_{i+1})$ that partition $[t_1, t_n]$.

Assume the property holds for degree $p - 1$. For degree $p$, using the recursive definition:

$$ \begin{aligned} \sum_{i=1}^{n} N_{i,p}(x) &= \sum_{i=1}^{n} [\omega_{1}^{i}(x) N_{i,p-1}(x) + \omega_{2}^{i}(x) N_{i+1,p-1}(x)] \\ &= \sum_{i=1}^{n} \frac{x - t_{i}}{t_{i+p} - t_{i}} N_{i,p-1}(x) + \sum_{i=1}^{n} \frac{t_{i+p+1} - x}{t_{i+p+1} - t_{i+1}} N_{i+1,p-1}(x) \end{aligned} $$

After collecting terms and using the induction hypothesis:

$$ \sum_{i=1}^{n} N_{i,p}(x) = \sum_{i=1}^{n} N_{i,p-1}(x) = 1 $$

The domain restriction $[t_p, t_{n-p}]$ ensures that all basis functions contributing to the sum are well-defined.

We now introduce the discrete space that will serve as the basis for choosing our test function. The test space $V_h$ is defined as:

$$ V_h = \operatorname{span}\{N_{i,p} : i = 1, \dots, n\} \qquad (2.5) $$

### 3. NUMERICAL EXPERIMENTS

In this section we move on to the numerical simulation of the model with our VPINNs network based on B-Splines. The analytical solution [25] of the problem is known for the case $\alpha = \gamma = 1$, $\beta = 4$ and will serve as a basis for comparison and calculation of model errors. It's defined by:

$$ u_{\text{exact}} = 11 + 15\tanh(\theta) - 15\tanh^2(\theta) - 15\tanh^3(\theta) $$

where $\theta = t - \frac{x}{2}$. The conditions imposed in the model remain unchanged in this section. We're using the PyTorch library with Python to set up our network, which comprises 2 inputs, 3 hidden layers with $n$ (variables) neurons and one output, all on a MacBook pro 2017, core i7, RAM 16Gb with a frequency of 3.1Ghz.

We use 20 test functions and a degree 3 B-Spline for our experiment. The number of quadratures is variable for calculating the integrals involved in the loss function. The simulation time is set to 2s.

For a first attempt, we have set a maximum number of iterations of 1000, as the aim is to obtain a very satisfactory result in a reduced number of iterations. A quadrature of 128 is chosen under a network with $n=16$ neurons per hidden layer. The optimizer used is Adam with a learning step of $10^{-3}$ and the runtime is 2 minutes 18 seconds. The results are as follows:
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Plots of exact (the first two) and predicted (the last two) solutions using **n=16** neurons per hidden layer and **Adam optimizer**, followed by time evolution.</figcaption></figure>

Initial analysis of the graphs reveals that our model is not optimally trained. To improve performance, we decided to proceed in stages, starting by increasing the number of neurons in the hidden layers. This one-off modification will enable us to assess whether the network's low capacity was the main cause of the unsatisfactory results. For this new experiment, we increased the number of neurons to 20 per hidden layer, keeping all other parameters unchanged. We also increase the quadrature to 200 for this part. The runtime here is 5 minutes 26 seconds.

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Plots of predicted solutions using **n=20** neurons per hidden layer and **Adam optimizer**.</figcaption></figure>

Comparison with the exact solution shows that the prediction, although improved in its temporal and spatiotemporal evolution, still shows significant deviations, particularly in terms of maximum and minimum amplitude.

This observation suggests that increasing the number of neurons and quadrature alone is not enough to achieve the desired accuracy. To further refine our model, we apply an additional optimization phase using the L-BFGS algorithm, retaining the network weights obtained with n=20 neurons as a starting point. This optimization is configured to run 20 learning epochs with 500 iterations per epoch and a learning rate of 10⁻³. The total runtime at this step is 9 minutes 55 seconds. This refinement approach aims to further minimize the loss function while preserving the initial learning gains.

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Plots of predicted solutions using **n=20** neurons per hidden layer with **Adam + L-BFGS optimizer**.</figcaption></figure>
Detailed analysis of the maximum and minimum values of the predictions shows excellent agreement with the exact solution, confirming the convergence of the model and demonstrating the quality of its training.

We now present the losses of our neural models to see for ourselves how effective they really are.

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Loss function: 16 neurons with Adam (left); 20 neurons with Adam and L-BFGS (right)</figcaption></figure>

With our first algorithm, we obtain a minimum loss of the order of $10^{-3}$ and with the second and third attempts, we obtain a minimum loss of the order of $10^{-3}$ approximately.

We present the temporal evolution of the absolute, relative, $L^2$ and $L_\infty$ errors between the predicted solution and the exact solution.

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Errors plots</figcaption></figure>

## 4. DISCUSSION

At the end of our study, it became clear that VPINNs are particularly robust for handling nonlinear and high-order problems. The various approaches implemented in our model—such as the use of B-Spline test functions and the variational formulation—have allowed us to achieve remarkable convergence and near conformity with the expected exact solution.

This performance attests to the ability of VPINNs to accurately capture the complex dynamics of the system modeled by (1.2), even under conditions where classical methods might struggle. In particular, obtaining an error order below $10^{-2}$ demonstrates that, despite the intrinsic complexity of the problem, our model is able to faithfully reproduce the overall trend and behavior of the phenomenon under study.

From an optimization standpoint, we are considering several avenues for improvement. One major direction concerns reducing the model's complexity with the goal of lowering computational cost while maintaining the approach's accuracy. This reduction could be achieved by simplifying the neural network or by optimizing the number and placement of collocation points, without compromising the model's ability to represent the system dynamics.
Moreover, the results emphasize the importance of judiciously balancing the dif-
ferent terms in the loss function. By carefully adjusting the coefficients associated
with the PDE residual, boundary conditions, and regularization, it is possible to
steer the optimization toward a solution that not only converges rapidly but also
adheres to the expected physical properties.

Finally, this work paves the way for new perspectives, notably regarding the
integration of adaptive approaches for the placement of collocation points and the
use of advanced optimization techniques, such as gradually transitioning to L-BFGS-
type optimizers. These prospects are particularly promising for tackling even more
complex problems and further reducing computational costs while ensuring high
fidelity in reproducing the physical phenomena.

In summary, this study highlights the relevance and power of VPINNs in ad-
dressing highly complex problems, while also underlining the potential optimization
avenues that could make these methods even more performant and efficient in a
wide range of practical applications.

CONCLUSION

Our study demonstrates that VPINNs, combined with B-Spline test functions, provide a robust and accurate method for solving the generalized Kuramoto-Sivashinsky equation. The numerical experiments confirm the model's ability to reproduce the complex dynamical behavior of the equation, achieving high accuracy with an error order below $10^{-2}$. The combination of neural network optimization techniques, including Adam and L-BFGS, ensures both convergence and stability of the solution. Additionally, the use of a variational approach effectively captures the equation's underlying physics while reducing computational cost.

Despite these promising results, some areas remain open for further improvement.
Future work will focus on optimizing the computational complexity by refining the
selection of test functions and improving the efficiency of the optimization process.
Moreover, adaptive collocation strategies could enhance the accuracy of the model
while reducing training time. These advancements will further solidify VPINNs
as a powerful tool for solving complex nonlinear PDEs in various scientific and
engineering applications.
REFERENCES

[1] Stéphane Balac. La transformée de fourier vue sous l'angle du calcul numérique. *HAL Open Science*, 2018.

[2] R. Bartels, J. C. Beatty, and B. A. Barsky. *An Introduction to Splines for Use in Computer Graphics and Geometric Modeling*. Morgan Kaufmann Publishers, San Francisco, 1987.

[3] D. J. Benney. Long waves on liquid films. *Journal of Mathematics and Physics*, 45(1):150–155, 1966.

[4] H. Chang, J. M. Hyman, and B. Nicolaenko. Chaotic dynamics of the kuramoto-sivashinsky equation. *Physical Review Letters*, 68(22):3217–3220, 1992.

[5] M. C. Cross and P. C. Hohenberg. Pattern formation outside of equilibrium. *Reviews of Modern Physics*, 65(3):851–1112, 1993.

[6] J.H.P. Dawes. The kuramoto-sivashinsky equation. *Lecture Notes*, 2006.

[7] K.-J. Engel and R. Nagel. *One-Parameter Semigroups for Linear Evolution Equations*, volume 194 of *Graduate Texts in Mathematics*. Springer, Berlin, 2000.

[8] G. W. Gibbons and S. W. Hawking. Galilean invariance and the principle of relativity. *Physical Review D*, 76(4):044020, 2007.

[9] G. H. Hardy, J. E. Littlewood, and G. Pólya. Inequalities. 1988.

[10] G. M. Homsy and A. A. Nepomnyashchii. Stability of liquid films on inclined planes. *Journal of Fluid Mechanics*, 61(2):221–240, 1974.

[11] Dominique Hulin. Intégration. 2021.

[12] James M. Hyman, Basil Nicolaenko, and Stéphane Zaleski. The kuramoto-sivashinsky equation: A bridge between pdes and dynamical systems. *Physica D: Nonlinear Phenomena*, 18(1-3):113–126, 1986.

[13] George E. Karniadakis, Ioannis G. Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. *Physics-informed machine learning: theory and applications to computational science*. Springer Nature, 2021.

[14] Ioannis G. Kevrekidis, Basil Nicolaenko, and James C. Scovel. A computer assisted study of the kuramoto-sivashinsky equation. *SIAM Journal on Applied Mathematics*, 50(3):760–790, 1990.

[15] Ehsan Kharazmi, Z. Zhang, and George Em Karniadakis. Variational physics-informed neural networks for solving partial differential equations. *ArXiv*, abs/1912.00873, 2019.

[16] Yoshiki Kuramoto. Diffusion-induced chaos in reaction systems. *Progress of Theoretical Physics Supplement*, 64:346–367, 1978.

[17] G. Lumer and R. S. Phillips. A theorem on the generation of strongly continuous semigroups. *Proceedings of the American Mathematical Society*, 12(2):243–246, 1961.

[18] A. Oron, S. H. Davis, and S. G. Bankoff. Long-wave instabilities of thin liquid films. *Reviews of Modern Physics*, 69(3):931–980, 1997.

[19] Maziar Raissi, Paris Perdikaris, and George E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. *Journal of Computational Physics*, 378:686–707, 2019.

[20] Y. Shen and L. Wang. Stabilization of the kuramoto-sivashinsky equation with boundary feedback control. *SIAM Journal on Control and Optimization*, 34(2):676–693, 1996.

[21] Gregory I. Sivashinsky. Nonlinear analysis of hydrodynamic instability in laminar flames—i. derivation of basic equations. *Acta Astronautica*, 4(11-12):1177–1206, 1977.

[22] Gregory I. Sivashinsky. Instabilities, pattern formation, and turbulence in flames. *Annual Review of Fluid Mechanics*, 15(1):179–199, 1983.

[23] Roger Temam. *Infinite-dimensional dynamical systems in mechanics and physics*. Springer-Verlag, 1988.

[24] Leslie S. Tuckerman and Dwight Barkley. Steady states of the kuramoto-sivashinsky equation. *Physica D: Nonlinear Phenomena*, 33(1-3):109–120, 1990.

[25] Kirill Zubov, Zoe McCarthy, Yingbo Ma, Francesco Calisto, Valerio Pagliarino, Simone Azeglio, Luca Bottero, Emmanuel Luján, Valentin Sulzer, Ashutosh Bharambe, et al. Neuralpde: Automating physics-informed neural networks (pinns) with error approximations. *arXiv preprint arXiv:2107.09443*, 2021.
Authors institutions

Université Nationale des Sciences, Technologies, Ingénierie et Mathématiques (UN-
STIM - Bénin)

École Nationale Supérieure de Génie Mathématique et Modélisation (ENSGMM)

*Corresponding author: charbelzeusmamlankou@gmail.com*