A logarithmically improved regularity criterion for the 3D MHD equations in Morrey-Campanato space

1.   Introduction

Nonlinear partial differential equations (PDEs) are crucial for studying nonlinear physical phenomena in various models across different branches. Analyzing qualitative parameters is necessary to comprehend the underlying facts of a model. Nonlinear evolution equations have significant implications in applied sciences, such as optics, crystal lattice theory, fluid dynamics, and plasma physics, making them an excellent source of motivation in recent decades. There are several analytical methods for solving nonlinear partial differential equations (NLPDEs), including the Jacobi elliptic expansion ^[1,2,3], the tanh-function method ^[4,5,6], the improved Q-expansion technique ^[7], the generalized direct algebraic method ^[8], the extended generalized $\left(G^{\prime}/G\right)$-expansion method ^[9,10], and others ^{[11,12,13,14,15,16,17,18,19,20,21]}. In 1995, Gilson and Pickering presented a third-order nonlinear differential equation explaining wave propagation in plasma physics^[22]. This equation is widely recognized:

There are four constants, $\alpha, {\mu}, \gamma, \beta$, in the Gilson-Pickering equation, all of which are nonzero. The equation has three different types: Rosenau-Hyman equation when $\alpha = 0$, ${\mu} = 1$, $\gamma = 0$, $\beta = 3$; Fuchssteiner-Fokas-Camassa-Holm equation when $\alpha = 1$, ${\mu} = -3$, $\beta = 2$; and Fornberg-Whitham equation when $\alpha = 1$, ${\mu} = -1$, $\beta = 3$, $\gamma = 0.5$. Recently, there have been several numerical schemes developed for solving NLPDEs. These include finite element, finite differences ^{[23,24,25,26,27]}, the adaptive moving mesh technique ^{[28,29,30,31,32,33]}, and the parabolic Monge-Ampère (PMA) method ^{[34,35,36,37]}. Lately, various computational methods have been employed to explore the Gilson-Pickering equation, including Chen et al. ^[38] who performed a qualitative analysis of this equation. Clarkson et al. investigated the symmetry reductions of Eq (1.1) using both the classical Lie method and the non-classical method ^[39], the ($1/G^\prime$) method ^[40], and the ($G^\prime/G^2$) method ^[41]. Additionally, Baskonus provides several complex soliton solutions using the Bernoulli sub-equation function method ^[42]. There are many researchers who investigated the Gilson-Pickering equation by different methods ^{[43,44,45,46,47]}. The PMA approach mentioned in this context has yet to be used for Eq (1.1). This motivated us significantly to apply it. It is widely acknowledged that the wave solutions in this situation have areas with rapid spatial variations, such as steep front structures. However, it is important to note that this method's effectiveness in minimizing errors is specifically limited to regions with rapid spatial variations, such as steep front structures. However, this method is effective in minimizing errors in these specific regions.

The following organization is presented in this article: Section 2 provides exact solutions to the Gilson-Pickering equation using the F-expansion method and the ($1/G^\prime$)-expansion method. Section 3 introduces a numerical solution for the same equation on the uniform mesh utilizing implicit finite difference with a study of the stability, error analysis, and convergence of the numerical scheme. Section 4 on the moving mesh method applies the PMA technique to the Gilson-Pickering equation, while the results and discussion are in Section 5. Ultimately, the conclusions are in Section 6.

2.   Traveling wave solution

We consider a nonlinear evolution equation with some physical fields $u(x, t)$ in two variables $x$ and $t$ as follows:

Step 1. We look for the traveling-wave solutions of Eq (2.1) that are formed as follows:

where $k, c$ are constant, and $c$ is the wave speed.

Step 2. Using Eqs (2.1) and (2.2) we directly reduce to

To convert Eq (1.1) into an ordinary differential equation (ODE), substitute the wave transformation Eq (2.2) into it.

where $\prime$ denotes the derivative with respect to $\xi$. We can obtain the following result by integrating once for $\xi$ and putting the constant of integration to zero:

2.1.   The F-expansion method

The F-expansion procedure ^[7] introduces the traveling wave solution of Eq (2.3) of the form

where $a_0, a_j$, and $b_j$ are constants to be determined. The function $F(\xi)$ represents a solution of the following differential equation:

where $A_0$, $A_1$, and $A_2$ are constant ^[7]. By balancing the second derivative of $U$ with the nonlinear term $U^2$ in Eq (2.5) we required that $N+2 = 2N$. Then, we get the value of $N$ as $N = 2$, and Eq (2.6) becomes

