ABSTRACT

The Method Of Lines (MOL) and Scheifele’s G-functions in the design of algorithms adapted for the numeric integration of parabolic Partial Differential Equations (PDE) in one space dimension are applied. The semi-discrete system of ordinary differential equations in the time direction, obtained by applying the MOL to PDE, is solved with the use of a method of Adapted Series, based on Scheifele’s G-functions. This method integrates exactly unperturbed linear systems of ordinary differential equations, with only one G-function. An implementation of this algorithm is used to approximate the solution of two test problems proposed by various authors. The results obtained by the Dufort-Frankel, Crank-Nicholson and the methods of Adapted Series versus the analytical solution, show the results of mistakes made.

Keywords:

Series Method, Numerical Solutions, Parabolic Initial-Boundary Value Problems, Method of Lines

1. Introduction

Many physical phenomena related to processes in which there are net transfers or transportation of energy, material, etc., are modeled using Partial Differential Equations (PDEs). Parabolic initial-boundary value problems in one dimension, such as $u_t\left(x,t\right)=a{u}_{xx}\left(x,t\right)+f\left(x,t,u\left(x,t\right),u_x\left(x,t\right)\right)$ with $a>0$ , are subjected to conditions with initial values whose solution satisfies the conditions of known boundaries. The positive constant a, is characteristic of each physical location, x is the position variable, $u\left(x,t\right)$ is the field that corresponds to the transport phenomenon in question.

The field $u\left(x,t\right)$ can represent the distribution of voltage in the transmission of a cable or the concentration in distinct processes in the diffusion of liquids. When $u\left(x,t\right)$ represents the temperature distribution of a rod in position x, in instant t, is called the heat equation. The heat equation comes from the study of the flow of electricity in a long cable or a transmission line, in this case, known as the transmission equation. The heat equation is sometimes known as the one-dimensional diffusion equation of a solute in a solvent, as the diffusion of dissolved substances in solution is analagous to heat flow in a solid. Constant a, is in this case called the diffusivity coefficient of the material.

In this article an application of a series method based on Scheifele’s G-functions [1] [2] , adapted to the resolution of a semi-discrete system of ordinary differential equations in the time direction, obtained through the application of the method of lines (MOL) [3] [4] to the given diffusion equation is shown.

Runge-Kutta methods [5] - [10] , are used for resolving this type of differential equations linear systems. These methods don’t need the calculus of partial derivatives of the perturbation function, but they have the problem of a lot of intermediate calculations with a computational cost increase. It is preferred that the numerical methods used to resolve this type of systems verify the property of integrating exactly the homogeneous problem. The methods based on the Scheifele G-functions show this good property [1] [2] [11] [12] [13] . Runge-Kutta methods also have this property but it is necessary to do a transformation in the differential equations. This process is complicated because we need to change the integration order, furthermore the characteristic equation roots have to be known. The G-function series methods are based on a refinement of Taylor expansion [14] [15] and they are applied to problems where the solution has a near sinusoidal behaviour which is not possible to be effectively resolved by Fourier series and Taylor’s expansions.

To illustrate the good performance of the Adapted Series Method in the resolution of the models which describe these types of phenomena, two test problems proposed by various authors are resolved.

Finally, the results obtained with the Dufort-Frankel, Crank-Nicholson and the Adapted Series methods are compared to the analytical solution, presented with the results of the absolute errors committed.

2. Adapted Algorithm of the Numerical Integration of Parabolic PDEs

Let’s consider the equation in partially derivative parabolics (PDEs):

$u_t=a{u}_{xx}+f\left(x,t,u,u_x\right),x\in \left[0,1\right],a>0$ , (1)

subject to the following boundary conditions (BVP):

$\begin{array}{l}{\alpha }_0\left(t\right)u+{\beta }_0\left(t\right)u={\gamma }_0\left(t\right)\text{in}x=0\\ {\alpha }_1\left(t\right)u+{\beta }_1\left(t\right)u={\gamma }_1\left(t\right)\text{in}x=0\end{array}$ , (2)

