Parabolic partial differential equation, numerical methods

Methods for solving parabolic partial differential equations on the basis of a computational algorithm. For the solution of a parabolic partial differential equation numerical approximation methods are often used, using a high speed computer for the computation. The grid method (finite-difference method) is the most universal.

As an example, the grid method is considered below for the heat equation

$$ \tag{1 } \frac{\partial u }{\partial t } = \ \frac{\partial ^ {2} u }{\partial x ^ {2} } + f( x, t),\ \ t > 0,\ 0 < x < l, $$

with boundary conditions of the first kind:

$$ u( 0, t) = \mu _ {1} ( t),\ \ u( l, t) = \mu _ {2} ( t),\ \ u( x, 0) = u _ {0} ( x). $$

One introduces a uniform grid of nodes $ ( x _ {i} , t _ {n} ) $,

$$ x _ {i} = ih,\ \ i = 0, \dots, N; \ \ hN = l,\ \ t _ {n} = n \tau ,\ \ n = 0, 1, \dots $$

and the notation

$$ y _ {i} ^ {n} = \ y( x _ {i} , t _ {n} ),\ \ y _ {t,i} = \frac{( y _ {i} ^ {n+ 1} - y _ {i} ^ {n} ) } \tau , $$

$$ y _ {\overline{x} ,i } ^ {n} = \frac{( y _ {i} ^ {n} - y _ {i- 1} ^ {n} ) }{h} ,\ y _ {x,i} ^ {n} = \frac{( y _ {i+ 1} ^ {n} - y _ {i} ^ {n} ) }{h} , $$

$$ y _ {\overline{x} x,i } ^ {n} = \frac{( y _ {i+ 1} ^ {n} - 2y _ {i} ^ {n} + y _ {i- 1} ^ {n} ) }{h ^ {2} } . $$

The grid method consists of approximating equation (1) by a system of linear algebraic equations (by a difference scheme)

$$ \tag{2 } \left . \begin{array}{c} y _ {t,i} ^ {n} = \ \sigma y _ {\overline{x} x,i } ^ {n+ 1} + ( 1 - \sigma ) y _ {\overline{x} x,i } ^ {n} + \phi _ {i} ^ {n} ,\ \ i = 1, \dots, N- 1 , \\ y _ {0} ^ {n} = \mu _ {i} ( t _ {n} ),\ \ y _ {N} ^ {n} = \mu _ {2} ( t _ {n} ),\ \ y _ {i} ^ {0} = u _ {0} ( x _ {i} ) , \\ \end{array} \right \} $$

where $ \sigma $ is a numerical parameter and the $ \phi _ {i} ^ {n} $ form a grid approximation of the function $ f( x, t) $, for example $ \phi _ {i} ^ {n} = f( x _ {i} , t _ {n} + 0.5 \tau ) $.

The system of equations (2) is solved in layers, i.e. for each $ n = 0, 1 \dots $ from the known values $ y _ {i} ^ {n} $ and $ \phi _ {i} ^ {n} $, the new values $ y _ {i} ^ {n+ 1} $ are calculated $ ( i = 1, \dots, N- 1) $. If $ \sigma = 0 $ (an explicit scheme), then the $ y _ {i} ^ {n+ 1} $ are expressed in an explicit manner in terms of the $ y _ {i} ^ {n} $ and $ \phi _ {i} ^ {n} $. If $ \sigma \neq 0 $ (an implicit scheme), then, with respect to the $ y _ {i} ^ {n+ 1} $, $ i = 1, \dots, N- 1 $, there arises a system of equations having a tridiagonal matrix. This system of equations is solved by the double-sweep method. The disadvantage of the explicit scheme is its strong restriction on the step $ \tau $, originating from the stability condition, namely $ \tau \leq 0.5 h ^ {2} $. The implicit schemes for $ \sigma \geq 0.5 $ are absolutely stable, that is, stable for arbitrary steps $ h $ and $ \tau $. Other grid methods for equation (1) are also known (see [1], [2]).