To solve for the values of $a_0$, $a_1, a_2$, $b_1, b_2$, $k$, and $c$, we can substitute Eq (2.8) into Eq (2.4) and use Eq (2.7). To create a system of algebraic equations, we collect coefficients of $F(\xi)^{l}$, $l = (-6, -5, \cdots, 5, 6)$. Then, solving this system using Mathematica 13.2 software will give us the desired values as follows:

when $A_0 = 0$, $A_1 = 1$, $A_2 = -1.$

The solution to Eq (1.1) can be determined through Case I and is as follows:

Based on Case II, the solution to Eq (1.1) can be determined:

Using Case III, we can determine the solution for Eq (1.1):

According to Case IV, we can find the solution to Eq (1.1):

when $A_0 = 0, A_1 = -1, A_2 = 1.$

The solution to Eq (1.1) can be determined through Case I and is as follows:

Based on Case II, the solution to Eq (1.1) can be determined:

Using Case III, we can determine the solution for Eq (1.1):

According to Case IV, we can find the solution to Eq (1.1):

2.2.   The (1/G$^\prime$)-expansion method

This subsection presents the important steps of the ($1/G^\prime$)-expansion method. Assume that Eq (2.3) has the solution of the form

The scalars $a_0, a_1, \cdots, a_N$ need to be determined, and $N$ represents the balance terms. Additionally, $G = G(\xi)$ must satisfy the given second-order linear ordinary differential equation

where $\vartheta$ and $\nu$ are constants. The solution of Eq (2.18) is given by

We can use the following format to express the algebraic equation provided in Eq (2.19) as a trigonometric function

After applying the balance principle to Eq (2.5) we required that $N+2 = 2N$, and then we obtain $N = 2$. Equation (2.17) takes the following form:

By inserting Eq (2.21) and its necessary derivatives, along with Eq (2.18) into Eq (2.3), we can obtain the polynomial of $\left(\frac{1}{G^\prime(\xi)}\right)^l$, $l \geq 0$. By summing the coefficients of $\left(\frac{1}{G^\prime(\xi)}\right)^l$ with the same power and setting the resulting summation to zero, we can derive a system of algebraic equations for $a_0, a_1$ and $a_2$. Solving this system using Mathematica 13.2 software yet gives us the values of the variables, which in turn provides the solutions to Eq (1.1).

Family 1:

We can obtain the solution for Eq (1.1) as follows:

In Figures 1–3, we show the exact solutions.

Family 2:

The solution for the traveling wave of Eq (1.1) can be expressed as

where

Family 3:

The solution of Eq (1.1) can be determined as follows:

Figures 4 and 5 show the wave behavior of the solutions $u_9(x, t)$ and $u_{10}(x, t)$, respectively, when changing specific parameters and fixing others. The amplitude increases as $\alpha$ and $\gamma$ increase. Figure 4 shows that when $\alpha$ is positive, the wave is in the positive direction, while in Figure 5, the wave is in the negative direction. However, parameter $\gamma$ has the opposite effect of parameter $\alpha$ on the wave direction between the solutions $u_9(x, t)$ and $u_{10}(x, t)$.

3.   Implicit finite difference method

To obtain the numerical solutions of Eq (1.1) in this section, over the domain $a \leq x \leq b$, $0 \leq t\leq t_f$, where $t_f$ is a certain time, these coordinates are covered by a rectangular grid of points

where $N_x$ represents the nummber of grid points, $\Delta_x = (b-a)/N_x$ represents the distance between any two consecutive points, and $\Delta_t$ denote the temporal increases.

To start, we write the variable $v$ in the following form:

Thus, Eq (1.1) is reformed as

Then, we use the approximation solutions $U^n_m$ to the exact solution $u(x_m, t_n) = u$. To calculate spatial derivatives, we utilize finite difference operators while keeping the temporal differentiation continuous. The weighted average or ($\theta$-method) approximation ^[48] for Eq (3.2) is as follows:

where

We will assume we are using an average with nonnegative weights so that $0 \leq \theta \leq 1$. We will analyze two implicit methods: the fully implicit scheme when $\theta = 1$ and the well-known and popular Crank-Nicolson scheme when $\theta = \frac{1}{2}$. The boundary conditions for Eq (3.2) are $u_{t, a} = u_{t, b} = 0$. We set $t = 0$ in Eq (2.22) to obtain the initial condition. The above system is solved utilizing an ODE solver in FORTRAN 95 called the DASPK solver ^[49,50]. To determine the accuracy of the numerical scheme (3.3), we utilize Taylor expansions to analyze the truncation error. Specifically, we apply Taylor's series expansion for all terms in Eq (3.3) around $U^n_m$, where $u$ represents the exact solution, then we get