and initial conditions given as $u\left(x,0\right)$ , with $0\le x\le 1$ .

Using a grid in the range [0, 1], the partial derivatives of the second part of the Equation (1), are substituted for the corresponding approximations by means of the differences. Once these approximations are made, the problem becomes just resolving a system of first order ordinary differential equations in time.

That is to say, the spatial semi-discretisation used in the numerical method of lines (MOL), leads to a system that will adapt the method of series based on G-functions for time integration.

For this purpose, a uniform grid of the spatial interval $x\in \left[0,1\right]$ is generated $\left\{x_0=0,x_1,\cdots ,x_{n-1},x_n=1\right\}$ , where $x_{i+1}=x_i+r$ , $r=1/n$ , with $0\le i\le n-1$ , so that, if we consider the boundary conditions:

$\begin{array}{l}{\alpha }_0\left(t\right)u\left(0,t\right)+{\beta }_0\left(t\right)u_x\left(0,t\right)={\gamma }_0\left(t\right)\\ {\alpha }_1\left(t\right)u\left(1,t\right)+{\beta }_1\left(t\right)u_x\left(1,t\right)={\gamma }_1\left(t\right)\\ \text{with}{\beta }_0\left(t\right)=0\text{and}{\beta }_1\left(t\right)=0,\end{array}$ (3)

we get:

$\begin{array}{l}{u}_0\left(t\right)=u\left(x_0,t\right)=u\left(0,t\right)={\gamma }_0\left(t\right)/{\alpha }_0\left(t\right),\\ u_i\left(t\right)=u\left(x_i,t\right)=u\left(ir,t\right)\text{with}1\le i\le n-1,\\ u_n\left(t\right)=u\left(x_n,t\right)=u\left(1,t\right)={\gamma }_1\left(t\right)/{\alpha }_1\left(t\right).\end{array}$ (4)

Making the spatial semi-discretisation:

$\begin{array}{l}{u}_x\approx \left(u_{i+1}-u_{i-1}\right)/\left(2r\right),\\ u_{xx}\approx \left(u_{i-1}-2u_i+u_{i+1}\right)/r^2,\end{array}$ (5)

the resulting system of ordinary differential equations is:

$\begin{array}{l}{u^{\prime }}_i\left(t\right)=a\left(u_{i-1}-2u_i+u_{i+1}\right)/r^2+f\left(x_i,t,u_i,\left(u_{i+1}-u_{i-1}\right)/\left(2r\right)\right),\\ \text{with}1\le i\le n-1,\end{array}$ (6)

considering the following equations separately:

${u^{\prime }}_0\left(t\right)=u^{\prime }\left(0,t\right)={\left({\gamma }_0\left(t\right)/{\alpha }_0\left(t\right)\right)}^{\prime },$ (7)

${u^{\prime }}_n\left(t\right)=u^{\prime }\left(1,t\right)={\left({\gamma }_1\left(t\right)/{\alpha }_1\left(t\right)\right)}^{\prime }.$ (8)

Like

$\begin{array}{c}{u^{\prime }}_1\left(t\right)=a\left(u_0-2u_1+u_2\right)/r^2+f\left(x_1,t,u_1,\left(u_2-u_0\right)/\left(2r\right)\right)\\ =a\left({\gamma }_0\left(t\right)/{\alpha }_0\left(t\right)-2u_1+u_2\right)/r^2\\ +f\left(r,t,u_1,\left(u_2-{\gamma }_0\left(t\right)/{\alpha }_0\left(t\right)\right)/\left(2r\right)\right)\end{array}$ (9)

and

$\begin{array}{c}{u^{\prime }}_{n-1}\left(t\right)=a\left(u_{n-2}-2u_{n-1}+u_n\right)/r^2+f\left(x_{n-1},t,u_{n-1},\left(u_n-u_{n-2}\right)/\left(2r\right)\right)\\ =a\left(u_{n-2}-2u_{n-1}+{\gamma }_1\left(t\right)/{\alpha }_1\left(t\right)\right)/r^2\\ +f\left(\left(n-1\right)r,t,u_{n-1},\left({\gamma }_1\left(t\right)/{\alpha }_1\left(t\right)-u_{n-2}\right)/\left(2r\right)\right),\end{array}$ (10)

