 # NTIM solution of the fractional order parabolic partial differential equations

• Nasir Ali , Rashid Nawaz , Anwar Saeed , Taza Gul and
From the journal Open Physics

## Abstract

In this article, natural transform iterative method has been used to find the approximate solution of fractional order parabolic partial differential equations of multi-dimensions together with initial and boundary conditions. The method is applicable without any discretization or linearization. Three problems have been taken as test examples and the results are summarized through plots and tables to show the efficiency and reliability of the method. By practice of a few iterations, we observe that the approximate solution of the parabolic equations converges to the exact solution. The fractional derivatives are considered in the Caputo’s sense.

## 1 Introduction

Fractional calculus was first studied in the seventeenth century and has recently received a lot of interest. Scientists have discovered that fractional calculus may describe memory and hereditary features of various problems in science and engineering due to fractional-order derivatives. As a result, we observe the fractional calculus in many domains like signal processing, diffusion, physics, fluid mechanics, biology, chemistry, economics, polymer rheology etc. It has its significance almost in every field of science and technology. It models natural phenomena in more suitable way than classical calculus. Due to the rapid development in the day-to-day life, fractional calculus has also played important role in engineering, biosciences, and finance [1,2,3,4]. Fractional calculus is the generalization of classical calculus. The differential equations arising in fractional calculus are termed as fractional differential equations (FDEs). FDEs have a non-integer order derivative and can be solved through derivative and integral operators related to fractional calculus. Different operators have been defined by mathematicians for the solution of FDEs [5,6,7].

In many instances exact solution to differential equations is not always possible and the approximate solutions to these equations are obtained by using different numerical and analytical methods. In order to obtain approximate solutions to differential equations, different numerical methods were developed over time. However, the numerical solutions were not enough to determine the overall properties of certain systems of differential equations which leads us to the development of some new analytical and semi-analytical methods. These methodologies have revolutionized numerical analysis, allowing us to present difficult issues with both qualitative and quantitative analysis. For the solution of linear and non-linear FDEs, several analytical, numerical, and homotopy-based approaches have been used [8,9,10,11]. Some analytical techniques in combination with transformations have also been applied to handle FDEs more suitably. The homotopy analysis Sumudu transform Method is a combination of homotopy analysis method (HAM) and Sumudu transform used by Singh et al. for handling the fractional Caudrey–Dodd–Gibbon equations . Dubey et al. used the local fractional natural HAM which is the combination of the HAM and local fractional natural transform for solving partial differential equations (PDEs) of fractional order . Supriya Yadav et al. also used the q-HASTM for solving the fractional reaction-diffusion equations . Recently Jagdev Singh presented the composite fractional derivative to analyze the fractional blood alcohol model . The derivative is considered in the Caputo’s sense due to its most applicability and popularity to the FDEs and can be handled easily by the proposed method.

One of the relevant algorithms we have used in this research is the natural transform iterative method (NTIM), which is based on new iterative method (NIM) and the natural transform [16,17]. The usage of natural transform and NIM make it easier and more appropriate to handle FDEs. For the investigation of PDEs of integer order and of fractional order, this method is free of discretization and linearization. NIM has been used by a number of other researchers to solve fractional order PDEs [18,19,20,21]. The proposed method’s convergence may be demonstrated, as demonstrated by Bhalekar et al. . This work looks into the fractional order parabolic differential equations of fourth-order with variable coefficients. Fractional order parabolic PDEs have the general form as reported in ref. .

(1) D τ β c ζ + μ ( η , υ , ψ ) ζ η η η η + λ ( η , υ , ψ ) ζ υ υ υ υ + δ ( η , υ , ψ ) ζ ψ ψ ψ ψ = 0 , 1 < β 2 , τ 0 ,

where ζ = ζ ( η , υ , ψ , τ ) , and μ , λ , δ are the variable functions of ( η , υ , ψ ) . Eq. (1) is subjected to the initial conditions given as

(2) ζ ( η , υ , ψ , 0 ) = f 1 ( η , υ , ψ , 0 ) , ζ τ ( η , υ , ψ , 0 ) = f 2 ( η , υ , ψ ) .

and boundary conditions as

(3) ζ ( a , υ , ψ , τ ) = g 1 ( υ , ψ , τ ) , ζ τ ( b , υ , ψ , τ ) = g 2 ( υ , ψ , τ ) ,

(4) ζ ( η , a , ψ , τ ) = k 1 ( η , ψ , τ ) , ζ τ ( η , b , ψ , τ ) = k 2 ( η , ψ , τ ) ,

(5) ζ ( η , υ , a , τ ) = h 1 ( η , υ , τ ) , ζ τ ( η , υ , b , τ ) = h 2 ( η , υ , τ ) ,

(6) ζ η η ( a , υ , ψ , τ ) = g ¯ 1 ( υ , ψ , τ ) , ζ η η ( b , υ , ψ , τ ) = g ¯ 2 ( υ , ψ , τ ) ,

(7) ζ υ υ ( η , a , ψ , τ ) = k ¯ 1 ( η , ψ , τ ) , ζ υ υ ( η , b , ψ , τ ) = k ¯ 2 ( η , ψ , τ ) ,

(8) ζ ψ ψ ( η , υ , a , τ ) = h ¯ 1 ( η , υ , τ ) , ζ ψ ψ ( η , υ , b , τ ) = h ¯ 2 ( η , υ , τ ) ,

where f i , g i , k i and h i , i = 1 , 2 are continuous functions of ( η , υ , ψ ) . The fourth-order initial and boundary value problems have their significance in many fields of sciences. For example, fourth-order PDEs are used to model floor systems, bridge slabs, aviation wings, and window glasses [24,25]. A lot of problems in the fields of electrochemistry, electromagnetics, material science, diffusion processes, and chaotic dynamics can be modeled by FDEs [26,27,28]. The parabolic PDEs play a vital role in the study of viscoelastic and inelastic flows, beam deformation, and layer deflection [29,30]. Many researchers in the literature have solved the parabolic fourth-order PDEs using various methods in the literature [31,32,33]. In this work, we look at non-perturbation analytical approaches for solving fractional order PDEs. The goal is to solve FDEs with analytical approaches that provides accuracy, convergence, and computational efficiency.