Since $u$ is the exact solution of Eq (1.1), for Crank-Nicolson scheme $\theta = \frac{1}{2}$, the first and second brackets are equal to zero. As a result, in both time and space, the scheme is second-order. Therefore, the truncation error of the numerical method becomes $\mathcal{O}(\Delta_x^2, \Delta_t^2)$. With the fully implicit scheme of $\theta = 1$, the first bracket equals zero, and the truncation error becomes $\mathcal{O}(\Delta_x^2, \Delta_t)$. We will examine the stability of the numerical solution. To begin, we will express Eq (1.1) in the following form:

Our numerical scheme's stability is displayed by using von-Neumann stability analysis. The von-Neumann analysis, further called Fourier analysis, is only used when the scheme is linear. The differential Eq (3.5) must be written in the linearized form:

where $\hat{\text{w}} = \max |u|.$ The general form of the implicit scheme for Eq (3.6) can be written as:

Assume that

When using von-Neumann, $|\lambda|\leq 1$ is required for stability. The following result was obtained by substituting Eq (3.8) into Eq (3.7), and canceling the common factor $U^n_m$, and doing some operations

Let

Therefore, Eq (3.9) can be written as

Hence,

As a result, $\lambda$ can be expressed explicitly as

Let us analyze the stability of two schemes: the Crank-Nicolson scheme and the fully implicit scheme. Starting with the Crank-Nicolson scheme, we can determine its stability by plugging in $\theta = \frac{1}{2}$ into Eq (3.12)

Then, using Eq (3.13), we find that $|\lambda| = 1$, meaning that the Crank-Nicolson scheme is linearly unconditionally stable. Moving on to the fully implicit scheme, we substitute the value of $\theta = 1$ into Eq (3.12) to get Eq (3.14)

By analyzing Eq (3.14), we find that $|\lambda| < 1$, indicating that the fully implicit scheme is also linearly unconditionally stable. To discuss the convergence of our numerical method, we have some initial data and perform a sequence of computations using two refined meshes, where $\Delta_x$ and $\Delta_t$ both approach zero. If, for a given point in the rectangular domain $[a, b]\times[0, t_f]$, with coordinates $(x^*, t^*)$, we have $x_m$ and $t_n$ approaching $x^*$ and $t^*$, respectively, this means that the numerical scheme is convergent. In other words, $U^n_m = u(x^*, t^*)$. Assume that

at the point $(x_m, t_{n+1})$, the approximation solution is represented by $U^{n+1}_m$, and the analytical solution is represented by $u(x_m, t_{n+1})$. The error at this point is denoted by $e^{n+1}_m$

We insert Eq (3.16) into Eq (3.3) to obtain

We have

We can identify using the provided initial data $E^0 = 0$. As a result, the inequality is expressed as

Then,

Hence, Eq (3.3) is convergent in both the Crank-Nicolson scheme and in the fully implicit method, as $\Delta_x, \Delta_t \rightarrow 0.$

Figures 6 and 7 illustrate the behavior of the analytical and numerical solutions of Eq (2.22). The parameter values are taken as $\alpha = 1.2$, $\gamma = -1.2$, $k = 0.5$, $\nu = 1.8$, $\vartheta = -1.2$, $\epsilon = 0.1$, $b = 1.5$, and $x = -50 \rightarrow 20$. The numerical solution developed for Eq (2.23) aligns with the exact solution of the equation. Figures 8 and 9 demonstrate that the exact and numerical solutions for Eq (2.24) behave similarly when $\alpha = 1$, $\gamma = 2$, $\nu = -0.05$, $\vartheta = 1$, $\epsilon = -1.5$, $b = 0.2$, and $x = -50 \rightarrow 10$.

4.   The PMA equation