the system of first order differential equations in time (6), can be expressed in matrix form:

$\begin{array}{l}{u}^{\prime }=Au+f\left(u\left(t\right),t\right),\\ u\left(0\right)={\left(\begin{array}{ccc}{u}_1\left(0\right)& \cdots & u_{n-1}\left(0\right)\end{array}\right)}^t\\ ={\left(\begin{array}{ccc}u\left(x_1,0\right)& \cdots & u\left(x_{n-1},0\right)\end{array}\right)}^t,\end{array}$ (11)

where

(12)

To integrate this system, let us assume that the function $f\left(u\left(t\right),t\right)$ defined by the solution $u\left(t\right)$ , to the IVP (11) in an interval [0, T] admits an expansion:

$f\left(u\left(t\right),t\right)={\displaystyle \sum _{k=0}^{\infty }{c}_{k+1}}{t}^k/k!=g\left(t\right)=\left(g_1\left(t\right),\cdots ,g_n\left(t\right)\right),$ (13)

thus, the solution to IVP (11) coincides with the solution of the auxiliary IVP:

$u^{\prime }=Au+{\displaystyle \sum _{k=0}^{\infty }{c}_{k+1}}{t}^k/k!,u\left(0\right)={\left(\begin{array}{ccc}u\left(x_1,0\right)& \cdots & u\left(x_{n-1},0\right)\end{array}\right)}^t=u_0.$ (14)

For the resolution of this IVP (14) the following IVPs are considered:

$\begin{array}{l}{\left({\underset{\_}{u}}_j^k\right)}^{\prime }=A{\underset{\_}{u}}_j^k+{\underset{\_}{e}}^k{t}^j/j!,{\underset{\_}{u}}_j^k\left(0\right)=\underset{\_}{0},\\ j=0,1,2,\cdots \text{,}k=1,2,\cdots ,n-1,\end{array}$ (15)

where ${\underset{\_}{u}}_j^k\in M\left(n,1;ℝ\right)$ , ${\underset{\_}{e}}^k={\left(0,\cdots ,\stackrel{k)}{1},\cdots ,0\right)}^t$ , and by combining their solutions linearly with the coefficients $c_k$ .

In order to simplify the notation we now introduce the sequence $U_j$ where $j=0,1,2,\cdots$ of $\left(n-1\right)\times \left(n-1\right)$ -matrices having the functions ${\underset{\_}{u}}_j^k$ as columns: $U_j=\left({\underset{\_}{u}}_j^1,{\underset{\_}{u}}_j^2,\cdots ,{\underset{\_}{u}}_j^n\right)$ .

The $\left(n-1\right)\times \left(n-1\right)$ -matrices $U_j$ are solutions of the next IVP:

${U^{\prime }}_j=A{U}_j+t^j/j!I,U_j\left(0\right)=\underset{\_}{0},j=0,1,\cdots$ (16)

The solutions to Equation (16) are the Scheifele G-functions [2] and will be denoted as $G_{j+1}\left(t\right)=U_j\left(t\right)$ , $j=0,1,2,\cdots$

$G_0\left(t\right)$ is defined as the solution of IVP to the initial value problem $U^{\prime }=AU$ , $U\left(0\right)=I$ .

The solutions of IVP (11), in terms of G-functions, is given by

$u\left(t\right)=G_0\left(t\right)u_0+{\displaystyle \sum _{k=1}^{\infty }{G}_k\left(t\right)c_k}$ . (17)

Acknowledging that the solution of Equation (11) can be expressed as

$u\left(t\right)={\displaystyle \sum _{k=0}^{\infty }{a}_k{t}^k/k!}$ (18)

and substituting Equation (18) into Equation (11) and identifying coefficients it is obtained:

$a_{k+1}=A{a}_k+c_{k+1}\text{with}k\ge 0\text{and}{a}_0=u_0.$ (19)

What is more, substituting Equation (19) into Equation (17), it is possible to express the solution $u\left(t\right)$ as

