New Generalized Jacobi Galerkin Operational Matrices of Derivatives: An Algorithm for Solving Multi-Term Variable-Order Time-Fractional Diffusion-Wave Equations

: The current study discusses a novel approach for numerically solving MTVO-TFDWEs under various conditions, such as IBCs and DBCs. It uses a class of GSJPs that satisfy the given conditions (IBCs or DBCs). One of the important parts of our method is establishing OMs for Ods and VOFDs of GSJPs. The second part is using the SCM by utilizing these OMs. This algorithm enables the extraction of precision and efficacy in numerical solutions. We provide theoretical assurances of the treatment’s efficacy by validating its convergent and error investigations. Four examples are offered to clarify the approach’s practicability and precision; in each one, the IBCs and DBCs are considered. The findings are compared to those of preceding studies, verifying that our treatment is more effective and precise than that of its competitors.


Introduction
Fractional calculus has emerged as a subject of great interest among researchers from diverse fields in recent decades.The universal perspective that fractional operators offer for comprehending system dynamics is what primarily motivates this interest.Fractional derivatives offer a more precise and comprehensive description of various physical phenomena compared to traditional integer-order derivatives [1][2][3][4].Consequently, there exists a wide range of definitions for fractional differentiation in the literature (for further details, refer to [5][6][7]).A list of abbreviations used in the paper is given in Abbreviations Section.
Researchers have made worthy progress in the field of fractional calculus, emphasizing the extension of the framework to encompass VOFDs.This extension facilitates an understanding of diverse dynamic systems.The authors in [8] conducted a specific study in which they investigated the characteristics of VOFD operators [9,10].
VOFC is a strong framework that obtains the nonlocal properties of different systems very well [11][12][13][14][15][16][17][18].For instance, in study [19], the utilization of VOFD operators enabled the modeling of the microscopic structure of materials.As shown in an additional publication [20], continuum elasticity also made use of the Riesz-Caputo fractional derivative of space-dependent order.In [21], the authors recommend using collocation and tau spectral techniques to solve the time-fractional heat equation numerically.Studies [22,23] specifically focused on understanding and analyzing how the fractional order, which characterizes the system's dynamics, evolves over time and influences its viscoelastic properties.By examining the time-dependent changes in the fractional order, valuable insights were gained into the complex behavior of such systems.Overall, VOFC has proven to be a valuable tool in engineering mechanics, providing a versatile approach to accurately describe and analyze a wide range of dynamic systems.
The pursuit of analytical solutions for FDEs poses a significant challenge, often necessitating reliance on numerical approximations.As a result, various numerical techniques have addressed the complexities associated with FDEs, enabling researchers to effectively tackle problems that would otherwise be challenging to solve analytically.Previous research has witnessed the utilization of a variety of approaches to create numerical solutions for FDEs through the use of both orthogonal and non-orthogonal polynomials.There are a number of papers [24][25][26][27][28][29][30] that talk about different OMs of JPs and some of their special cases.These papers use different spectral methods by utilizing these OMs to solve different kinds of DEs and FDEs numerically, subject to different kinds of IBCs.Additionally, researchers in [31][32][33][34][35][36] follow the same methodology to solve VOFDEs numerically.Additionally, in [37,38], Bernstein polynomials were employed to approximate solutions for FDEs.Another notable contribution in [39] proposed a numerical scheme based on Fourier analysis for solving FDEs.Furthermore, finite difference approximations were discussed in [40] as a means to construct numerical schemes for FDEs.In [41], the author uses a class of modified JPs to introduce a novel method for numerically solving MTVO-FDEs with initial conditions.
Many researchers [42][43][44] have shown that JPs have properties that make them very useful for solving different kinds of DEs, especially when spectral methods are used.Some of these properties are that they are orthogonal, have exponential accuracy, and have two parameters that allow us to shape approximate solutions in different ways.These inherent properties of JPs make them highly suitable for effectively solving a wide range of diverse problems encountered in various fields of study.
It is worth noting that a significant class of fractional partial differential equations that has received considerable attention in recent years is the TFDWE.This equation arises from the classical diffusion-wave equation by replacing the second-order time derivative term with a fractional derivative of order 1 < α < 2 [45].The TFDWE describes important physical effects seen in many different types of systems, such as colloidal, amorphous, glassy, and porous materials, as well as dielectrics and semiconductors, comb structures, polymers, random and disordered media, biological systems, and geophysical and geological processes (see [46] and references therein).
Furthermore, it is worth mentioning that the TFDWE serves as an accurate model for many universal electromagnetic, acoustic, and mechanical responses [47,48].A single TFDWE may not be able to fully describe the underlying processes in some real-life situations.This is why an MTFDWE was created, as shown in [46,49].This MTFDWE formulation offers a more comprehensive representation of complex systems and their dynamic behavior, allowing for more accurate modeling of the underlying processes.
We considered the general form of MTVO-TFDE in our investigation as follows: subject to the IBCs: or the DBCs: where (i = 1, 2, . . ., m) are the VOFDs defined in the Caputo sense, as given in Section 2. The significance and difficulty of proving the existence, unicity, and dependency of parameters should be emphasized.The significance of determination of ν and ν j in (1) is evident in theory and practice due to the strong relationship between them and the heterogeneity and associated physical features of media [50].For instance, in the situation of a single term case, Cheng et al. [51] first demonstrated the uniqueness of ν using the boundary condition data that were provided.
We present a new Galerkin OM for Ods and new OMs for VOFDs of GSJPs in the sense of Caputo.This was made in order to find a new way to solve the problem shown by (1) and the conditions (2) or (3).These operational matrices are specifically tailored for the basis vectors of GSJPs.Leveraging these OMs, we have established a powerful tool that enables the accurate computation of numerical solutions using the SCM to solve a wide range of MTVO-TFDEs.This novel method opens up new ways to solve this type of FDE numerically and more effectively.
To summarize, the main article's contributions are as follows: (i) We introduce two classes of GSJPs to satisfy the given IBCs and DBCs (see Section 3.2).
(ii) We establish Galerkin OMs for the Ods and for VOFDs of the introduced GSJPs in the sense of Caputo (see Sections 4 and 5).(iii) We address the presented MTVO-TFDE using the proposed GSJPs and their constructed OMs in conjunction with the SCM (see Section 6).(iv) We present a study of convergence and error analysis for the numerical solution obtained through the proposed scheme (see Section 7).
The paper has the following outline: Section 2 provides comprehensive coverage of the needed concepts of VOFC.Section 3 highlights the necessary attributes of JPs and GSJPs.Sections 4 and 5 emphasize the development of new Galerkin OMs for Ods and VOFDs of GSJPs.These OMs are intended to address Equation (1) when IBCs (2) or DBCs (3) are considered.Section 6 explores the selection of newly generated OMs in the SCM to address the aforementioned issue.Section 7 analyzes convergence and error estimation.Section 8 includes four examples to clarify the approach's practicability and precision, and the findings are compared to those of preceding studies, verifying that our treatment is more effective and precise than its competitors.Finally, Section 9 summarizes the key outcomes, implications derived from our investigation, and the scope of future work.