To find an approximate solution to a PDE numerically, we divide the domain into small parts or nodes, creating a mesh. This mesh should be fine enough to capture areas of the solution where the values change rapidly, such as steep fronts or shock-like structures. The step size, $\Delta_x$, representing the distance between adjacent nodes, should be small enough to resolve these features accurately across the entire domain. While this method is effective, it comes at the cost of intensive computation and high expense. Additionally, it only applies to problems that do not change over time. If the solution changes with time, the mesh must constantly be redistributed. Therefore, a more efficient and alternative approach is needed to reduce the error in regions where the solution has a high curvature. In the PMA method, mesh points are redistributed in time to adapt to changes in the solution. This involves generating moving mesh nodes using the potential function's gradient as the mesh generator. The PMA equation has the advantage of requiring one fewer equation to solve for the mesh in two dimensions, and the resulting meshes are typically regular and free of tangles compared with the moving mesh partial differential equation (MMPDE) method. The MMPDE method is another technique for generating meshes in numerical simulations, particularly for problems involving moving interfaces or regions with high solution gradients. To examine the numerical solution of Eq (3.2) on the moving mesh, the PMA approach is applied. We will now explore a new transformation

The physical and computational coordinates are represented by $x$ and $\eta$, respectively. The solution for $u$ is derived from this:

We now divide the physical domain into equal subintervals as follows:

Therefore, the moving mesh is rewritten as

To obtain the physical coordinate $x$, we can use the gradient of the mesh potential $P(\eta, t)$, which is achieved using the one-dimensional PMA mesh equation:

The monitor function is known as $\Lambda(P_\eta, t)$. Additionally, $\tau$ indicates a relaxation parameter. When the value of $\tau$ is smaller, the mesh responds more quickly to changes in $\Lambda(x, t)$. Conversely, the mesh moves more slowly when the value of $\tau$ is extensive. $\rho \in \mathbb{R}$ is a positive smoothing parameter. It is subject to the boundary condition given by

and the initial condition is chosen to be

where $\eta_m \in [0, 1]$ is the computational coordinate. We get the following using the chain rule:

Thus, the expression of Eq (3.2) is given by

where $m = 1, \cdots, N_x+1$. The boundary conditions

and the initial condition are determined based on the solution of Eq (1.1) at $t = 0$. To obtain the coordinate transformation $x(\eta)$, we use the 1D PMA mesh equation

where $\Delta_\eta$ represents the step size of the computational coordinate.

4.1.   Discretization

The computational coordinate $\eta$, in this problem, is determined by

where

Then, the physical coordinate $x$ is defined by $x_m = x(\eta_m, t)$, where the boundary grids are forced to be $x_1 = a$ and $x_{N_x+1} = b$. Thus, the location of the grids $x_m$ is determined as follows:

$P_m = P(\eta_m, t)$ is the mesh potential obtained by solving Eq (4.10). Therefore, the semi-discretisation of Eq (4.9) is defined by

For a successful moving mesh method, it is essential to select a monitor function carefully. This function helps increase mesh density in areas of high interest in the PDE solution, such as propagating fronts or singularities. The monitor function should be stable under numerical discretizations, meaning that the numerical solution should not exhibit excessive oscillations or other problems. The monitor function should be computationally efficient to evaluate, as it is typically called at each iteration of the PMA algorithm. Furthermore, using two monitor functions together is common when employing adaptive methods to solve NPDEs ^[29].

For example, one monitor function could track the error in the solution, while another could track the smoothness of the solution. The grid would then be refined in areas where the error is significant or the smoothness is poor. In general, using two monitor functions can be a valuable tool for improving the performance of adaptive methods. By tracking different aspects of the solution, the adaptive method can ensure that the grid is refined in the areas most needed, leading to more accurate and efficient solutions.

Curvature monitor function:

where $\hat{\rho}$ is the user-specified parameter. If $u$ is not smooth, the discretized monitor function may change abruptly and cause slow computation. It is common practice in most standard moving mesh methods to smooth the monitor function to achieve a smoother mesh. A simple yet effective smoothing scheme involves weighted averaging as its essential ingredient, as cited in ^[29,51]. Weighted averaging can help to smooth out the monitor function and reduce oscillations in the mesh. This can lead to a more accurate and stable solution to the underlying PDE. In addition to weighted averaging, other smoothing schemes can be used for monitor functions, such as filtering and interpolation. The choice of smoothing scheme depends on the specific application and the desired properties of the mesh ^[32].

From Figures 10–12, adaptive methods are crucial in resolving intricate spatial scales and rapidly changing temporal behaviors in various applications. The fundamental idea behind the PMA equation method is to increase the number of points in regions with higher curvature. An important aspect of implementing this method is the transformation between the computational and physical domains.