If the scheme (2) is stable and the $ \phi _ {i} ^ {n} $ approximate $ f( x, t) $, then for $ h, \tau \rightarrow 0 $ the solution $ y _ {i} ^ {n} $ of the difference problem converges to the solution $ u( x _ {i} , t _ {n} ) $ of the original problem (see [1]). The order of accuracy depends on the parameter $ \sigma $. Thus, a symmetric scheme ( $ \sigma = 0.5 $, $ \phi _ {i} ^ {n} = f( x _ {i} , t _ {n} + 0.5 \tau ) $) has second-order accuracy with respect to $ \tau $ and $ h $, i.e. for each $ n $ there is an estimate

$$ \max _ {0 < i < N } | y _ {i} ^ {n} - u( x _ {i} , t _ {n} ) | \leq M( \tau ^ {2} + h ^ {2} ) , $$

with $ M $ a constant not depending on $ h $ or $ \tau $. For $ \sigma = 0.5 - h ^ {2} /( 12 \tau ) $ and a special choice of the $ \phi _ {i} ^ {n} $ the scheme (2) has accuracy $ O( \tau ^ {2} + h ^ {4} ) $. For other values of $ \sigma $ the accuracy is $ O( \tau + h ^ {2} ) $.

For the solution of a parabolic partial differential equation on large intervals of time one essentially uses the asymptotic stability of the difference scheme. The solution of the equation (1) with $ f( x, t) = 0 $, $ \mu _ {1} ( t) = \mu _ {2} ( t) = 0 $ behaves, as $ t \rightarrow \infty $, like $ e ^ {- \lambda _ {1} t } $, where $ \lambda _ {1} = \pi ^ {2} /l ^ {2} $. Not every difference scheme for equation (1) possesses this property. For example, a symmetric scheme $ ( \sigma = 0.5 ) $ is asymptotically stable under the condition $ \tau \leq lh / \pi $. Schemes that have second-order accuracy and are asymptotically stable for arbitrary values of $ \tau $ and $ h $ have been constructed; however, they do not belong to the system (2) (see [3]).

Boundary conditions of the third kind,

$$ \frac{\partial u }{\partial x } ( 0, t) - \beta u ( 0, t) = \mu ( t), $$

are approximated up to order $ O( h ^ {2} ) $ by the difference equations

$$ \sigma ( y _ {x,0} ^ {n+ 1} - ( \beta y _ {0} ^ {n+ 1} + \mu ^ {n+ 1} )) + ( 1 - \sigma )( y _ {x,0} ^ {n} - ( \beta y _ {0} ^ {n} + \mu ^ {n} )) = $$

$$ = \ 0.5 h ( y _ {t,0} ^ {n} - \phi _ {0} ^ {n} ),\ \mu ^ {n} = \mu ( t _ {n} ). $$

For the heat equation in spherical or cylindrical coordinates,

$$ \frac{\partial u }{\partial t } = \ \frac{1}{r ^ {m} } \frac \partial {\partial r } \left ( r ^ {m} \frac{\partial u }{ \partial r } \right ) ,\ \ m = 1, 2,\ 0 < r < R, $$

one introduces the grid $ ( r _ {i} , t _ {n} ) $, where $ t _ {n} = n \tau $, $ n = 0, 1, \dots $

$$ r _ {i} = \left ( i + \frac{m}{2} \right ) h,\ \ i = 0, \dots, N; \ \ h = \frac{R}{N + 0.5 m } ,\ m = 1, 2, $$

and considers the difference equations (see [4])

$$ y _ {t,i} ^ {n} = \Lambda _ {r} ^ {( m)} ( \sigma y _ {i} ^ {n+ 1} + ( 1 - \sigma ) y _ {i} ^ {n} ) , $$

where

$$ \Lambda _ {r} ^ {( m)} y _ {0} ^ {n} = \ \frac{1}{r _ {0} ^ {m} } \left ( \overline{r} _ {1} ^ {m} \frac{y _ {1} ^ {n} - y _ {0} ^ {n} }{h} \right ) , $$