The rest of this article is set out as follows. The first section contains some basic fractional calculus definitions. The basic concept of NTIM is explained in Section 2. The application of NTIM to parabolic equations is covered in Section 3. In Section 4, some numerical results are presented. Finally, Section 5 gives the conclusion.

## 2 Fractional calculus

Some definitions are presented from the fractional calculus.

### 2.1 Definition

The fractional integral in Riemann–Liouville’s (R–L) sense of a function f(ϕ) is defined as

I x β = 1 Γ ( β ) 0 x ( x ϕ ) β 1 f ( ϕ ) d γ , if β > 0, t > 0 , f ( ϕ ) , if β = 0,

where Γ is the gamma function defined as

Γ ( w ) = 0 e x x w 1 d x , w .

### 2.2 Definition

The fractional order derivative of a function f (ϕ) in the Caputo’s sense is given as

for n N , t 1 , x > 0 and ϕ t :

D t β φ ( x ) = I n β n t n f ( ϕ ) , if n 1 < β n , n N , d β d t β ( f ( x ) ) , if β N .

### 2.3 Definition

Relationship of the Caputo’s derivative and R–L integral is given as

For n N , and n 1 < β n

I x β [ D x β φ ( χ ) ] = φ ( χ ) + k = 0 n 1 φ ( m ) ( γ ) ( χ γ ) k Γ ( k + 1 ) , x > 0 .

### 2.4 Definition

Natural transform of ϕ ( t ) is 

+ ( ϕ ( τ ) ) = ( s , u ) = 1 u 0 e s τ u ( ϕ ( τ ) ) d τ , s , u > 0 ,

u and s are the transformation variables.

### 2.5 Definition

The inverse of the natural transform ( s , u ) is defined as

[ ( s , u ) ] = ϕ ( τ ) = 1 2 π i c i c + i e s τ u ( ( s , u ) ) d s ,

s = a + b i is the complex plan for executing the integral along s = c , where c R .

### 2.6 Definition

If ϕ ( τ ) is a function and ϕ n ( τ ) is its nth derivative, then natural transform of ϕ ( τ ) is

+ ( φ n ( τ ) ) = n ( s , u ) = s n u n ( s , u ) k = 0 n 1 s n ( k + 1 ) u n k ( φ n ( 0 ) ) , n 1 .

### 2.7 Theorem

If k ( τ ) and h ( τ ) are defined on a set A and have the natural transforms K(s,u) and

H(s,u), respectively, then,

[ k h ] = u K ( s , u ) H ( s , u ) ,

where [h * k] is convolution of h and k.

### 2.8 Remark

A few important natural transformations of some functions are given below.

 ℕ + [ 1 ] = 1 s ℕ + [ e k τ ] = 1 s − k u ℕ + [ τ ] = u s 2 ℕ + [ sin ( τ ) ] = u s 2 − u 2 ℕ + τ k − 1 Γ ( k ) = u k - 1 s k ℕ + [ cos ( τ ) ] = s s 2 + u 2 ℕ + τ k β Γ ( k β + 1 ) = u k β s k β + 1 ℕ + τ k β − 1 Γ ( k β ) = u k β − 1 s k β

## 3 NTIM 

Consider the FDE of the form

(9) D t β ( ζ ( y ̲ , τ ) ) = g ( y ̲ , τ ) + L ( ζ ( y ̲ , τ ) ) + N ( ζ ( y ̲ , τ ) ) , 0 y ̲ , τ , k 1 < β k ,

where D t β is the Caputo’s fractional derivative operator of order β, k ∈ N and for this work we take k = 2 and y ̲ = y 1 , y 2 , , y k . L and N represent the linear and nonlinear functions, respectively. g ( y ̲ , τ ) is a known function. The initial conditions are

(10) ζ ( y ̲ , 0 ) = f 1 ( y ̲ ) , ζ τ ( y ̲ , 0 ) = f 2 ( y ̲ ) ,

Taking the natural transform of Eq. (9), we have

(11) + [ D t β ( ζ ( y ̲ , τ ) ) ] = + [ g ( y ̲ , τ ) ] + + [ L ( ζ ( y ̲ , τ ) ) + N ( ζ ( y ̲ , τ ) ) ] ,

by applying the differentiation property of natural transform (given in definition 2.6) to Eq. (11) we have

(12) s β u β + [ ζ ( y ̲ , τ ) ] s β 1 u β ζ ( y ̲ , 0 ) s β 2 u β 1 ζ τ ( y ̲ , 0 ) = + [ g ( y ̲ , τ ) ] + + [ L ζ ( y ̲ , τ ) + N ζ ( y ̲ , τ ) ] ,

Using the initial conditions and rearranging Eq. (12) we obtain

(13) + [ φ ( y ̲ , τ ) ] = f 1 ( y ̲ ) s + u f 2 ( y ̲ ) s 2 τ + u β s β ( + [ g ( y ̲ , τ ) ] ) + u β s β ( + [ L ( φ ( y ̲ , τ ) ) + N ( φ ( y ̲ , τ ) ) ] ) ,

Where the linear term L ζ ( y ̲ , t ) can be written in the form of infinite series as

(14) L m = 0 ζ m ( y ̲ , τ ) = m = 0 L ( ζ m ( y ̲ , τ ) ) ,

and N ( ζ ( y ̲ , τ ) ) , the non-linear term is proposed as in Eq. (15)

(15) N m = 0 ζ m ( y ̲ , τ ) = N ( ζ 0 ( y ̲ , τ ) ) + m = 1 N j = 0 i ζ j ( y ̲ , τ ) N j = 0 m 1 ζ j ( y ̲ , τ ) ,