Basic Definition of Caputo VOFDs
This section serves to introduce and discuss important definitions and necessary assets that lay the groundwork for the development of our proposed technique.These tools are crucial building blocks that provide the necessary foundation for effectively handling the MTVO-TFDE under consideration.Definition 1 ([39,45]).Suppose h(x, t) is a twice differentiable function.The Caputo derivative operator with the variable-order ν(x, t) is defined as: The Caputo VOFD exhibits the characteristics: Remark 1.The reader who is interested may find several definitions and additional attributes of VOFDs in [4] (pp. [35][36][37][38][39][40][41][42].

An Overview of the Shifted JPs and Their Generalized Ones
Introducing the fundamental features of the JPs and their shifting form is the main objective of this section.Additionally, a group of GSJPs is introduced.

An Overview of the Shifted JPs
The orthogonal JPs, J (α, β) n (x), α, β > −1, satisfy the following relationship [52]: The shifted JPs, denoted as J , are in accordance with: where w α, β The fundamental expansions that will be used in this paper are [53] (Section 11.3.4): 1.
The power form representations of J (α, β) ℓ,n (t) are as follows: where Alternatively, the expressions for t k and (ℓ − t) k in relation to J (α, β) ℓ,r (t) have the forms: and b

Introducing GSJPs
In this section, it is helpful to talk about two polynomials, ϕ T ,j (t) j≥0 and ψ (α, β) ℓ,j (t) j≥0 .These are needed to satisfy the homogeneous form of the given IBCs (2) and DBCs (3): Subsequently, these polynomials satisfy the orthogonality relations as follows: where w α, β

Two OMs for Ods and VOFDs of ϕ
T ,j (t) In this part, we introduce two OMs related to ODEs and VOFDs of ϕ (α, β) T ,j (t).To facilitate this, we commence with Theorem 1: T ,i (t)) for all i ≥ 0 can be computed as follows: where ϵ n,i (t) = n i! (−1) i ( β + 1) i t n−1 , and Then, the two desired OMs of can be computed as follows: Corollary 1 ([41]).The general derivative of Φ (α, β) T ,N (t) has the form: with η (m) T ,i (t) for all i ≥ 0, has the form: which leads to: where D (ν(x,t)) = (d i,j (ν(x, t))) is a matrix of order (N + 1) × (N + 1), which can be expressed explicitly as: where and Proof.By considering (7) and applying (6b), we obtain: By utilizing (9), expanding and collecting similar terms, and using useful computational formulae of the Pochhammer symbol and the Gamma function (see [54] (p.758)), it is possible to convert (24) to (19), which is represented as follows: and this representation enables us to declare that ( 20) is proved and that the theorem is completely proved.□ As an application of Theorem 2, for N = 4, α = − β = 1/2, and ν(x, t) = x t, the OM D (ν(x,t)) has the form:

Two OMs for Ods and VOFDs of ψ (α, β)
ℓ,i (t) In this part, we establish two OMs for Ods and VOFDs of ψ (α, β) ℓ,i (t).In order to achieve this, we commence by proving the subsequent lemma: Proof.We have In view of ( 7) and ( 9), we obtain: consequently, through expansion, the accumulation of similar terms, and using useful computational formulae of the Pochhammer symbol and the Gamma function (see [54] (p.758)), it becomes evident that: Substitution of (31) and J into (29) leads to (27), and then the lemma is proved completely.□ ℓ,i (t), for all i ≥ 0, has the form: where Proof.We have Using Theorem 1 and Lemma 1 leads to (32), and the proof of the theorem is complete.

□
Then, the desired two OMs of can be computed as follows: with ξ where ) is a matrix of order (N + 1) × (N + 1), which can be expressed explicitly as: where the matrix elements are defined as follows: and the coefficients Λ i,j (ν(x, t)) have the form where and the functions ε ℓ,i (x, t) have the form Proof.By considering (7) and applying (6b), it is easy to see that which may be expressed as: where and Using (9), through expansion, the accumulation of similar terms, and some algebraic manipulation (see [54] (p.758)), it becomes evident that Also, substitution of the form c s in (49)-after some simple algebraic manipulationgives Now, it remains that we show that i+1 = 0-after some algebraic manipulation-leads to: By expanding and collecting similar terms, we obtain: where Again, using (9), expanding and collecting similar terms leads to: Using the expressions of c s+1 and b (r) j after some algebraic manipulation leads to: then In view of Equations ( 46), ( 50), (51), and ( 58), the proof of ( 37) is complete, which can be expressed as follows: this gives (38), and the theorem is completely proved.□ As an application of Corollary 2 and Theorem 4, for N = 4, α = − β = 1/2, and ν(x, t) = x t, the OMs H and D(ν(x,t)) have the forms: and (61)
We can compute a i,j (where i, j = 0, 1, . . ., N) by using the right solver to work through a set of (N + 1) 2 Equation (68).Achieving the target numerical solution relies heavily on these coefficients.

Non-Homogeneous IBCs and DBCs
Changing Equation (1) along with the non-homogeneous conditions (2) or (3) into a similar form with homogeneous conditions is an important part of creating the proposed algorithm.The following transformation accomplishes this change: u(x, t) = y(x, t) − q(x, t), q(x, t) = q 1 (x, t) + q 2 (x, t), where the two functions q 1 (x, t) and q 2 (x, t) are defined as follows: In the Case of IBCs: In the Case of DBCs: As a result, the current issue may be simplified by solving the following updated equation: where subject to the homogeneous IBCs or the homogeneous DBCs Then, Remark 2. Section 8 explains the algorithm that was used to solve various numerical examples.An Intel®Core™ i9-10850 CPU running at 3.60 GHz, with 10 cores and 20 logical processors and 64.0 GB RAM, was used to do the calculations on a computer system equipped with Mathematica 13.3.1.0with GSJCOPMM, and the following algorithmic steps may be described to solve the MTVO-TFDWE (Algorithms 1 and 2): Algorithm 1 GSJCOPMM algorithm to solve (1) subject to IBCs.
Use Mathematica's built-in numerical solver to solve the system obtained in [Output 5].Stage 6.