$$ \Lambda _ {r} ^ {( m)} y _ {i} ^ {n} = \frac{1}{r _ {i} ^ {m} } \frac{1}{h} \left ( \overline{r} _ {i+ 1} ^ {m} \frac{y _ {i+ 1} ^ {n} - y _ {i} ^ {n} }{h} - \overline{r} {} _ {i} ^ {m} \frac{y _ {i} ^ {n} - y _ {i- 1 } ^ {n} }{h} \right ) , $$

$ \overline{r} _ {i} = 0.5( r _ {i} + r _ {i+ 1} ) $ for $ m = 1 $, $ \overline{r} _ {i} = \sqrt {r _ {i} r _ {i- 1} } $ for $ m = 2 $,

$$ i = 1, \dots, N- 1 . $$

For solutions of parabolic partial differential equations with variable coefficients, homogeneous conservative difference schemes, expressing on the grid the principles of conservation inherent in the initial differential equation, are employed (see [1], [5]). For the construction of difference schemes for parabolic differential equations with variable coefficients the following methods are employed: balance, variational and finite-element methods. For example, for the equation

$$ \frac{\partial u }{\partial t } = \ \frac \partial {\partial x } \left ( k( x, t) \frac{\partial u }{\partial x } \right ) + f( x, t),\ \ 0 < x < l,\ t > 0, $$

difference schemes with weights are utilized:

$$ y _ {t,i} ^ {n} = ( a( \sigma y ^ {n+ 1} + ( 1 - \sigma ) y ^ {n} ) _ \overline{ {x}} ) _ {x,i} + \phi _ {i} ^ {n} , $$

where

$$ a = a _ {i} ^ {n} = k( x _ {i} - 0.5 h, t _ {n} + 0.5 \tau ). $$

Convergence has been proved and error estimates have been obtained for homogeneous conservative difference schemes for parabolic differential equations in the case of discontinuous coefficients as in the case of continuous coefficients (see [1]).

For the solution of systems of parabolic differential equations containing one spatial variable one uses the same difference schemes as for a single equation. The solution vector $ y ^ {n+ 1} $ for the new layer is obtained by the matrix factorization method.

For the solution of quasi-linear parabolic partial differential equations one uses implicit absolutely-stable difference schemes. The solution $ y ^ {n+1} $ of the difference equations is obtained by iteration using double-sweep. For example, for the equation

$$ \tag{3 } \frac{\partial u }{\partial t } = \frac \partial {\partial x } \left ( k( u) \frac{\partial u }{\partial x } \right ) ,\ \ k( u) > 0, $$

one uses the purely implicit difference scheme

$$ y _ {t,i} ^ {n} = ( a( y _ {i} ^ {n+ 1} ) y _ {\overline{x} } ^ {n+ 1} ) _ {x,i} , $$

$$ a( y _ {i} ^ {n+ 1} ) = 0.5 ( k( y _ {i} ^ {n+ 1} ) + k( y _ {i- 1} ^ {n+ 1} )), $$

which is solved by the iteration method ( $ s $ is the number of the iteration):

$$ y _ {i} ^ {( s+ 1) } = \ y _ {i} ^ {n} + \tau ( a( y _ {i} ^ {( s)} ) y _ {\overline{x} } ^ {( s+ 1)} ) _ {x,i } ,\ \ s = 0, \dots, m, $$

$$ y _ {i} ^ {( 0)} = y _ {i} ^ {n} ,\ y _ {i} ^ {( m+ 1)} = y _ {i} ^ {n+ 1} . $$

The solution $ y _ {i} ^ {( s+ 1)} $ for each iteration is obtained by the double-sweep method. For the solution of non-linear parabolic partial differential equations, an application has also been found for schemes based on the Runge–Kutta method. Thus, for equation (3) one uses the two-step method

$$ y _ {i} ^ {n+ 1/2} = \ y _ {i} ^ {n} + 0.5 \tau ( a( y _ {i} ^ {n} ) y _ {\overline{x} } ^ {n+ 1/2} ) _ {x,i} , $$

$$ y _ {i} ^ {n+ 1} = y _ {i} ^ {n} + 0.5 \tau ( a( y _ {i} ^ {n+ 1/2} ) ( y ^ {n+ 1} + y ^ {n} ) _ {x} ) _ {x,i} . $$