Using Eqs. (14) and (15) in Eq. (13) we obtain

(16) + i = 1 ζ i = f 1 ( y ̲ ) s + u f 2 ( y ̲ ) s 2 τ + u β s β ( + [ g ( y ̲ , τ ) ] ) + u β s β + m = 0 L ( ζ m ) + N ( ζ 0 ) + m = 1 N j = 0 m ζ j N j = 0 m 1 ζ j ,

The recursive relation of Eq. (16) by the use of natural transform is

(17) + [ ζ 0 ( y ̲ , τ ) ] = u β s β + [ g ( y ̲ , τ ) ] , + [ ζ 1 ( y ̲ , τ ) ] = u β s β + [ ( ζ 0 ) + N ( ζ 0 ) ] , + [ ζ 2 ( y ̲ , τ ) ] = u β s β + [ L ( ζ 1 ) + N ( ζ 0 + ζ 1 ) ( ζ 0 ) ] + [ ζ i + 1 ( y ̲ , τ ) ] = u β s β + [ L ( ζ i ) + ( ζ 0 + ζ 1 + + ζ i ) N ( ζ 0 + ζ 1 + + ζ i 1 ) ] , i 0 ,

Applying the inverse natural transform to Eq. (17) the solution component can be obtained as

(18) ζ 0 ( y ̲ , τ ) = f 1 ( y ̲ ) s + u f 2 ( y ̲ ) s 2 τ + u β s β + [ g ] ζ 1 ( y ̲ , τ ) = u β s β + [ L ( ζ 0 ) + N ( ζ 0 ) ] ζ 2 ( y ̲ , τ ) = u β s β + [ L ( ζ 1 ) + N ( ζ 0 + ζ 1 ) N ( ζ 0 ) ] ζ i + 1 ( y ̲ , τ ) = u β s β + [ L ( ζ i ) + N ( ζ 0 + ζ 1 + + ζ i ) N ( ζ 0 + ζ 1 + + ζ i 1 ) ] , i 0 ,

The n terms’ approximate solution of Eqs. (9) and (10) by the proposed method is obtained by adding the components as

(19) ζ ( y ̲ , τ ) = ζ 0 ( y ̲ , τ ) + ζ 1 ( y ̲ , τ ) + ζ 2 ( y ̲ , τ ) + + ζ n 1 ( y ̲ , τ ) .

### 3.1 Convergence of NTIM 

#### 3.1.1 Theorem

If N is analytic in a neighborhood of ζ 0 and

N m ( ζ 0 ) = Sup { N m ( ζ 0 ) ( b 1 , b 2 , . . b n ) / b k 1 , 1 k m } l ,

for any m and for some real l > 0 & ζ k M < 1 e k = 1 , 2 , , then the series m = 0 G m is absolutely convergent and moreover

G m l M m e m 1 ( e 1 ) , m = 1 , 2 , .

To show the boundedness of ζ k , for every 1 k m , the conditions on N ( j ) ( ζ 0 ) are given which are sufficient to guarantee convergence of the series.

Sufficient condition for convergence is given in the following theorem.

#### 3.1.2 Theorem

If N is C and N m ( ζ 0 ) M e 1 m , then the series m = 0 G m is absolutely convergent.

The above mentioned are the required conditions for the convergence of the series ζ j .

## 4 Applications

### 4.1 Problem 1

Consider the (1 + 1) dimension parabolic equation of the form 

(20) D τ β c ζ ( η , τ ) + 1 η + η 4 120 ζ η η η η ( η , τ ) = 0 , 1 < β 2 , τ 0 ,

with initial conditions

(21) ζ ( η , 0 ) = 0 , ζ τ ( η , 0 ) = 1 + η 5 120 ,

and boundary conditions as

(22) ζ 1 2 , τ = 1 + ( 1 / 2 ) 5 120 sin ( τ ) , ζ ( 1 , τ ) = 1 + 121 120 sin ( τ ) , ζ η 1 2 , τ = 1 6 1 2 3 sin ( τ ) , ζ η ( 1 , τ ) = 1 6 sin ( τ ) .

Rearranging Eq. (20) as

(23) D t β c ζ ( η , τ ) = 1 η + η 4 120 ζ η η η η ( x , τ ) ,

and applying the natural transform, we get

(24) + [ D t β c ζ ( η , τ ) ] = + 1 η + η 4 120 ζ η η η η ( η , τ ) ,

Using the differentiation property of natural transform, we have

(25) s β u β + [ D t β c ζ ( η , τ ) ] s β 1 u β ζ ( η , 0 ) s β 2 u β 1 u ζ τ ( η , 0 ) = + 1 η + η 4 120 ζ η η η η ( η , τ ) ,

which after rearranging and simplification we have

(26) + [ D t β c ζ ( η , τ ) ] = ζ ( η , 0 ) s + u ζ τ ( η , 0 ) s 2 + u β s β + 1 η + η 4 120 ζ η η η η ( η , τ ) ,

Applying the inverse natural transformation, we obtain

(27) ζ ( η , τ ) = ζ ( η , 0 ) s + u ζ τ ( η , 0 ) s 2 + u β s β + 1 η + η 4 120 ζ η η η η ( η , τ ) ,

Using the initial conditions, we have

(28) ζ ( η , τ ) = 1 + η 5 120 τ + u β s β + 1 η + η 4 120 ζ η η η η ( η , τ ) ,

Using the recursive relation, as ζ ( η , τ ) is in series form, we have

