Numerical Simulation for Fractional-Order Bloch Equation Arising in Nuclear Magnetic Resonance by Using the Jacobi Polynomials

: In the present paper, we numerically simulate fractional-order model of the Bloch equation by using the Jacobi polynomials. It arises in chemistry, physics and nuclear magnetic resonance (NMR). It also arises in magnetic resonance imaging (MRI) and electron spin resonance (ESR). It is used for purity determination, provided that the molecular weight and structure of the compound is known. It can also be used for structural determination. By the study of NMR, chemists can determine the structure of many compounds. The obtained numerical results are compared and simulated with the known solutions. Accuracy of the proposed method is shown by providing tables for absolute errors and root mean square errors. Di ﬀ erent orders of the time-fractional derivatives results are illustrated by using ﬁgures.


Introduction
The Bloch equation is a system of differential equations. It is mainly valuable for studying expensive biological samples like RNA, DNA, proteins and nucleic acids. It has many real-life applications like process control, liquid media, petrochemical plants and process optimization in oil refineries. Surface magnetic resonance is based on the principle of NMR, and the measurements can be used to indirectly estimate the water content of the saturated and unsaturated zones. The standard system of Bloch equations is given as follows: with the initial conditions N x (0) = a 1 , N y (0) = a 2 and N z (0) = a 3 .
Here N x (t), N y (t) and N z (t) denote the system magnetisation in x, y and z components, respectively; w 0 is the resonant frequency given by the Larmor relationship w 0 = γM 0 , where M 0 is the static magnetic field in z−component, N 0 is the equilibrium magnetisation, T 1 and T 2 are the spin-lattice and spin-spin relaxation time, respectively, and a 1 , a 2 and a 3 are real constants. The analytical solution of Equation (1) is given by Real-life applications of fractional calculus are in subjects like biology [1], viscoelasticity [2][3][4], signal processing [5], control theory [6], fluid dynamics [7]. For more details, the reader should refer to [8]. There are many magnetic resonance systems which are modelled by the fractional-order Bloch equation, and it is well known that the fractional-order derivatives are non-local in nature. Therefore, we will replace the integer-order Bloch equation by the fractional-order Bloch equation with a view to further understand the resulting magnetic resonance systems. Therefore, we replace the integer-order time-derivative by the non-integer-order time-derivative: where 0 < α, β, γ ≤ 1.
In this paper, we are considering that β ∈ (0, 1); therefore, we will take l = 1. The time-fractional derivatives play a key role upsetting the spin dynamics defined by the Bloch equations in Equation (3) (see [9,10]). The magnetic resonance components of the magnetisation are identified in the initial state of the system, and hence, these should be visibly predictable. The physical meaning of the non-integer order Bloch equations can be understood in the basic preparation of the non-integer-order Schrödinger equation.
Bloch equations in NMR can be simulated numerically and analytically (see, for details, [11][12][13][14][15][16]). The time-fractional order Bloch equation having fractional derivative in Caputo sense is solved in [17]. Recently, Kumar et al. [18] solved fractional-order Bloch equation by using homotopy perturbation method (HPM). Use has been made of operational matrix method with Legendre polynomials in [19] and with the Laguerre polynomials in [20] for the solution of this equation. In [21], this equation was solved numerically by using the iterative method. Furthermore, in [22], by using numerical approximation, a special class of this equation, namely the fuzzy time-fractional Bloch equation, was solved. In this paper we propose to solve the fractional-order Bloch equation by using the Jacobi polynomials. Some developments on orthogonal approximations can be found in [23][24][25][26][27][28][29][30]. Some introductory overview and recent development of fractional calculus can be seen in [31]. In this method, we get unknown coefficients for the approximated parameter in the model and, by the use of these coefficients, we obtain approximate solutions of the fractional-order Bloch equation in NMR.