$u\left(t\right)=G_0\left(t\right)a_0+{\displaystyle \sum _{k=1}^{\infty }{G}_k\left(t\right)\left(a_k-A{a}_{k-1}\right)},$ (20)

defining a new series of coefficients

$\begin{array}{l}{b}_0=a_0,\\ b_k=a_k-A{a}_{k-1},k\ge 1,\end{array}$ (21)

the solution of Equation (11) can be expressed in the form

$u\left(t\right)=G_0\left(t\right)b_0+{\displaystyle \sum _{k=1}^{\infty }{G}_k\left(t\right)b_k}.$ (22)

Once the coefficients $a_k$ with $k\ge 1$ are calculated and the step-size h has been carried out, the approximation of the solution of Equation (11) at the point $t=h$ , is given by

$w_1=G_0\left(t\right)b_0+{\displaystyle \sum _{k=1}^{\infty }{G}_k\left(t\right)b_k}.$ (23)

Once the approximation of the solution of Equation (11) at the point $t=ih$ is calculated, denoted by $w_i$ , verifying

$\begin{array}{l}{u}^{\prime }=Au+f\left(u,t\right),\\ u\left(t_i\right)=\left(\begin{array}{ccc}{u}_1\left(t_i\right)& \cdots & u_{n-1}\left(t_i\right)\end{array}\right)=w_i,\end{array}$ (24)

to obtain an approximation of the solution at the point $t=\left(i+1\right)h$ , a change to the independent variable $\tau =t-t_i$ is made, which transforms Equation (24) into

$\begin{array}{l}{u}^{\prime }=Au+f\left(u,\tau +t_i\right),\\ u\left(0\right)=\left(\begin{array}{ccc}{u}_1\left(0\right)& \cdots & u_{n-1}\left(0\right)\end{array}\right)=w_i,\end{array}$ (25)

which leads to the initial situation.

Calculating the coefficients of the development of:

$f\left(u\left(\tau \right),\tau +t_i\right)={\displaystyle \sum _{k=0}^{\infty }\frac{{\tau }^k}{k!}{c}_{k+1}}\text{with}{c}_{k+1}={\text{d}}^{k}g\left(0\right)/\text{d}{\tau }^k={\text{d}}^{k}g\left(ih\right)/\text{d}{t}^k,$ (26)

the following integration, that is to say the approximation of the solution at the point $t=\left(i+1\right)h$ , is obtained using the following algorithm:

$\begin{array}{l}{a}_0=w_i,\\ a_{k+1}=A{a}_k+c_k,k\ge 0,\\ b_0=a_0,\\ b_{k+1}=a_{k+1}-A{a}_k,k\ge 0,\\ w_{i+1}=G_0\left(t\right)b_0+{\displaystyle \sum _{k=1}^m{G}_k\left(t\right)b_k},\end{array}$ (27)

which constitutes an adaptation of the methods of series based on the G-functions for the numerical integration of parabolic PDEs in one space dimension.

3. Test Problems

We present the results of the implementation of the adapted algorithm to approximate the solution of two test problems selected from the literature.

3.1. Problem 1

The process of subjecting a bar to temperature change, is modeled using the equation in one-dimensional parabolic type partial derivatives (PDE), as it is assumed that the lateral temperature change is insignificant compared to the horizontal direction x.

In order to determine the temperature $u\left(x,t\right)$ at instant t and at point x, it is necessary to solve the following PDE

$\lambda u_t\left(x,t\right)=D{u}_{xx}\left(x,t\right)+F\left(x,t\right)$ (28)

where $\lambda$ is the heat capacity and D is the material diffusivity.

It is considered that the length of the bar is 1 unit, when $\lambda =1$ and $D=10$ , the following problem, proposed in [3] [16] , is obtained,

$u_t=10u_{xx}+10{\text{π}}^2\sin \left(\text{π}x\right),0<x<1,$ (29)

subject to the boundary conditions

$u\left(0,t\right)=u\left(1,t\right)=1,\text{0}<t\le 1,$ (30)

where initially:

$u\left(x,0\right)=1,\text{0}\le x\le \text{1,}$ (31)

whose exact solution is:

$u\left(x,t\right)=1+\left(1-{\text{e}}^{-10{\text{π}}^{2}t}\right)\sin \left(\text{π}x\right).$ (32)

It is considered a uniform mesh of space interval [0, 1], where $r=1/11$ and a uniform mesh in the time interval, where step size $h=0.005$ .

Given the Equation (3), then:

$\begin{array}{l}{\alpha }_0\left(t\right)=1,{\beta }_0\left(t\right)=0,{\gamma }_{\text{0}}\left(t\right)=1,\\ {\alpha }_1\left(t\right)=1,{\beta }_1\left(t\right)=0,{\gamma }_{\text{1}}\left(t\right)=1.\end{array}$ (33)

Applying the techniques shown in Equations (5)-(11), we get:

$\begin{array}{l}{u^{\prime }}_0\left(t\right)=u^{\prime }\left(0,t\right)={\left({\gamma }_0\left(t\right)/{\alpha }_0\left(t\right)\right)}^{\prime }=0,\\ {u^{\prime }}_{11}\left(t\right)=u^{\prime }\left(1,t\right)={\left({\gamma }_1\left(t\right)/{\alpha }_1\left(t\right)\right)}^{\prime }=0,\end{array}$ (34)

and the system of differential equations:

$\left(\begin{array}{c}{u^{\prime }}_1\\ {u^{\prime }}_2\\ ⋮\\ {u^{\prime }}_9\\ {u^{\prime }}_{10}\end{array}\right)=A\left(\begin{array}{c}{u}_1\\ u_2\\ ⋮\\ u_9\\ u_{10}\end{array}\right)+\left(\begin{array}{c}10\left({11}^2+{\text{π}}^2\sin \left(\text{π}/11\right)\right)\\ 10{\text{π}}^2\sin \left(2\text{π}/11\right)\\ ⋮\\ 10{\text{π}}^2\sin \left(9\text{π}/11\right)\\ 10\left({11}^2+{\text{π}}^2\sin \left(10\text{π}/11\right)\right)\end{array}\right),$ (35)

where

$A=\{\begin{array}{ll}-2\cdot {10}^2\cdot {11}^2,\hfill & i=j\hfill \\ 10\cdot {11}^2,\hfill & \left|i-j\right|=1,i,j=1,\cdots ,10\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}.$ (36)

In a more concise form:

$\begin{array}{l}{u}^{\prime }=Au+f\left(u,t\right),\\ u\left(0\right)=\left(\begin{array}{c}1\\ 1\\ ⋮\\ 1\\ 1\end{array}\right).\end{array}$ (37)

To calculate the approximations only two G-functions are used.

Given $u\left(t\right)$ is the solution to problem (37), then:

$u\left(t\right)=\left(\begin{array}{c}{u}_1\\ ⋮\\ u_{10}\end{array}\right)=\left(\begin{array}{c}{\displaystyle \sum _{k=0}^{\infty }{a}_{k,1}{t}^k/k!}\\ ⋮\\ {\displaystyle \sum _{k=0}^{\infty }{a}_{k,10}{t}^k/k!}\end{array}\right).$ (38)

Substituting Equation (38) into Equation (37) we get:

$\begin{array}{l}{a}_0=u\left(0\right),\\ a_1=\left(\begin{array}{c}{a}_{1,1}\\ a_{1,2}\\ ⋮\\ a_{1,9}\\ a_{1,10}\end{array}\right)=A\cdot \left(\begin{array}{c}{a}_{0,1}\\ a_{0,2}\\ ⋮\\ a_{0,9}\\ a_{0,10}\end{array}\right)+f\left(u,h\right),\end{array}$ (39)

which allows us to define the following coefficients:

$\begin{array}{l}{b}_0=a_0,\\ b_1=\left(\begin{array}{c}{b}_{1,1}\\ b_{1,2}\\ ⋮\\ b_{1,9}\\ b_{1,10}\end{array}\right)=f\left(u,h\right).\end{array}$ (40)