(29) ζ 0 ( η , τ ) = 1 + η 5 120 τ , ζ 1 ( η , τ ) = u β s β + 1 η + η 4 120 ζ 0 η η η η ( η , τ ) = 1 + η 5 120 τ β + 1 Γ ( β + 2 ) , ζ 2 ( η , τ ) = u β s β + 1 η + η 4 120 ζ 1 η η η η ( η , τ ) = 1 + η 5 120 τ 2 β + 1 Γ ( 2 β + 2 ) , ζ 3 ( η , τ ) = u β s β + 1 η + η 4 120 ζ 2 η η η η ( η , τ ) = 1 + η 5 120 τ 3 β + 1 Γ ( 3 β + 2 ) , ζ 4 ( x , τ ) = u β s β + 1 η + η 4 120 ζ 3 η η η η ( η , τ ) = 1 + η 5 120 τ 4 β + 1 Γ ( 4 β + 2 ) ,

Adding the solution components, we get the approximate solution as

(30) ζ ( η , τ ) = ζ 0 ( η , τ ) + ζ 1 ( η , τ ) + ζ 2 ( η , τ ) + ζ 3 ( η , τ ) + ζ 4 ( η , τ ) + = 1 + η 5 120 τ τ β + 1 Γ ( β + 2 ) + τ 2 β + 1 Γ ( 2 β + 2 ) τ 3 β + 1 Γ ( 3 β + 2 ) + τ 4 β + 1 Γ ( 4 β + 2 ) + ,

For β = 2, Eq. (30) reduces to

(31) ζ ( η , τ ) = 1 + η 5 120 τ τ 3 3 ! + τ 5 5 ! τ 7 7 ! + τ 9 9 ! ,

which converges to the exact solution as

(32) ζ ( η , τ ) = 1 + η 5 120 sin ( τ ) .

### 4.2 Problem 2

Consider the (2 + 1) dimension parabolic equation of the form 

(33) D τ β c ζ ( η , υ , τ ) + 2 1 η 2 + η 4 6 ! ζ η η η η ( η , υ , τ ) + 2 1 υ 2 + υ 4 6 ! ζ υ υ υ υ ( η , υ , τ ) = 0 , 1 < β 2 , τ 0 ,

with initial conditions

(34) ζ ( η , υ , τ ) = 0 , ζ τ ( η , υ , τ ) = 2 + η 6 6 ! + υ 6 6 ! ,

and boundary conditions as

(35) ζ 1 2 , υ , τ = 2 + ( 1 / 2 ) 2 6 ! + υ 6 6 ! sin ( τ ) , ζ 1 2 , υ , τ = 2 + ( 1 ) 6 6 ! + υ 6 6 ! sin ( τ ) , ζ η 1 2 , υ , τ = ( 1 / 2 ) 2 4 ! sin ( τ ) , ζ η 1 2 , υ , τ = 1 24 sin ( τ ) , ζ η η , 1 2 , τ = ( 1 / 2 ) 2 4 ! sin ( τ ) , ζ η η , 1 2 , τ = 1 24 sin ( τ ) ,

Rearranging Eq. (33) as

(36) c D τ β ζ ( η , υ , τ ) = 2 1 η 2 + η 4 6 ! ζ η η η η ( η , υ , τ ) 2 1 υ 2 + υ 4 6 ! ζ υ υ υ υ ( η , υ , τ ) ,

Applying the basic procedure of NTIM and using the initial conditions, we obtain the solution components as

(37) ζ 0 ( η , υ , τ ) = 2 + η 6 6 ! + υ 6 6 ! τ , ζ 1 ( η , υ , τ ) = 2 + η 6 6 ! + υ 6 6 ! τ β + 1 Γ ( β + 2 ) , ζ 2 ( η , υ , τ ) = 2 + η 6 6 ! + υ 6 6 ! τ 2 β + 1 Γ ( 2 β + 2 ) , ζ 3 ( η , υ , τ ) = 2 + η 6 6 ! + υ 6 6 ! τ 3 β + 1 Γ ( 3 β + 2 ) , ζ 4 ( η , υ , τ ) = 2 + η 6 6 ! + υ 6 6 ! τ 4 β + 1 Γ ( 4 β + 2 ) ,

Adding the solution components, we get the approximate solution as

(38) ζ ( η , υ , τ ) = ζ 0 ( η , υ , τ ) + ζ 1 ( η , υ , τ ) + ζ 2 ( η , υ , τ ) + ζ 3 ( η , υ , τ ) + ζ 4 ( η , υ , τ ) + = 2 + η 6 6 ! + υ 6 6 ! τ τ β + 1 Γ ( β + 2 ) + τ 2 β + 1 Γ ( 2 β + 2 ) τ 3 β + 1 Γ ( 3 β + 2 ) + τ 4 β + 1 Γ ( 4 β + 2 ) + ,

For β = 2, Eq. (33) reduces to

(39) ζ ( η , υ , τ ) = 2 + η 6 6 ! + υ 6 6 ! τ τ 3 3 ! + τ 5 5 ! τ 7 7 ! + τ 9 9 ! ,

which converges to the exact solution given as

(40) ζ ( η , υ , τ ) = 2 + η 6 6 ! + υ 6 6 ! sin ( τ ) .

### 4.3 Problem 3

Consider the (3 + 1) dimension parabolic equation of the form 

(41) D τ β c ζ ( η , υ , ψ , τ ) + υ + ψ 2 cos ( η ) 1 ζ η η η η ( η , υ , ψ , τ ) + η + ψ 2 cos ( υ ) 1 ζ υ υ υ υ ( η , υ , ψ , τ ) + η + υ 2 cos ( ψ ) 1 ζ ψ ψ ψ ψ ( η , υ , ψ , τ ) = 0 , 1 < β 2 , τ 0 ,

with initial conditions

(42) ζ ( η , υ , ψ , 0 ) = η + υ + ψ ( cos ( η ) + cos ( ψ ) + cos ( ψ ) ) , ζ τ ( η , υ , ψ , 0 ) = ( cos ( η ) + cos ( ψ ) + cos ( ψ ) ) ( η + υ + ψ ) ,

and boundary conditions as