Preliminaries
The Jacobi polynomial of degree i on [0, 1] is given by [28] The orthogonal property of the Jacobi polynomials with respect to the weight function w (a,b) (t) = (1 − t) a t b is given by where δ mn is Kronecker delta function and v a,b n = can be expanded as follows: where (7), for finite dimensional approximation, is written in the following form: where C and q m (t) are (m + 1) × 1 matrices given by C = [c 0 , c 1 , . . . ., c m ] T and q m (t) = [σ 0 , σ 1 , . . . ., σ m ] T . Theorem 1. If q n (t) = [σ 0 , σ 1 , . . . ., σ n ] T denotes the shifted Jacobi vector and if v > 0, then I v σ i (t) = I (v) q n (t), where I (v) = (u(i, j)), is the (n + 1) × (n + 1) operational matrix of fractional integral of order v, and its (i, j)th entry is given by Proof. Please see [28].

Construction of Algorithm
In this section, we construct an algorithm to get the approximate solution of the Bloch equation. Using this algorithm, we can then obtain magnetisation in each direction.

of 18
From Equations (10) and (11), we can write Using Equations (10), (12), (13) and (14) in Equation (3), we get where I (α) , I (β) and I (γ) are the operational matrices of non-integer-order integration of order α, β and γ, respectively. Here I is an identity matrix. The simpler form for Equations (15)-(17) is given as follows: where The matrices W 1 , W 2 , W 3 , W 4 , W 5 , G 1 , G 2 and G 3 are given in terms of known values, so these matrices are known matrices.

of 18
Using Equations (29)- (31) in Equations (12)- (14), respectively, we get the system magnetisation N x (t), N y (t) and N z (t) for Bloch equations in NMR. 1], so the Taylor polynomial of d α N x dt α at t = 0 is given as follows:

Convergence Analysis
The upper bound of the error of the Taylor polynomial is given by where Taking m → ∞ in Equation (35), we get