The testing and verification of the properties of difference schemes for non-linear parabolic partial differential equations is carried out by means of comparison with an exact auto-modelled solution, such as, for example, the solution of the type of a heat wave (see [6]).

For the solution of a multi-dimensional parabolic partial differential equation one uses a variable-directions method, which allows one to reduce the many variable problem to a succession of one variable problems (see [1], [7]–[10]). A significant number of absolutely stable algorithms of variable directions have been proposed and examined. These algorithms are economic in the sense that the number of arithmetic operations necessary for the calculations on the new time layer $ t _ {n+ 1} $ has the order of the number of nodes of the spatial grid. Below an example is considered of a scheme of variable directions for the equation

$$ \tag{4 } \frac{\partial u }{\partial t } = \ \frac{\partial ^ {2} u }{\partial x _ {1} ^ {2} } + \frac{\partial ^ {2} u }{ \partial x _ {2} ^ {2} } ,\ \ ( x _ {1} , x _ {2} ) \in G ,\ \ t > 0, $$

in a rectangle $ G = \{ {0 < x _ \alpha < l _ \alpha } : {\alpha = 1, 2 } \} $. One introduces the grid

$$ x _ {ij} = ( x _ {1} ^ {( i)} , x _ {2} ^ {( j)} ),\ \ x _ {1} ^ {( i)} = ih _ {1} ,\ \ i = 0, \dots, N _ {1} ; \ \ h _ {1} N _ {1} = l _ {1} , $$

$$ x _ {2} ^ {( j)} = jh _ {2} ,\ j = 0, \dots, N _ {2} ; \ h _ {2} N _ {2} = l _ {2} , $$

and the notation

$$ y _ {ij} ^ {n} = y( x _ {ij} , t _ {n} ),\ \ y _ {ij} ^ {n+ 1/2} = y( x _ {ij} , t _ {n} + 0.5 \tau ), $$

$$ \Lambda _ {1} y _ {ij} ^ {n} = \frac{( y _ {i+ 1},j ^ {n} - 2y _ {i,j} ^ {n} + y _ {i- 1},j ^ {n} ) }{h _ {1} ^ {2} } , $$

$$ \Lambda _ {2} y _ {ij} ^ {n} = \frac{( y _ {i,j+ 1} ^ {n} - 2y _ {ij} ^ {n} + y _ {i,j- 1} ^ {n} ) }{h _ {2} ^ {2} } . $$

Equations (4) are solved by the following difference scheme:

$$ \frac{y _ {ij} ^ {n+ 1/2} - y _ {ij} ^ {n} }{0.5 \tau } = \ \Lambda _ {1} y _ {ij} ^ {n+ 1/2} + \Lambda _ {2} y _ {ij} ^ {n} , $$

$$ \frac{y _ {ij} ^ {n+ 1} - y _ {ij} ^ {n+ 1/2} }{0.5 \tau } = \Lambda _ {1} y _ {ij} ^ {n+ 1/2} + \Lambda _ {2} y _ {ij} ^ {n+ 1} . $$

These equations are valid for $ i = 1, \dots, N _ {1} - 1 $, $ j = 1, \dots, N _ {2} - 1 $, and are supplied by the appropriate boundary conditions. By a double-sweep in the $ x _ {1} $-direction, for each $ j = 1, \dots, N _ {2} - 1 $, one obtains $ y _ {ij} ^ {n+ 1/2} $ from the first equation, and then by a double-sweep in the $ x _ {2} $-direction, for each $ j = 1, \dots, N _ {1} - 1 $, one obtains $ y _ {ij} ^ {n+ 1} $. In such a manner the computational algorithm consists of successive applications of one-variable double-sweeps.

The theoretical basis for the construction and analysis of economic difference schemes for multi-dimensional parabolic partial differential equations is the method of total approximation (see [1], [8], [9]). Additive economic schemes with equations on graphs and vector schemes appeared (see [11]) as generalizations of the variable-directions method.