(43) ζ ( 0 , υ , ψ , τ ) = ( 1 + υ + ψ cos ( υ ) cos ( ψ ) ) e τ , ζ π 3 , υ , ψ , τ = 2 π 3 6 + υ + ψ cos ( υ ) cos ( ψ ) e τ , ζ ( η , 0 , ψ , τ ) = ( 1 + η + ψ cos ( η ) cos ( ψ ) ) e τ , ζ η , π 3 , ψ , τ = 2 π 3 6 + η + ψ cos ( η ) cos ( ψ ) e τ , ζ ( η , υ , 0 , τ ) = ( 1 + η + υ cos ( η ) cos ( υ ) ) e τ ,

Rearranging Eq. (41) as

(44) D τ β c ζ ( η , υ , ψ , τ ) = υ + ψ 2 cos ( η ) 1 ζ η η η η ( η , υ , ψ , τ ) η + ψ 2 cos ( υ ) 1 ζ υ υ υ υ ( η , υ , ψ , τ ) η + υ 2 cos ( ψ ) 1 ζ ψ ψ ψ ψ ( η , υ , ψ , τ ) ,

Applying the basic procedure of NTIM and using the initial conditions, we obtain the solution components as

(45) ζ 0 ( η , υ , ψ , τ ) = [ η + υ + ψ ( cos ( η ) + cos ( ψ ) + cos ( ψ ) ) ] ( 1 τ ) , ζ 1 ( η , υ , ψ , τ ) = [ η + υ + ψ ( cos ( η ) + cos ( ψ ) + cos ( ψ ) ) ] τ β Γ ( β + 1 ) τ β + 1 Γ ( β + 2 ) , ζ 2 ( η , υ , ψ , τ ) = [ η + υ + ψ ( cos ( η ) + cos ( ψ ) + cos ( ψ ) ) ] τ 2 β Γ ( 2 β + 1 ) τ 2 β + 1 Γ ( 2 β + 2 ) , ζ 3 ( η , υ , ψ , τ ) = [ η + υ + ψ ( cos ( η ) + cos ( ψ ) + cos ( ψ ) ) ] τ 3 β Γ ( 3 β + 1 ) τ 3 β + 1 Γ ( 3 β + 2 ) ,

Adding the solution components, we get the approximate solution as

(46) ζ ( η , υ , ψ , τ ) = ζ 0 ( η , υ , ψ , τ ) + ζ 1 ( η , υ , ψ , τ ) + ζ 2 ( η , υ , ψ , τ ) + ζ 3 ( η , υ , ψ , τ ) + = [ η + υ + ψ ( cos ( η ) + cos ( ψ ) + cos ( ψ ) ) ] 1 τ + τ β Γ ( β + 1 ) τ β + 1 Γ ( β + 2 ) + τ 2 β Γ ( 2 β + 1 ) τ 2 β + 1 Γ ( 2 β + 1 ) + τ 3 β Γ ( 3 β + 1 ) τ 3 β + 1 Γ ( 3 β + 1 ) + ,

For β = 2, Eq. (41) reduces to

(47) ζ ( η , υ , ψ , τ ) = [ η + υ + ψ ( cos ( η ) + cos ( ψ ) + cos ( ψ ) ) ] 1 τ + τ 2 2 ! τ 3 3 ! + τ 4 4 ! τ 5 5 ! + τ 6 6 ! τ 7 7 ! + ,

which converges to the exact solution given as

(48) ζ ( η , υ , ψ , τ ) = [ η + υ + ψ ( cos ( η ) + cos ( ψ ) + cos ( ψ ) ) ] e τ .

## 5 Numerical results and discussions

We have obtained the analytical approximate solution of the fourth-order parabolic time fractional PDEs by applying the natural transform in combination with the NIM. It is noted that the solution pattern converges to the exact solution after a few iterations which shows the reliability of our proposed method. All the solutions are approximated up to fourth-order for problems 1–3. Figures 1 and 2 show the approximate solution and exact solution, respectively, for problem 1. Figure 3 shows the absolute error for problem 1. Figures 4 and 5 show the absolute error and relative error by 2D plots. Figures 6 and 7 show the comparison of the fractional value of β as it converges to the exact solution when the value of β approaches to 2. Similarly, Figures 8 and 9 represent the approximate and exact solution in 3D plots for problem 2. The absolute error is shown by 3D plot in Figure 10 by keeping one parameter ղ constant. Figures 11 and 12 are the 2D graphs of the absolute and relative errors, respectively. The approximate solution is compared with the exact solution by giving different values to β in Figures 13 and 14. The approximate solution and exact solution for problem 3 have been shown through 3D plots by Figures 15 and 16, respectively. The absolute error is shown for the said problem in Figure 17. The absolute and relative errors are shown in Figures 18 and 19, respectively, with the help of 2D plots. The different fractional values are compared with exact solution by giving different values to β (Figures 20 and 21). The numerical values for different values of β are compared in Tables 1 and 2 for problem 1, Tables 3 and 4 for problem 2, and Tables 5 and 6 for problem 3. As the value of β approaches 2, which changes the FDE to a classical PDE, the approximate solution converges to the exact solution, as seen in the figures and tables. The parameters have been interpreted by keeping some parameters constant and others to expand through 2D and 3D plots. All the figures show that the method is convergent and has excellent degree of accuracy. Figure 1

Approximate solution of problem 1 at β = 2 . Figure 2

Exact solution of problem 1 at β = 2 . Figure 3

Absolute error of problem 1 for β = 2 . Figure 4

Absolute error of problem 1 for τ = 0.1 , β = 2 . Figure 5

Relative error of problem 1 for τ = 0.1 , β = 2 . Figure 6

Approximate solution of problem 1 for different values of β at η = 1 . Figure 7

Approximate solution of problem 1 for different values of β at τ = 1 . Figure 8

Approximate solution of problem 2 at υ = 1 , β = 2 . Figure 9