Stage 3.
Calculate the matrices: Define R N (x, t) as in Equation (67).
Use Mathematica's built-in numerical solver to solve the system obtained in [Output 5].Stage 6.

Convergence and Error Analysis
Here, we look at the suggested method's convergence and error estimations.The space S N , defined as follows, is our primary area of interest, T ,j (t) : i, j = 0, 1, . . ., N , forDBCs.
Additionally, the error between y(x, t) and its approximation y N (x, t) can be defined by The numerical scheme's error is examined in the paper via the use of the L ∞ norm error estimation In the following, Theorem 5 discusses the error estimate when the solution satisfies IBCs, while Theorem 6 discusses it when DBCs are considered.
Proof.Following the same procedures in the proof of Theorem 5, we obtain (91).□ The next corollary demonstrates how quickly the resulting errors are convergent.
Corollary 3.For all N > q − 1, the following estimate holds: The stability of error, or the process of estimating the propagation of error, is the focus of the subsequent theorem.Theorem 7. Given any two iterative estimates of y(x, t), the result is: Proof.We have By considering (92), we can obtain (93).□ Remark 3.While relative error can be a useful measure in certain situations, the choice to utilize absolute error in our study was driven by considerations of stability, consistency with established literature, practical interpretation, and the characteristics of the problem at hand.We believe that using absolute error allows for a reliable and meaningful assessment of the accuracy of our proposed method in solving MTVO-TFDWEs.

Numerical Simulations
In order to show how the approach presented in Section 6 works and how efficient it is, four examples are provided in this section.The assessment is based on MAE between the precise and approximate solutions.We show in Example 1 that, for certain simple cases of the functions ν(x, t) and ν j (x, t), j = 1, 2, . . ., m, we can use GSJCOPMM to find the precise solution to the given numerical problem, whether it involves IBCs or DBCs, and that it has a polynomial solution of degree N.This solution can be found by combining T ,j (t), i, j = 0, 1, . . ., N − 2. Otherwise, some numerical solutions are obtained with high accuracy.Additionally, the choice of examples in our study is motivated by the desire to cover different conditions, compare with existing results, showcase accuracy and efficiency, and demonstrate practical relevance.These examples collectively provide a comprehensive assessment of the performance and applicability of our proposed method for solving MTVO-TFDWEs.
In addition, Tables 1-4 display the calculated errors that were obtained to produce numerical solutions y N (x, t) using GSJCOPMM with N = 0, 1, . . ., 14.In these tables, excellent computational results are obtained.The comparisons between GSJCOPMM and other techniques in [57,58] are presented in Tables 5 and 6.The tables show that, compared to the other methods, GSJCOPMM gives the most accurate results.Furthermore, as can be seen in Figures 1, 2, 3, 4a, 5a, 6, 7, 8a, 9a, and 10, the precise and approximate solutions in Examples 1-4 have a high level of agreement.The absolute and log errors in Figures 4b, 5b, and 9 serve to illustrate the convergence and stability of the solutions to the given Problems 2 and 3 when applying GSJCOPMM and using different N, α, and β.The choice of parameter combinations enables us to fine-tune our method, which leads to the most accurate results.Fractal Fract.2024, 1, 0 19 of 26 subject to IBCs: Fractal Fract.2024, 1, 0 19 of 26 subject to IBCs: subject to IBCs: y(x, 0) = 0, y t (x, 0) = 0, y(0, t) = 0, y(1, t) = 0, or DBCs: y(x, 0) = 0, y(x, 1) = 0, y(0, t) = 0, y(1, t) = 0, (96) where g(x, t) is selected such that the solution to (94) is y(x, t) = 100 The application of GSJCOPMM to obtain approximated solutions with the two cases of IBCs (95) and DBCs (96) using N = 2, . . ., 7, ν(x, t) = 2 − 0.3 e −x t and ν 1 (x, t) = 2 − 0.6 e −x t gives acceptable accuracy, as shown in Table 1. Figure 1a,c demonstrates that the obtained solutions achieve an accuracy of 10 −21 at N = 7, aligning with the exact solution depicted in Figure 1b.Remark 4. It is worth noting that, although the precise solution to (94) is a polynomial, the obtained approximate solutions did not give this precise solution; this is because of the complexity form of the two functions ν(x, t), ν 1 (x, t).However, in other special cases, for instance, ν(x, t) = x t, ν 1 (x, t) = 2x t, for all available values of α and β, the precise solution is obtained using N = 2, where the expansion coefficients c i,j = 0, i, j = 0, 1, 2, have the form: Remark 5.Although errors smaller than 10 −16 may not have direct practical significance, increasing the value of N can still meaningfully improve accuracy in numerical computations.The choice of N depends on the specific requirements and desired level of accuracy for the problem under consideration.We emphasize that the relative improvement in accuracy should be considered when evaluating the meaningfulness of different values of N.
or DBCs: y(x, 0) = sin(x), y(x, 1) = sin(x − 1), y(0, t) = − sin(t), y(1, t) = sin(1 − t), where g(x, t) is selected such that the solution to (100) is y(x, t) = sin(x − t).The application of GSJCOPMM to obtain approximated solutions to (100) subject to IBCs (101) or DBCs (102), respectively, using N = 2, 4, 6, 8, 10, 12, ν(x, t) = 1.8 + 0.2 sin(x t) and ν j (x, t) = 1.8 − (0.1)j + 0.2 sin (x t), j = 1, 2, 3, 4, 5, gives acceptable accuracy for the solutions obtained, as shown in Table 3 and Figure 7a,c.These solutions reach an agreement of an accuracy of 10 −14 at N = 12, aligning with the approximate and exact solutions depicted in Figure 6.Again, the heat map graphs displayed in Figure 7b,d and the error graphs in Figures 8a and 9a confirm the effectiveness of GSJCOPMM.Remark 6.In view of the presented CPU time (in seconds), our approach has efficient performance.The calculations show that the memory consumption was excellent.For example, the calculated CPU time using N = 10 is 6% slower than N = 7 and, moreover, requires increasing by 20% the memory consumption of RAM compared to the N = 7 calculation.The numerical examples and comparisons provided in our paper highlight the superior accuracy and efficiency of our algorithm, solidifying its potential for solving MTVO-TFDEs effectively.When we compared the resource use of our method to that described in [57,58], we saw that those papers did not provide CPU time and memory usage.However, based on our analysis, our approach demonstrates better performance compared to the referenced methods.