Once the value of the two G-functions has been found, we can see from $w_1$ the approximation of the solution $u\left(h\right)$ is $w_1=G_0\left(t\right)b_0+G_1\left(t\right)b_1$ .

Assuming the calculation of the approximation of the solution $u\left(ih\right)$ , which is noted for $w_i$ is calculated, the following steps of integration can be completed following the algorithm

$\begin{array}{l}{a}_0=w_i,\\ a_1=A{a}_0+c_0,\\ b_0=a_0,\\ b_1=a_1-A{a}_0,\\ w_{i+1}=G_0\left(t\right)b_0+G_1\left(t\right)b_1.\end{array}$ (41)

In Table 1 the maximum norm of errors is presented from the Adapted Series

and Crank-Nicholson methods where ${\Vert u_{exact}-u_{aprox}\Vert }_{x_i,\infty }=\underset{0\le t\le 1}{\max }\left|u\left(x_i,t\right)-w_i\left(t\right)\right|$ and ${\Vert u_{exact}-u_{aprox}\Vert }_{\infty }=\underset{0\le i\le 10}{\max }\underset{0\le t\le 1}{\max }\left|u\left(x_i,t\right)-w_i\left(t\right)\right|$ .

Figure 1 shows the absolute errors from the three methods: Crank-Nicholson, Dufort-Frankel and Adapted Series to PDE.

In Figure 2 and Figure 3 the absolute errors from the Crank-Nicholson, Dufort-Frankel and Adapted Series to PDE methods are compared vs. the exact solution of this problem.

For a better visualization of Figure 1, Figure 2 and Figure 3, the time scale is shown from $t=0.25$ .

3.2. Problem 2

Considering the problem proposed in [3] [17] [18]

$u_t\left(x,t\right)=a\cdot u_{xx}\left(x,t\right)+f\left(x,t\right),$ (42)

where

$f\left(x,t\right)=\left(r_2+2r_1^{2}a\tanh \left(r_1\left(x-1\right)+r_{2}t\right)\right){\text{sech}}^2\left(r_1\left(x-1\right)+r_{2}t\right),$ (43)

subject to the boundary conditions

$\begin{array}{l}u\left(0,t\right)=\tanh \left(r_{2}t-r_1\right),\\ u\left(1,t\right)=\tanh \left(r_{2}t\right),\text{0}<t\le 1,\end{array}$ (44)

and initial values

$u\left(x,0\right)=\tanh \left(r_1\left(x-1\right)\right),\text{0}\le x\le \text{1,}$ (45)

whose exact solution is:

$u\left(x,t\right)=\tanh \left(r_1\left(x-1\right)+r_{2}t\right).$ (46)

The solution is a wave travelling in the negative direction of the OX axis, where $r_1$ and $r_2$ are positive constants which determine the steepness and rate of propagation.

Considering a uniform grid of space intervals [0, 1], where $r=1/31$ and a uniform grid on the time interval, with step size $h=0.00025$ .

Applying the techniques presented from Equation (5) to Equation (11), we get:

$\begin{array}{l}{u^{\prime }}_0\left(t\right)=r_2{\mathrm{sech}}^2\left(r_{2}t-r_1\right),\\ {u^{\prime }}_{31}\left(t\right)=r_2{\mathrm{sech}}^2\left(r_{2}t\right),\end{array}$ (47)

and the system of differential equations:

$\left(\begin{array}{c}{u^{\prime }}_1\\ {u^{\prime }}_2\\ ⋮\\ {u^{\prime }}_{29}\\ {u^{\prime }}_{30}\end{array}\right)=A\left(\begin{array}{c}{u}_1\\ u_2\\ ⋮\\ u_{29}\\ u_{30}\end{array}\right)+\left(\begin{array}{c}\frac{a}{r^2}\tanh \left(r_{2}t-r_1\right)+{\theta }_1{\mathrm{sech}}^2\left(r_1\left(r-1\right)+r_{2}t\right)\\ {\theta }_2{\mathrm{sech}}^2\left(r_1\left(2r-1\right)+r_{2}t\right)\\ ⋮\\ {\theta }_{29}{\mathrm{sech}}^2\left(r_1\left(29r-1\right)+r_{2}t\right)\\ \frac{a}{r^2}{r}_2{\mathrm{sech}}^2\left(r_{2}t\right)+{\theta }_{30}{\mathrm{sech}}^2\left(r_1\left(30r-1\right)+r_{2}t\right)\end{array}\right),$ (48)