Finding a suitable monitor function that works well with the underlying PDE is a crucial aspect of the success of a PMA method. Additionally, it is important to mention that the mesh redistribution is carried out by utilizing the monitor function, which oversees the mesh evaluation. The choice of monitor function can also significantly impact the PMA algorithm's convergence rate. Sometimes, it may be necessary to experiment with different monitor functions to find one that works well for a particular problem.

The choice of monitor function is ultimately a problem-specific decision. However, by understanding the key properties of a good monitor function, researchers can make informed choices that lead to more accurate and efficient solutions to their problems. It is helpful to understand the concept of a non-uniform mesh and how it can be used to improve the accuracy of numerical simulations. Figure 12 shows the associated time development of a non-uniform mesh. The mesh is initially uniform, but it becomes non-uniform over time. The solution changing rapidly between $-30$ and $-10$ highlights the importance of refining the mesh in areas with significant variation. This refinement allows the simulation to capture the details of the solution more accurately, which can be crucial for problems with steep gradients or sharp transitions.

5.   Results and discussion

We have successfully applied two analytical methods to extract traveling wave solutions of the Gilson-Pickering equation. We use the F-expansion method and the ($1/G^\prime$)-expansion strategy. Figure 4 presents the behavior of $u_9(x, t)$ when we change the values of $\gamma$ and $\alpha$, while the other parameters take fixed values $k = 1$, $\vartheta = -0.2$, $\nu = 0.01$, $\epsilon = 3.1$, and $b = 1$. Figure 5 illustrates the wave behavior of $u_{10}(x, t)$ with the values of parameters $k = 0.1$, $\vartheta = 1.01$, $\nu = -0.3$, $b = 2$, and $\epsilon = - 0.8$. As a result, these plots show that the value of $\alpha, \gamma$ affect the direction and the amplitude of the wave. Figures 4 and 5 are displayed with $N_x = 1000$ and at time $t = 1$. We verify the solutions with numerical results using two numerical techniques. First, we convert the underlying problems into a system of ODEs using the finite difference method while keeping time derivatives continuous. Then, we solve the resulting system of ODEs using the differential-algebraic system problem solver kit (DASPK) and the PMA approach, where the PMA technique reduces the error in the region where the solution has a significant variation. This method gives reliable and robust results. This can be seen in the graphical comparisons in the figures mentioned above.

In Figure 13, the numerical solutions for $u_9(x, t)$ differ from the exact solution for different values of $N_x$. It is evident from the figure that as $N_x$ becomes large, the numerical solutions of $u_9(x, t)$ converge towards the exact solutions. Specifically, there is a notable increase in error when employing $N_x = 200$. However, the numerical solutions closely approximate the exact solution (depicted by the blue line) as $N_x$ is increased to $1600$.

From Figure 14 and Table 1, the error of the numerical solution for $u_{11}(x, t)$ decreases as the value of $\Delta_x$ approaches zero, indicating a convergence towards zero error when we increase the number of points $N_x$.

6.   Conclusions

This article has concentrated on developing the exact and numerical solutions of the Gilson-Pickering equation, utilizing the F-expansion and the ($1/G^\prime$)-expansion methods to extract the analytical solutions. Moreover, we used the moving mesh method, the PMA approach for the Gilson-Pickering equation, and the implicit finite difference approach on a uniform mesh with the examination of the stability for two implicit scheme: the fully implicit method and the Crank-Nicolson method. We provided that the numerical schemes are unconditionally stable. We have successfully computed both the analytical and numerical solutions. In this work, it has been demonstrated that the solutions generated using the ($1/G^\prime$)-expansion method are novel compared to the ones determined in ^[40]. Additionally, the solutions derived through the F-expansion method are demonstrated to be more inclusive than previously noted in the literature. For example, $\beta$ can take any value and is not specified by a value. We utilized Mathematica 13.2 software to validate the correctness of our solutions by plugging them back into the original equation. One of the most excellent ways to ensure accurate solutions is to check the conformity of exact and numerical solutions. While some experts focus exclusively on obtaining exact solutions, this study employs a comparative approach, examining exact and numerical solutions to guarantee accuracy and correctness. The methodologies employed in this work can be applied to other NLPDEs encountered in the natural sciences.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors extend their appreciation to the Deanship of Research and Innovation at the Ministry of Education in Saudi Arabia for funding this research through project number 445-9-746.

Conflict of interest

The authors declare that they have no potential conflicts of interest in this article.