In the case of irregular spatial grids, difference schemes for parabolic partial differential equations have been constructed using the method of finite elements and the balance method (see [12]).

For the solution of the grid equations arising from the approximation of multi-dimensional parabolic partial differential equations by implicit difference schemes, the effective direct methods and iterative methods worked out for elliptic difference boundary value problems are employed: the method of separation of variables, the algorithm of the fast discrete Fourier transform, the method of cyclic reduction, the alternately triangular iteration method, and others (see [13]).

References

[1]	A.A. Samarskii, "Theorie der Differenzverfahren" , Akad. Verlagsgesell. Geest u. Portig K.-D. (1984) (Translated from Russian)
[2]	V.K. Saul'ev, "Integration of equations of parabolic type by the method of nets" , Pergamon (1964) (Translated from Russian)
[3]	A.A. Samarskii, A.V. Gulin, "Stability of difference schemes" , Moscow (1973) (In Russian)
[4]	I.V. Fryazinov, "Difference schemes for Poisson's equation in polar, cylindrical and spherical coordinate systems" USSR Comp. Math. Math. Phys , 11 : 5 (1971) pp. 153–165 Zh. Vychisl. Mat. i Mat. Fiz. , 11 : 5 (1971) pp. 1219–1228
[5]	G.I. Marchuk, "Methods of numerical mathematics" , Springer (1982) (Translated from Russian)
[6]	A.A. Samarskii, I.M. Sobol', "Examples of the numerical calculation of temperature waves" USSR Comp. Math. Math. Phys. , 3 : 4 (1963) pp. 945–970 Zh. Vychisl. Mat. i Mat. Fiz. , 3 : 4 (1963) pp. 702–719
[7]	N.N. Yanenko, "The method of fractional steps; the solution of problems of mathematical physics in several variables" , Springer (1971) (Translated from Russian)
[8]	A.A. Samarskii, "Homogeneous difference schemes for non-linear parabolic equations" USSR Comp. Math. Math. Phys. , 2 : 1 (1962) Zh. Vychisl. Mat. i Mat. Fiz. , 2 : 1 (1962) pp. 25–56
[9]	N.N. Yanenko, "Weak approximation of systems of differential equations" Sibirsk. Mat. Zh. , 5 : 6 (1964) pp. 1431–1434 (In Russian)
[10]	E.G. D'yakonov, "Difference schemes with a "disintegrating" operator for multidimensional problems" USSR Comp. Math. Math. Phys. , 2 : 4 (1962) pp. 581–607 Zh. Vychisl. Mat. i Mat. Fiz. , 2 : 4 (1962) pp. 549–568
[11]	A.A. Samarskii, I.V. Fryazinov, "Difference approximation methods for problems of mathematical physics" Russian Math. Surveys , 31 : 6 (1976) pp. 179–213 Uspekhi Mat. Nauk , 31 : 6 (1976) pp. 167–196
[12]	I.V. Fryazinov, "Balance method and variational-difference schemes" Differential Eq. , 16 : 7 (1980) pp. 853–862 Differentsial'nye Uravneniya , 16 : 7 (1980) pp. 1332–1343
[13]	A.A. Samarskii, E.S. Nikolaev, "Numerical methods for grid equations" , 1–2 , Birkhäuser (1989) (Translated from Russian)
[a1]	G.D. Smith, "Numerical solution of partial differential equations: finite difference methods" , Clarendon Press (1978)
[a2]	C. Johnson, "Numerical solution of partial differential equations by the finite element method" , Cambridge Univ. Press (1987)
[a3]	A.R. Mitchell, D.F. Griffiths, "The finite difference method in partial differential equations" , Wiley (1980)
[a4]	L. Lapidus, G.F. Pinder, "Numerical solution of partial differential equations in science and engineering" , Wiley (1982)
[a5]	A.R. Gourlay, "Hopscotch, a fast second order partial differential equation solver" J. Inst. Math. Appl. , 6 (1970) pp. 375–390
[a6]	A.W. Peaceman, H.H. Rachford, "The numerical solution of parabolic and elliptic differential equations" SIAM J. , 3 (1955) pp. 28–41