where

$\begin{array}{l}{\theta }_i=\left(r_2+2a{r}_1^2\tanh \left(r_1\left(i\cdot r-1\right)+r_{2}t\right)\right),\\ A=\{\begin{array}{ll}-2\cdot a\cdot {31}^2,\hfill & i=j\hfill \\ a\cdot {31}^2,\hfill & \left|i-j\right|=1\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array},\\ i,j=1,\cdots ,30.\end{array}$ (49)

In a more concise form $u^{\prime }=Au+f\left(u,t\right)$ , with initial conditions

$u\left(0\right)=\left(\begin{array}{c}\tanh \left(r_1\left(r-1\right)\right)\\ ⋮\\ \tanh \left(r_1\left(30r-1\right)\right)\end{array}\right).$ (50)

Two G-functions only are used in the calculation of the approximations and once the value of the G-functions is obtained, the integration is carried out in a similar way to the previous problem, using the algorithm described in Equation (27):

$\begin{array}{l}{a}_0=u\left(0\right),\\ a_1={\left(\begin{array}{ccccc}{a}_{1,1}& a_{1,2}& \cdots & a_{1,29}& a_{1,30}\end{array}\right)}^t\\ =A{\left(\begin{array}{ccccc}{a}_{0,1}& a_{0,2}& \cdots & a_{0,29}& a_{0,30}\end{array}\right)}^t+f\left(u,h\right),\\ b_0=a_0,\\ b_1={\left(\begin{array}{ccccc}{b}_{1,1}& b_{1,2}& \cdots & b_{1,29}& b_{1,30}\end{array}\right)}^t=f\left(u,h\right).\end{array}$ (51)

The approximation of the solution $u\left(h\right)$ , is $w_1=G_0\left(t\right)b_0+G_1\left(t\right)b_1$ .

Assuming the calculation of the approximation of the solution in $u\left(ih\right)$ , noted from $w_i$ , the following integration steps are completed by use of the algorithm for $r_1=1$ , $r_2=0.01$ and distinct values of a:

$\begin{array}{l}{a}_0=w_i,a_1=A{a}_0+c_0,\\ b_0=a_0,b_1=a_1-A{a}_0,\\ w_{i+1}=G_0\left(t\right)b_0+G_1\left(t\right)b_1.\end{array}$ (52)

Table 2 shows, according to the values of a, the maximum norm of errors

versus the exact solution, $\underset{0\le i\le 10}{\max }\underset{0\le t\le 1}{\max }\left|u\left(x_i,t\right)-w_i\left(t\right)\right|$ , of the Adapted Series

Method, of Crank-Nicholson and of Dufort-Frankel. It can see the error stability of the Adapted Method vs. the increase of numerical value of parameter a, moreover, the Adapted Method has a better performance to high values of this parameter.

4. Conclusion

A series method has been adapted for the numerical solution of PDEs. The Adapted Series Method to PDEs is based on the Method Of Lines (MOL) and on Scheifele’s G-functions. Numerical and computational algorithms have been designed. These numerical algorithms contain a simple algebraic procedure for the computation of the coefficients using recurrence relations which involve the perturbation function. Two test problems, proposed by various authors in different articles, have been resolved using the Adapted Series Method. The good performance of the Adapted Series Method is presented contrasting the results obtained himself and other known methods such as Crank-Nicholson and Dufort-Frankel, vs. exact solution.

Cite this paper

Cortés-Molina, M., Reyes, J.A. and García-Alonso, F. (2018) The G-functions Series Method Adapted to the Numerical Integration of Parabolic PDE. Journal of Applied Mathematics and Physics, 6, 161-173. https://doi.org/10.4236/jamp.2018.61016

References