Conclusions
This study presents a generalized form of shifted JPs that fulfill homogeneous IBCs or DBCs.Then, by making use of the OMs derived in Sections 4 and 5 with the SCM, an approximation algorithm for the given MTVO-TFDWE is established.Three examples, including MTVO-TFDWE (1), were tested using the suggested technique, GSJCOPMM, to prove its high accuracy and efficiency.The examples were subjected to either the IBCs (2) or the DBCs (3).It would be great to see our results generalized to other kinds of initial and boundary conditions.Studying the system's behavior in different settings would provide interesting insights and make our conclusions even more applicable.Additionally, the theoretical results obtained in this study can be developed to deal with more complex versions of MTVO-TFDWEs, such as multi-dimensional multi-term Caputo's time-fractional mixed sub-diffusion and diffusion-wave equations.Additionally, incorporating adaptive strategies and parallel computing techniques could further enhance the efficiency and scalability of the algorithm.Overall, this research contributes to the advancement of numerical methods for MTVO-TFDWEs and opens up avenues for exploring their applications in various fields.
Funding: This research received no external funding.

Fractal
Fract.2024, 1, 0 16 of 26 solutions in Examples 1-4 have a high level of agreement.The absolute and log errors in Figures 4b, 5b, and 9 serve to illustrate the convergence and stability of the solutions to the given Problems 2 and 3 when applying GSJCOPMM and using different N, α, and β.The choice of parameter combinations enables us to fine-tune our method, which leads to the most accurate results.

Fractal
E 12 (x, t) associated with IBCs.(b) Heat map graph of E 12 (x, t) with IBCs.(c) E 12 (x, t) associated with DBCs.(d) Heat map graph of E 12 (x, t) with DBCs.
Absolute errors at t = 0.1.Graph of Log 10 E N against N.
Graph of Log 10 E N against N.

E 14 (
Absolute errors at t = 0.3 with IBCs.(b) Heat map graph of E 14 (x, t) with IBCs.

E 14 (
Absolute errors at t = 0.3 with DBCs.(d) Heat map graph of E 14 (x, t) with DBCs.

E 14 (
Absolute errors at t = 0.3 with IBCs.(b) Heat map graph of E 14 (x, t) with IBCs.

E 14 (
Absolute errors at t = 0.3 with DBCs.(d) Heat map graph of E 14 (x, t) with DBCs.