Exact solution of problem 2 at υ = 1 , β = 1 . Figure 10

Absolute error of problem 2 at υ = 1 , β = 2 . Figure 11

Absolute error of problem 2 for τ = 0.1 , υ = 1 , β = 2 . Figure 12

Relative error of problem 2 for τ = 0.1 , υ = 1 , β = 2 . Figure 13

Approximate solution of problem 2 for different values of β at η = υ = 1 . Figure 14

Approximate solution of problem 2 for different values of β at τ = υ = 1 . Figure 15

Approximate solution of problem 3 for υ = 1 , ψ = 1 , β = 2 . Figure 16

Exact solution of problem 3 for υ = 1 , ψ = 1 , β = 2 . Figure 17

Absolute error of problem 3 for υ = 1 , ψ = 1 , β = 2 . Figure 18

Absolute error of problem 3 for τ = 0.1 , υ = 1 , ψ = 1 , β = 2 . Figure 19

Relative error of problem 3 for τ = 0.1 , υ = 1 , ψ = 1 , β = 2 . Figure 20

Approximate solution of problem 3 for different values of β at η = υ = ψ = 1 . Figure 21

Approximate solution of problem 3 for different values of β at τ = υ = ψ = 1 .

Table 1

Numerical comparison of fourth-order approximate solution for problem 1 at η = 1

τ β = 1.5 β = 1.7 β = 1.9 β = 2.0 Exact Abs. error
0.2 0.196306 0.198551 0.199884 0.200325 0.200325 2.77556 × 10−17
0.4 0.373684 0.383362 0.390131 0.392663 0.392663 1.16573 × 10−15
0.6 0.525633 0.546467 0.562749 0.569348 0.569348 2.11164 × 10−13
0.8 0.649208 0.68249 0.711002 0.723334 0.723334 8.87512 × 10−12
1.0 0.743628 0.788051 0.829477 0.848483 0.848483 1.61161 × 10−10
1.2 0.809641 0.861515 0.914196 0.939806 0.939806 1.72071 × 10−9
1.4 0.849104 0.902788 0.962692 0.993662 0.993662 1.27334 × 10−8
1.6 0.864641 0.913114 0.974042 1.0079 1.0079 7.20455 × 10−8
1.8 0.859358 0.894854 0.948843 0.981963 0.981963 3.32043 × 10−7
2.0 0.836591 0.851263 0.889127 0.916874 0.916875 1.30162 × 10−6
Table 2

Numerical comparison of fourth-order approximate solution for problem 1 at τ = 1

η β = 1.5 β = 1.7 β = 1.9 β = 2.0 Exact Abs. error
0.2 0.737484 0.78154 0.822624 0.841473 0.841473 1.59829 × 10−10
0.4 0.737545 0.781604 0.822692 0.841543 0.841543 1.59842 × 10−10
0.6 0.73796 0.782044 0.823155 0.842016 0.842016 1.59932 × 10−10
0.8 0.739496 0.783672 0.824868 0.843769 0.843769 1.60265 × 10−10
1.0 0.743628 0.788051 0.829477 0.848483 0.848483 1.61160 × 10−10
1.2 0.752774 0.797744 0.83968 0.85892 0.85892 1.63143 × 10−10
1.4 0.770535 0.816565 0.859491 0.879185 0.879185 1.66992 × 10−10
1.6 0.801924 0.84983 0.894504 0.915000 0.915000 1.73795 × 10−10
1.8 0.853609 0.904602 0.952155 0.973972 0.973972 1.84996 × 10−10
2.0 0.934144 0.989948 1.04199 1.06586 1.06586 2.02449 × 10−10
Table 3

Numerical comparison of fourth-order approximate solution for problem 2 at η = 1

τ β = 1.5 β = 1.7 β = 1.9 β = 2.0 Exact Abs. error
0.2 0.389908 0.394367 0.397014 0.397891 0.397891 9.99201 × 10−16
0.4 0.74222 0.761443 0.774889 0.779918 0.779918 2.10221 × 10−12
0.6 1.04403 1.08541 1.11775 1.13085 1.13085 1.81609 × 10−10
0.8 1.28948 1.35558 1.41221 1.4367 1.4367 4.29227 × 10−9
1.0 1.47703 1.56525 1.64753 1.68528 1.68528 4.98537 × 10−8
1.2 1.60821 1.71117 1.8158 1.86667 1.86667 3.69378 × 10−7
1.4 1.68681 1.79318 1.91213 1.97364 1.97364 2.00653 × 10−6
1.6 1.71828 1.8138 1.93469 2.00193 2.00192 8.68357 × 10−6
1.8 1.70936 1.77785 1.8847 1.95043 1.9504 3.15864 × 10−5
2.0 1.66774 1.69208 1.76625 1.82122 1.82112 1.00171 × 10−4
Table 4

Numerical comparison of fourth-order approximate solution for problem 2 at τ = 1

η β = 1.5 β = 1.7 β = 1.9 β = 2.0 Exact Abs. error
0.2 1.47601 1.56416 1.64639 1.68411 1.68411 4.98191 × 10−8
0.4 1.47601 1.56417 1.64639 1.68412 1.68412 4.98193 × 10−8
0.6 1.47605 1.56421 1.64644 1.68417 1.68417 4.98207 × 10−8
0.8 1.47627 1.56445 1.64669 1.68442 1.68442 4.98282 × 10−8
1.0 1.47703 1.56525 1.64753 1.68528 1.68528 4.98537 × 10−8
1.2 1.47906 1.5674 1.6498 1.6876 1.6876 4.99224 × 10−8
1.4 1.48372 1.57234 1.65499 1.69291 1.69291 5.00794 × 10−8
1.6 1.49319 1.58237 1.66555 1.70372 1.70372 5.03992 × 10−8
1.8 1.51084 1.60108 1.68525 1.72386 1.72386 5.09950 × 10−8
2.0 1.54156 1.63363 1.71951 1.75891 1.75891 5.20318 × 10−8
Table 5