Numerical Results and Discussion
In this section, we will numerically simulate our results with known results. For each numerical simulation, we will consider i. c. N x (0) = 0, N y (0) = 100 and N z (0) = 0. In Figures 1 and 2, we have presented 3D and 2D plots of the numerical solutions of the Bloch equation for integer order, respectively. These figures show the dynamics of N x , N y and N z for integer-order relaxation. In Figure 1, the entire trajectory of magnetisation is shown in 3D for integer order starting at i. c. N x (0), N y (0), N z (0) and returning to N 0 . From Figure 2, it is clear that the magnetisation N x in x−direction increases with time, and the magnetisation N y in y−direction decreases with time.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 18  In Figures 3 and 4, we have presented 3D and 2D plots of the numerical solutions of the Bloch equation for fractional order ( = = = 0.9), respectively. These figures show the dynamics of , and for fractional-order relaxation. In Figure 3, the entire trajectory of magnetisation is shown in 3D for fractional order ( = = = 0.9) starting at i. c.    In Figures 3 and 4, we have presented 3D and 2D plots of the numerical solutions of the Bloch equation for fractional order ( = = = 0.9), respectively. These figures show the dynamics of , and for fractional-order relaxation. In Figure 3, the entire trajectory of magnetisation is shown in 3D for fractional order ( = = = 0.9) starting at i. c.
(0), (0), (0) and returning to . From Figure 4, it is clear that the magnetisation in −direction increases with time, and the magnetisation in −direction decreases with time.
In Figures 3 and 4, we have presented 3D and 2D plots of the numerical solutions of the Bloch equation for fractional order (α = β = γ = 0.9), respectively. These figures show the dynamics of N x , N y and N z for fractional-order relaxation. In Figure 3, the entire trajectory of magnetisation is shown in 3D for fractional order (α = β = γ = 0.9) starting at i. c. N x (0), N y (0), N z (0) and returning to N 0 . From Figure 4, it is clear that the magnetisation N x in x−direction increases with time, and the magnetisation N y in y−direction decreases with time.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 18  In Figures 5 and 6, we have presented the 3D and 2D plots of the numerical solutions of the Bloch equation for fractional order ( = = = 0.8), respectively. These figures show the dynamics of , and for the fractional-order relaxation. In Figure 5, the entire trajectory of   In Figures 5 and 6, we have presented the 3D and 2D plots of the numerical solutions of the Bloch equation for fractional order ( = = = 0.8), respectively. These figures show the dynamics of , and for the fractional-order relaxation. In Figure 5, the entire trajectory of In Figures 5 and 6, we have presented the 3D and 2D plots of the numerical solutions of the Bloch equation for fractional order (α = β = γ = 0.8), respectively. These figures show the dynamics of N x , N y and N z for the fractional-order relaxation. In Figure 5, the entire trajectory of magnetisation is shown in 3D for fractional order (α = β = γ = 0.8) starting at i. c. N x (0), N y (0), N z (0) and returning to N 0 . From Figure 6, it is clear that the magnetisation N x in x−direction increases with time, and the magnetisation N y in y−direction decreases with time. magnetisation is shown in 3D for fractional order ( = = = 0.8) starting at i. c.
(0), (0), (0) and returning to . From Figure 6, it is clear that the magnetisation in −direction increases with time, and the magnetisation in −direction decreases with time.   −direction increases with time, and the magnetisation in −direction decreases with time.   In Figures 7 and 8, we have presented the 3D and 2D plots of the numerical solutions of the Bloch equation for fractional order (α = β = γ = 0.7), respectively. These figures show the dynamics of N x , N y and N z for the fractional-order relaxation. In Figure 7, the entire trajectory of magnetisation is shown in 3D for fractional order (α = β = γ = 0.7) with the starting initially N x (0), N y (0), N z (0) and returning to N 0 . From Figure 8, it is clear that the magnetisation N x in x−direction increases with time, and the magnetisation N y in y−direction decreases with time.  In Figures 9 and 10, we have presented the 3D and 2D plots of the numerical solutions of the Bloch equation for fractional order ( = 1.0, = 0.9, = 0.8), respectively. These figures show the dynamics of , and for the fractional-order relaxation. In Figure 9, the entire trajectory of magnetisation is shown in 3D for fractional order ( = 1.0, = 0.9, = 0.8) starting at i. c.   In Figures 9 and 10, we have presented the 3D and 2D plots of the numerical solutions of the Bloch equation for fractional order ( = 1.0, = 0.9, = 0.8), respectively. These figures show the dynamics of , and for the fractional-order relaxation. In Figure 9, the entire trajectory of magnetisation is shown in 3D for fractional order ( = 1.0, = 0.9, = 0.8) starting at i. c.
In Figures 9 and 10, we have presented the 3D and 2D plots of the numerical solutions of the Bloch equation for fractional order (α = 1.0, β = 0.9, γ = 0.8), respectively. These figures show the dynamics of N x , N y and N z for the fractional-order relaxation. In Figure 9, the entire trajectory of magnetisation is shown in 3D for fractional order (α = 1.0, β = 0.9, γ = 0.8) starting at i. c. N x (0), N y (0), N z (0) and returning to N 0 . From Figure 10, it is clear that the magnetisation N x in x−direction increases with time, and the magnetisation N y in y−direction decreases with time.
(0), (0), (0) and returning to . From Figure 10, it is clear that the magnetisation in −direction increases with time, and the magnetisation in −direction decreases with time.      In Figures 11 and 12 we have presented 3D and 2D plots of numerical solutions of Bloch equation for fractional order (α = 0.9, β = 0.9, γ = 1.0), respectively. These figures show the dynamic of N x , N y and N z for fractional order relaxation. In Figure 11, the entire trajectory of magnetisation is shown in 3D for fractional order (α = 0.9, β = 0.9, γ = 1.0) starting at i.c. N x (0), N y (0), N z (0) and returning to N 0 . From Figure 12, it is clear that the magnetisation N x in x−direction increases with time, and the magnetisation N y in y−direction decreases with time.
with time, and the magnetisation in −direction decreases with time.    In Figures 13 and 14, we have presented the 3D and 2D plots of the numerical solutions of the Bloch equation for fractional order (α = 1.0, β = 1.0, γ = 0.9), respectively. These figures show the dynamics of N x , N y and N z for the fractional-order relaxation. In Figure 13, the entire trajectory of magnetisation is shown in 3D for fractional order (α = 1.0, β = 1.0, γ = 0.9) starting at the initial level N x (0), N y (0), N z (0) and returning to N 0 . From Figure 14, it is clear that the magnetisation N x in x−direction increases with time, and the magnetisation N y in y−direction decreases with time. magnetisation is shown in 3D for fractional order ( = 1.0, = 1.0, = 0.9)starting at the initial level (0), (0), (0) and returning to . From Figure 14, it is clear that the magnetisation in −direction increases with time, and the magnetisation in −direction decreases with time.   in −direction increases with time, and the magnetisation in −direction decreases with time.
In Figures 15 and 16, we have shown the numerical simulation of analytical and numerical solutions of the Bloch equation for N x (t) and N y (t) for integer order, respectively. From Figure 16, it is clear that the solution has periodic behaviour at low frequency. This solution varies periodically for N x (t) and N y (t) at low frequency.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 18 In Figures 15 and 16, we have shown the numerical simulation of analytical and numerical solutions of the Bloch equation for ( ) and ( ) for integer order, respectively. From Figure 16, it is clear that the solution has periodic behaviour at low frequency. This solution varies periodically for ( ) and ( ) at low frequency.  In Figures 17-19, we have shown the absolute errors for ( ), ( ) and ( ), respectively, at different values of = 3, 6 and 9. In Figures 17-19, the absolute errors are denoted by , and for = 3, 6 and 9, respectively. In all these figures, and are multiplied by 10 and 10 , respectively.
In Figures 15 and 16, we have shown the numerical simulation of analytical and numerical solutions of the Bloch equation for ( ) and ( ) for integer order, respectively. From Figure 16, it is clear that the solution has periodic behaviour at low frequency. This solution varies periodically for ( ) and ( ) at low frequency.  In Figures 17-19, we have shown the absolute errors for ( ), ( ) and ( ), respectively, at different values of = 3, 6 and 9. In Figures 17-19, the absolute errors are denoted by , and for = 3, 6 and 9, respectively. In all these figures, and are multiplied by 10 and 10 , respectively.
In Figures 17-19, we have shown the absolute errors for N x (t), N y (t) and N z (t), respectively, at different values of m = 3, 6 and 9. In Figures 17-19, the absolute errors are denoted by E 1 , E 2 and E 3 for m = 3, 6 and 9, respectively. In all these figures, E 2 and E 3 are multiplied by 10 4 and 10 5 , respectively. Appl. Sci. 2020, 10, x FOR PEER REVIEW 14 of 18       . Errors for N z (t) at m = 3, 6 and 9, with parameters: From these figures, we can see that the errors decrease with the increase of m. In Figures 20-22, we have shown the behaviour of the solutions of N x (t), N y (t) and N z (t) at different values of α, β and γ, respectively. In Figures 20-22, exact solution means the analytical solution for integer order (α = β = γ = 1) Bloch equation as given by Equation (2).     From these figures it is clear that the solution varies consistently from non-integer order to integer order.
We have calculated these errors for integer order by taking the exact solution as given by Equation (2). ( ) = = Figure 21. Behaviour of the approximate solution of N y (t) at β = 0.6, 0.7 , 0.8, 0.9 and 1, with parameters:   From these figures it is clear that the solution varies consistently from non-integer order to integer order.
We have calculated these errors for integer order by taking the exact solution as given by Equation (2).  From these figures it is clear that the solution varies consistently from non-integer order to integer order.
In Table 1, we have listed the maximum absolute errors (l ∞ ) and the root-mean-square errors (l 2 ) for two different values of m = 4 and 8 at α = β = γ = 1, w 0 = 1, T 1 = 1, T 2 = 20, a = b = 0.8. We have calculated these errors for integer order by taking the exact solution as given by Equation (2). Table 1. Comparison of (l ∞ ) and (l 2 ) errors at a = b = 0.8, m = 5, 8 for integer order solution. From Table 1, it is detected that the errors decrease with the increase of m.

Conclusions and Future Scope
In this paper, we have presented the numerical solution and the simulation for fractional-order and integer-order Bloch equations. Mathematical model for NMR allows us to explore and define magnetisation for spin dynamics at resonance frequency in a static magnetic field. Implementation of our proposed technique is easy in comparison to the existing methods because the operational matrices are easy to construct. The numerical section shows how the solution given by the used technique varies consistently at different values of non-integer-order time-derivatives. Moreover, for integer order, the solution by the used technique is identical to the exact solution for the Bloch equation. The error table shows the accuracy of the proposed method. For future work, we can construct operational matrices for different polynomials in order to attain better exactness.