Numerical comparison of fourth-order approximate solution for problem 3 at η = 1

τ β = 1.5 β = 1.7 β = 1.9 β = 2.0 Exact Abs. error
0.2 1.19041 1.15741 1.13645 1.12911 1.12911 3.86358 × 10−14
0.4 1.06156 0.993212 0.943736 0.924434 0.924434 3.84487 × 10−11
0.6 0.962711 0.864368 0.788285 0.756862 0.756862 2.17861 × 10−9
0.8 0.883437 0.760652 0.661903 0.619666 0.619666 3.80252 × 10−8
1.0 0.818041 0.675718 0.558569 0.50734 0.50734 3.48164 × 10−7
1.2 0.762877 0.605227 0.473684 0.415373 0.415375 2.11990 × 10−6
1.4 0.715275 0.546021 0.40366 0.34007 0.34008 9.74103 × 10−6
1.6 0.672942 0.495657 0.345645 0.278398 0.278434 3.64289 × 10−5
1.8 0.633511 0.452083 0.29732 0.227846 0.227963 1.16411 × 10−4
2.0 0.594103 0.413349 0.25673 0.186311 0.18664 3.28612 × 10−4
Table 6

Numerical comparison of fourth-order approximate solution for problem 3 at τ = 1

η β = 1.5 β = 1.7 β = 1.9 β = 2.0 Exact Abs. error
0.2 0.0826461 0.0682673 0.0564318 0.0512562 0.0512562 1.75874 × 10−8
0.4 0.236281 0.195173 0.161336 0.146539 0.146539 5.02815 × 10−8
0.6 0.411698 0.34007 0.281112 0.25533 0.25533 8.76108 × 10−8
0.8 0.606632 0.501089 0.414216 0.376226 0.376226 1.29093 × 10−7
1.0 0.818041 0.675718 0.558569 0.50734 0.50734 1.74082 × 10−7
1.2 1.04223 0.8609 0.711646 0.646378 0.646378 2.21790 × 10−7
1.4 1.27498 1.05316 0.870575 0.79073 0.79073 2.71321 × 10−7
1.6 1.51176 1.24874 1.03225 0.937575 0.937576 3.21708 × 10−7
1.8 1.74784 1.44375 1.19345 1.08399 1.08399 3.71947 × 10−7
2.0 1.97855 1.63432 1.35098 1.22708 1.22708 4.21044 × 10−7

## 6 Conclusion

The proposed method is tested upon the time fractional parabolic equations of fourth order. The method shows its convergence as the solution pattern after a few iterations converges to the exact solution. The algorithm of the proposed method is easy to apply without any discretization. The approximate solutions are found to be in excellent agreement with the exact solution. The numerical values of approximate solution for different fractional values have been compared in tables, and the results have been shown in 3D and 2D graphs. We found that for β = 2, the suggested method’s approximate solution converges to the precise solution of the problems, ensuring NTIM’s dependability and correctness. The results show that the NTIM successfully delivers accurate, faster converging solutions while using fewer computer resources than other approaches in the literature.

1. Funding information: The authors state no funding involved.

2. Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.

3. Conflict of interest: The authors state no conflict of interest.

## References

 Alshammari M, Iqbal N, Ntwiga DB. A comparative study of fractional-order diffusion model within Atangana-Baleanu-Caputo operator. J Funct Spaces. 2022 Apr 30;2022(9):1–12.10.1155/2022/9226707Search in Google Scholar

 Borhanifar A, Ragusa MA, Valizadehaz S. High-order numerical method for two-dimensional Riesz space fractional advection-dispersion equation. arXiv Prepr arXiv:200604111. 2020 Jun 710.3934/dcdsb.2020355Search in Google Scholar

 Özdemir ME. New Refinements of Hadamard Integral inequlaity via k-Fractional Integrals for p-convex function. Turkish J Sci. 2021;6(1):1–5.Search in Google Scholar

 Mainardi F, Raberto M, Gorenflo R, Scalas E. Fractional calculus and continuous-time finance II: The waiting-time distribution. Phys A: Stat Mech its Appl. 2000 Dec 1;287(3–4):468–81.10.1016/S0378-4371(00)00386-1Search in Google Scholar

 Sabatier JA, Agrawal OP, Machado JT. Advances in fractional calculus. Dordrecht: Springer; 2007.10.1007/978-1-4020-6042-7Search in Google Scholar

 Hilfer R, editor. Applications of fractional calculus in physics. Singapore: World scientific; 2000 Mar 2.10.1142/3779Search in Google Scholar

 Phuong ND, Tuan NA, Kumar D, Tuan NH. Initial value problem for fractional Volterra integrodifferential pseudo-parabolic equations. Math Model Nat Phenom. 2021;16:27.10.1051/mmnp/2021015Search in Google Scholar

 Momani S, Odibat Z. Numerical comparison of methods for solving linear differential equations of fractional order. Chaos Solitons & Fractals. 2007 Mar 1;31(5):1248–55.10.1016/j.chaos.2005.10.068Search in Google Scholar

 Zada L, Nawaz R. Solution of time-fractional order RLW equation using optimal homotopy asymptotic method. In AIP Conference Proceedings. Vol. 2116. 1. AIP Publishing LLC; 2019 Jul 24. p. 300005.10.1063/1.5114305Search in Google Scholar

 Odibat Z, Momani S. Numerical methods for nonlinear partial differential equations of fractional order. Appl Math Model. 2008 Jan 1;32(1):28–39.10.1016/j.apm.2006.10.025Search in Google Scholar

 Nawaz R, Zada L, Khattak A, Jibran M, Khan A. Optimum solutions of fractional order Zakharov–Kuznetsov equations. Complexity. 2019 Dec 10;2019:1741958.10.1155/2019/1741958Search in Google Scholar

 Singh J, Gupta A, Baleanu D. On the analysis of an analytical approach for fractional Caudrey-Dodd-Gibbon equations. Alex Eng J. 2022 Jul 1;61(7):5073–82.10.1016/j.aej.2021.09.053Search in Google Scholar

 Dubey VP, Kumar D, Alshehri HM, Dubey S, Singh J. Computational analysis of local fractional LWR model occurring in a fractal vehicular traffic flow. Fractal Fract. 2022;6(8):426.10.3390/fractalfract6080426Search in Google Scholar

 Yadav S, Kumar D, Nisar KS. A reliable numerical method for solving fractional reaction-diffusion equations. J King Saud Univ Sci. 2021;33(2):101320.10.1016/j.jksus.2020.101320Search in Google Scholar

 Singh J. Analysis of fractional blood alcohol model with composite fractional derivative. Chaos Solitons & Fractals. 2020 Nov 1;140:110127.10.1016/j.chaos.2020.110127Search in Google Scholar

 Nawaz R, Ali N, Zada L, Nisar KS, Alharthi MR, Jamshed W. Extension of natural transform method with Daftardar-Jafari polynomials for fractional order differential equations. Alex Eng J. 2021 Jun 1;60(3):3205–17.10.1016/j.aej.2021.01.051Search in Google Scholar

 Ali N, Nawaz R, Zada L, Mouldi A, Bouzgarrou SM, Sene N. Analytical Approximate Solution of the Fractional Order Biological Population Model by Using Natural Transform. J Nanomaterials. 2022 Mar 19;2022:6703086.10.1155/2022/6703086Search in Google Scholar

 Nawaz R, Ali N, Zada L, Shah Z, Tassaddiq A, Alreshidi NA. Comparative analysis of natural transform decomposition method and new iterative method for fractional foam drainage problem and fractional order modified regularized long-wave equation. Fractals. 2020 Nov 23;28(7):2050124.10.1142/S0218348X20501248Search in Google Scholar

 Ali N, Nawaz R, Zada L, Nisar KS, Ali Z, Jamshed W, et al. Numerical investigation of generalized perturbed Zakharov–Kuznetsov equation of fractional order in dusty plasma. Waves Random Complex Media. 2022 Feb 15;1–20.10.1080/17455030.2021.2004332Search in Google Scholar

 Ali N, Yassen MF, Asiri SA, Nawaz R, Zada L, Alam MM, et al. New iterative method for solving a coupled system of fractional-order Drinfeld–Sokolov–Wilson (FDSW) and Fractional Shallow Water (FSW) equations. J Nanomaterials. 2022 Apr 7;2022:8370107.10.1155/2022/8370107Search in Google Scholar

 Bhalekar S, Daftardar-Gejji V. New iterative method: application to partial differential equations. Appl Math Comput. 2008 Sep 15;203(2):778–83.10.1016/j.amc.2008.05.071Search in Google Scholar

 Bhalekar S, Daftardar-Gejji V. Convergence of the new iterative method. Int J Differ Equ. 2011 Jan 1;2011:989065.10.1155/2011/989065Search in Google Scholar

 Naeem M, Azhar OF, Zidan AM, Nonlaopon K, Shah R. Numerical analysis of fractional-order parabolic equations via Elzaki transform. J Funct Spaces. 2021 Sep 1;2021:3484482.10.1155/2021/3484482Search in Google Scholar

 Yadeta DM, Gizaw AK, Mussa YO. Approximate analytical solution of one-dimensional Beam equations by using time-fractional reduced differential transform method. J Appl Math. 2020 Dec 22;2020:7627385.10.1155/2020/7627385Search in Google Scholar

 Khalid N, Abbas M, Iqbal MK. Non-polynomial quintic spline for solving fourth-order fractional boundary value problems involving product terms. Appl Math Computation. 2019 May 15;349:393–407.10.1016/j.amc.2018.12.066Search in Google Scholar

 Tariq H, Akram G. Quintic spline technique for time fractional fourth‐order partial differential equation. Numer Methods Partial Differ Equ. 2017 Mar;33(2):445–66.10.1002/num.22088Search in Google Scholar

 Hamaidi M, Naji A, Taik A. Solving parabolic and hyperbolic equations with variable coefficients using space-time localized radial basis function collocation method. Model Simul Eng. 2021 Feb 8;2021:6688806.10.1155/2021/6688806Search in Google Scholar

 Almuqrin MA, Goswami P, Sharma S, Khan I, Dubey RS, Khan A. Fractional model of Ebola virus in population of bats in frame of Atangana-Baleanu fractional derivative. Results Phys. 2021 Jul 1;26:104295.10.1016/j.rinp.2021.104295Search in Google Scholar

 Gorman DJ. Free vibration analysis of beams and shafts (Book). Research supported by the National Research Council of Canada. Vol. 395. New York: Wiley-Interscience; 1975. p. 1975.Search in Google Scholar

 Wazwaz AM. On the solution of the fourth order parabolic equation by the decomposition method. Int J Comput Math. 1995 Jan 1;57(3–4):213–7.10.1080/00207169508804424Search in Google Scholar

 Khaliq AQ, Twizell EH. A family of second order methods for variable coefficient fourth order parabolic partial differential equations. Int J Comput Math. 1987 Jan 1;23(1):63–76.10.1080/00207168708803608Search in Google Scholar

 Dehghan M, Manafian J. The solution of the variable coefficients fourth-order parabolic partial differential equations by the homotopy perturbation method. Z für Naturforsch A. 2009 Aug 1;64(7–8):420–30.10.1515/zna-2009-7-803Search in Google Scholar

 Andrade C, McKee S. High accuracy ADI methods for fourth order parabolic equations with variable coefficients. J Comput Appl Math. 1977 Mar 1;3(1):11–4.10.1016/0771-050X(77)90019-5Search in Google Scholar