Application of the galerkin method to the solution of burgers' equation

COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING 14 (1978) 305-321 0 NORTH-HOLLAND PUBLISHING COMPANY

APPLICATION OF THE GALERKIN METHOD TO THE SOLUTION OF BURGERS’ EQUATION A.M. DAVIES Institute of Oceanographic Sciences, Bidston Laboratory, Birkenhead, Merseyside L43 7RA, England

Received 22 August 1977 The Galerkin method, with an expansion of B-splines is used to solve Burgers’ equation from an initial state involving a discontinuity in the space domain. In an initial series of calculations a fixed spatial distribution of knots is used to define the B-splines, and the effect upon accuracy of changes in the order of the spline and the number of terms in the expansion is examined. In a second series of calculations the spatial distribution of knots is varied as the solution advances through time so that an area of tine knot resolution is maintained in the region where the solution of Burgers’ equation changes most rapidly. Results of high accuracy are obtained using this method, which avoids the necessity of having a fine resolution knot distribution over the whole spatial domain in order to achieve an accurate solution.

1. Introduction In a previous paper [ 11 solutions of Burgers’ equation for a number of initial conditions corresponding to waves of differing wavelength were calculated using the Galerkin method. The basis functions used were B-splines defined upon a uniformly spatially distributed set of knots, and the effect upon accuracy of changing the order of the spline and the knot spacing was examined. In this paper Burgers’ equation is again used as a one-dimensional model of the shallow water equations, althoughnow the initial condition is of the form of a shock. The solution of Burgers’ equation for this initial condition, represents a shock which progresses through space with a uniform velocity. The width of the shock is determined by the coefficient of diffusivity specified in the equation. Burgers’ equation is a simple representative model for nonlinear shallow water equations, having both an advective term and a dispersive term. The example considered here of a progressive shock wave provides a rigorous test of how accurately the Galerkin method with a basis of B-splines can be expected to calculate the steepening of a water wave moving into shallow water. In this motion the wave front gradually steepens, due to the nonlinear term in the shallow water equations, and energy is transferred to shorter wavelengths. Eventually the wavefront becomes vertical, and a bore is formed. The numerical solutions calculated here with the Galerkin method, considering a range of dispersive coefficients, are compared with the analytical solution of Lighthill 121. A number of calculations were performed in which the order of spline and knot spacing were varied; the effect of these changes upon accuracy is examined. Results of the calculations showed that the accuracy of the solution can be improved if a finer knot spacing is used in the region of the shock.

A.M. Davies, Application

306

of the Galerkin method to the solution of Burgers’ equation

A method has been developed which permits the area of the finer knot spacing to move along with the region of the shock, the former thereby continually embracing the latter. Results using this method are presented and compared with those obtained with a uniform knot distribution.

2. Numerical solution of Burgers’ equation Burgers’ equation is a simple one-dimensional model of the Navier-Stokes into account both convection and diffusion and is given by

au

gl=_ctl

is-*+ u ax

ax2

’

x
0,

x>o

taking

v>o,

(1)

where v is the coefficient of diffusivity, For the initial conditions A,

equations

u is the velocity

in the x direction,

and t is the time.

(2)

U(X, 0) =

an exact solution

of equation

(1) is given by Lighthill

[ 21, namely

This solution was used in determining the errors arising in applying the Galerkin solution of Burgers’ equation. In the numerical examples considered later the boundary conditions

method

to the

~(0, t) = u,, (t) = A ,

(4) u(l,t)=u,(t)=O were used. Numerical solutions of eq. (1) were derived using the Galerkin splines. An expansion for u(x, t) of the form u(x, t) =

5

method

with a basis set of B-

a,(t) M,,(x)

i=l

was used, where M,,(x) is a B-spline of order ~1,a,(t) are time-dependent coefficients, and m is the number of basis functions. Substituting eq. (5) into (1) and using the Galerkin method gives the set of coupled differential equations

A.M. Davies, Application

ai i

= ’ fi

$

of the Galerkin method to the solution of Burgers’ equation

(“nitx)) Mnk(X) diX

9

k=

1,2,...m.

307

(6)

In the numerical examples presented later the order y1of the spline took values 2,4 and 6. For y1= 2 the B-spline is a linear function; thus the second derivative is zero. Hence, in order to solve eq. (6) with second-order splines as basis functions, the integral involving the second derivative is integrated by parts. Thus

s’

o

-d2 dx*

(,Mi(X))M,(X)dX=-&(Mi(X))M*(X)l10

’

-J

d

a”i(x)$Mk(X)dY’

(7)

0

In this equation the order y1of the spline has been dropped from the notation. Substituting eq. (7) into (6) and expressing (6) in matrix form gives C*=vAa-vEa-Da, dt

(8)

where C, A, E and D are m X m matrices with elements

M,tx) Mitx) dX

cik = j

7

Ajk =“ktX)

0

’

E, = 1

d

z

M,(X)

& M,(x) do

7

D, =

0

2 Mitx)ll 0 )

,c

(9)

ai i

M,(x)

0

where a is the column vector {a, a, . . . a,,, }.

3. Form of the basis functions

The ease with which the integrals given by (9) can be evaluated depends upon the nature of the basis set. Spline functions can be readily integrated and differentiated because they are piecewise polynomials, and since they have a basis with small support, many of the integrals that occur are zero. The basis functions used here are B-splines define upon an increasing set of N interior knots A. Also 2n additional knots outside the region 0 < x < 1 are required for support, positioned at

x,_,Gl,_,c..
(104

308

A.M. Davies, Application

of the Galerkin method to the solution of Burgers’ equation

and

(lob)

1 = hiv+i < AN+2 . . . < hN+n .

Thus m, the number of B-splines in expansion (S), is equal to N + ~1,and for a uniformly spatially distributed set of knots AX the knot separation is for the region 0 < x < 1 considered here equal to I/h, whereh =N+ 1. A general B-spline M,(x) of order n is defined upon knots hi_,,, &++r, . . . hi with the restricted support property

M,(X)

>o,

hi_n
= 0,

x <

(11)

xi_,

)

x > hi .

Hence Mi(x) is only nonzero over y1knot intervals. The spline functions were evaluated using’the numerically stable method given by Cox [31. This method has the advantage that the knot spacing is completely arbitrary and allows for multiple knots. Rather than evaluating the integrals given in (9) directly, the splines were expanded in terms of Chebyshev polynomials of the first kind, and the integrals were evaluated from these expansions. Details of the method have been given previously [ 11 and will not be discussed again here. Initial values of the coefficients a&O) required in thedolution of eq. (8) were calculated by expanding the initial condition u(x, 0) given in eq. (2) in terms of B-splines using (5) and integrating over the spatial domain to give

,$ai

i

Mi(X) M,(X) do = i

0

0

U(X,0) M,(X) do

7

k= 1, 2, . . . m .

(12)

For the initial conditions used, all the integrals arising in (12) were evaluated analytically. Using a basis of B-splines, boundary conditions can be readily implemented. Specifying knots having a multiplicity of n - 1 at x = 0 and x = 1, all the B-splines in expansion (5) except the first and last are zero at these points. Thus eq. (5), incorporating the boundary conditions given in (4), is given by

U(X,

m-l

110(0

t) = I_ 4

M,(X)

+

C

i=2

aiMi

u,(t)

+A

,

(13)

m

where A, = M, (0), A, = Mm (1), and the unknown coefficients a2, as, . .. a, _-1are determined as described previously using eq. (8), though now with k = 2, 3, ,.. m - 1. Before solving the set of coupled differential equations.given in (8), using an explicit time stepping technique, the matrix C is inverted using any standard numerical technique. The object of the results presented later is to show how the accuracy of the numerical solution is effected by the order of the spline, the number of basis functions and the knot separation. Hence, to keep errors arising

A.M. Davies, Application of the Galerkin method to the solution of Burgers’equation

309

from solving these equations at discrete time steps from masking those arising using different orders of spline, a fourth-order Runge Kutta method was used. An initial series of calculations using both second and sixth order B-splines with k = 25 was performed in which At took values 0.10, 0.05, 0.02, 0.01 and 0.005. Numerical solutions calculated using time steps 0.02, 0.01 and 0.005 agreed to seven significant figures, and a value of 0.02 was subsequently used in all calculations in order to minimize computer time, though still retaining seven-figure accuracy.

4. Numerical calculations using a time independent

knot spacing

In this section results are presented coming from the numerical solution of Burgers’ equation involving a number of values of the diffusion coefficient v. Various orders of spline and numbers of terms m in the expansion were used. Initially a uniform knot spacing was employed, although in a later series of calculations a much finer knot spacing was used in the vicinity of the shock.

Fig. 1. Comparison of analytical and numerical solutions at time t = 0 and t = 0.5: tion (h = 25, n = 2, v = 0.0075).

analytical solution, ---

numerical solu-

The analytical solution of Burgers’ equation given previously (3) represents a shock wave progressing along the x axis, the shock front travelling with a velocity U = 0.5. The analytical solution has a maximum value of unity behind the shock and a minimum value of zero in front of the shock (see plot of analytic solution in fig. 1). The numerical approximation of the initial state in terms of an expansion of the B-splines considered here oscillates about the true initial state (see plot of numerical approximation of the initial state in fig. 1). However, as Burgers’ equation is integrated forward in time, this oscillation should diminish in magnitude, and by the time t = 0.5 the initial wave form should have evolved into a steady state profile moving in the positive x direction with a mean shock speed U = 0.5. The width of the shock at time t = 0.5 depends upon the value of the coefficient of diffusion specified in the solution of Burgers’ equation. The higher the value of this coefficient the greater the width. Comparison of the width, on the one hand calculated numerically and on the other hand derived from the analytic solution, is therefore particularly important as

A.M. Davies, Application of the Galerkin method to the solution of Burgers’ equatiorz

310

indicating whether diffusive properties are being correctly reproduced in the numerical solution (Roache 141). Roache [5, 61 has shown for the space-centred finite difference method that if the local cell Reynolds number Re, = UAx/v > 2 (where U is the speed of propagation of the shock, in this case U = 0.5 and Ax is the grid spacing in the finite difference scheme), then a spurious oscillation of the solution occurs behind the shock front, and the numerical value of this solution can exceed unity. Taylor, Ndefo and Masson [ 71 in an extensive study of the application of various finite difference schemes to the solution of Burgers’ equation gave a number of examples in which these oscillations were present when Re, exceeded 2. Hence, in comparing various numerical results the magnitude of this oscillation is an important parameter. In the tables of numerical results presented later three important parameters are used to compare the accuracy of the solutions, namely: (i) the maximum value of the solution (the extent to which this exceeds unity, gives the error due to the spurious oscillation of the solution behind the shock front), (ii) the speed of propagation of the shock (the deviation from 0.5 gives errors in the propagation speed of the shock), (iii) the width of the shock profile (changes in its value may be due to artificial viscosity (Roache [4] )). In the numerical calculation an initial state given by i

1

3

0 < x < 0.4

U(X, 0) =

(14) 0,

0.4
1.0

was used together with boundary conditions given by (4), and Burgers’ equation was integrated forward through time to t = 0.5. By choosing this initial state and period of integration, the shock travelled from x = 0.4 to x = 0.65, and hence throughout the integration period it was a reasonable distance from any effect upon the solution imposed by the boundaries at x = 0 and x = 1. As stated previously, when this initial condition is approximated using the expansion of B-splines given in (S), the resulting approximation exhibits an oscillatory behaviour on both sides of the discontinuity (see the numerical approximation (h = 25, y1= 2) of the initial state in fig. 1). The maximum difference between the numerical approximation (for h = 25) and the true function was 0.1340 for second-order splines, 0.1009 for fourth-order splines and 0.0963 for sixth-order splines. Although the maximum of the spurious oscillation was diminished for the higher order splines, the region over which the oscillation was present was larger than for the lower order splines. The presence of this initial oscillation is useful in that if its magnitude diminishes as Burgers’ equation is integrated forward in time, the numerical solution will probably be meaningful. An increase in magnitude of the spurious oscillation with time yielding in any long time period integration to a solution in which this oscillation masks the true solution. In an initial series of numerical calculations a uniform knot distribution was used, and distance Ax between adjacent interior knots was decreased by increasing the number of knots h in the region 0 G x Q 1. The calculations were performed on an IBM 370/165 computer using double precision arithmetic (16 figures of accuracy). Numerical solutions were calculated for time t = 0.5 covering a range of values of v and orders of spline n. In order to determine the accuracy of the numerical solution, values of u were calculated at

A.M. Davies, Application of the Galerkin method to the solution of Burgers’ equation Table 1. Numerical solutions of Burgers’ equation calculated with a range of values of v at time t = 0.5 with a uniform knot interval Order of spline

h

Maximum value of solution u

Width of shock in units of l/200

Speed of propagation U

Calculated results with v = 0.00375 2

10 15 20 25

1.2750 1.2524 1.2370 1.1600

16 16 15 12

0.5070 0.4912 0.4937 0.5017

4

10 15 20 25

1.1412 1.1364 1.1297 1.1148

16 14 14 12

0.5041 0.4950 0.4978 0.4988

6

10 15 20 25

1.1416 1.1300 1.1040 1.0828

16 14 13 12

0.5036 0.4959 0.4990 0.4997

1.0000

13

0.5000

Analytical solution

Calculated results with a~= 0.0075 2

10 15 20 25

1.0993 1.0542 1.0454 1.0200

24 24 25 25

0.4918 0.4958 0.4989 0.499 1

4

10 15 20 25

1.0599 1.0592 1.0316 1.0143

22 24 24 25

0.5020 0.4974 0.4984 0.4990

6

10 15 20 25

1.0633 1.0493 1.0071 1.0071

22 22 25 25

0.5010 0.4977 0.4994 0.4994

1.0000

27

0.5000

Analytical solution

Calculated results with v = 0.01 2

10 15 20 25

1.0914 1.0203 1.0047 1.0031

39 33 34 35

0.494 1 0.4974 0.4992 0.4994

4

10 15 20 2s

1.0262 1.0199 1.0126 1.0041

33 33 38 38

0.5010 0.4978 0.4991 0.4996

6

10 15 20 25

1.0174 1.0187 1.0106 1.0029

34 34 38 37

0.5008 0.4986 0.4992 0.4997

1.0000

31

0.5000

Analytical solution

311

312

A.M. Davies, Application

of the Galerkin method to the solution of Burgers’ equation

time t = 0.5 at 200 equally spaced points in the region 0 S x < 1, and the maximum values of u are given in the tables. The difference between this numerical value and unity is the magnitude of the spurious oscillation behind the shock. The width of the shock was also calculated. In this investigation it was taken as the distance along the x axis (in units of l/200) over which the value of u diminished from 0.99 to 0.01. The propagation speed of the shock was also calculated. In table 1 numerical results are given at t = 0.5 for second, fourth and sixth order B-splines defined upon a uniform set of knots having four different knot spacings Ax determined by Ax = l/h. Three values of v were used, namely 0.00375, 0.0075 and 0.01. It is evident from the result in table 1 for the case v = 0.00375 that the maximum value of u calculated using second order splines, even with a fine knot spacing (h = 25) has a significant numerical error. Comparing results for a constant value of h, it is evident that an improvement in accuracy arises as the order of spline is increased, and a sixth-order spline with h = 25 accurately reproduces both the width of the shock and its speed of propagation although the magnitude of the spurious oscillation is still significant. The maximum value of u for all cases except y2= 6, h = 2.5 increased from its initial value. Increasing v from 0.00375 to 0.0075 reduces Re, by half for a given h, and it is evident from table 1 that the maximum value of u for a given h and n is also reduced. The maximum u is now less than its initial value in all cases. The improvement in accuracy obtained using higher order splines is no longer as significant as previously. Results obtained with v = 0.01 are very similar for all orders of spline although the higher order splines are still slightly more accurate. The maximum value of u has been considerably reduced, and it is evident that with h = 25 (Re, = 2) the magnitude of the spurious oscillation is below 0.005, and the speed of propagation and the width of the shock are accurately calculated. It is obvious from the results presented in table 1, that for a given v, as the knot separation is decreased, the magnitude of the spurious oscillation behind the shock diminishes, as is evident from the reduction in the maximum value of U. Also, for a given uniform knot distribution, increasing the viscosity diminishes the magnitude of the oscillation. It is evident from these results that the spurious oscillation will only diminish as the cell Reynolds number Re, (the cell length Ax was taken as the distance between adjacent knots) is reduced. A value of Re, < 2 appears necessary to remove the spurious oscillations. This restriction on Re, is consistent with that required using centred space finite difference techniques (Roache [ 5, 61). Zienkiewicz and Godbole [8] have shown that using the standard Galerkin method with a basis of piecewise functions results in expressions that are similar to those of central differences, and the same restriction upon Re, is required. The reduction in magnitude of the spurious oscillation for higher order splines and their increased accuracy (compare second and fourth order results with h = 25 in table 1) suggest that they are preferable as basis functions at low v values. From the work of Zienkiewicz and Godbole [8] and the results presented in table 1 it is evident that - provided v and h are such that Re, < 2 - accurate solutions can be obtained. However, for low v values this requires the use of a large number of spline fucntions, which is computationally uneconomic. The remainder of this paper deals with obtaining accurate solutions for low values of u (i.e. Y = 0.0075 and 0.00375) without the need to increase h (i.e. h never exceeds 25). One way of reducing Ax in the region of the shock is to use an unequal knot spacing, and in a second series of numerical calculations a nonuniform knot spacing was employed, in which a much finer knot resolution was used in the region river which the shock would propagate. The knot

A.M. Davies, Application of the Galerkin method to the solution of Burgers’ equation

313

Table 2. Numerical solutions of Burgers’ equation calculated with a range of values of u at time t = 0.5 with a nonuniform knot interval --Order of spline

u

Maximum value of solution u

Width of shock in units of l/200

Speed of propagation U

0.00375

1.046 1 1.0308 1.0293

10 10 12

0.4973 0.4976 0.499 1

0.0075

1.0004 1.0001 1.0008

28 28 26

0.4990 0.5003 0.4996

0.01

1.0005 1.0000 1.0003

37 38 36

0.4995 0.5000 0.5000

spacing in the x domain was 0.0 (0.1) 0.3 (0.05) 0.40 (0.025) 0.7 (0.05) 0.8 (0.1) 1.0, with extra knots outside the region 0 < x < 1 for support. Thus h = 2 1, although Ax = 0.025 in the region through which the shock propagated (i.e. 0.4 < x < 0.65) and had a maximum value Ax = 0.1 outside that region. The maximum values of U, together with the width and speed of propagation of the shock obtained from the numerical solution of Burgers’ equation, are given in table 2 for v = 0.00375, 0.0075 and 0.01 at time t = 0.5. It is evident, comparing the results presented in table 2 with those of table 1, that the maximum value of the spurious oscillation in the rear of the shock is reduced for a given value of u, and the width of the shock and its speed of propagation are also slightly improved. Again the accuracy of the solution improves for the case v = 0.00375 as the order of spline is increased. From the results it is evident that using a finer knot resolution in the vicinity of the shock produces results of greater accuracy with a smaller number of basis functions (i.e. compare results in table 2 with those in table 1 for h = 25) than those obtained with a uniform knot resolution throughout the whole region. Obviously, in any long time integration it would be necessary to move the fine resolution region through space as the shock front advances in order to keep the area of finer resolution over the region of the shock. Details of a method of accomplishing this along with the results obtained are given in the next section.

5. Numerical calculations using a time-dependent

knot distribution

In meteorology a number of finite difference techniques using finite difference grids of different resolution (Hill [ 91, Phillips and Shukla [ 101, Harrison and Elsberry [ 111) have been developed in connection with hurricane prediction, in which a fine mesh positioned over the centre of the hurricane moves through a coarser mesh. By this technique the centre of the hurricane can be accurately resolved without the need to have a fine grid over the whole region. A similar idea can be applied here by changing the knot spacing as the shock progresses along the x axis.

A.M. Davies, Application

314

of the Galerkin method to the solution ofBurgers’

equation

Fig. 2. Distribution C (see table 3) of knots at time t = T, at which time the shock has moved by two knot spaces (arrow indicates position of centre of shock).

In the series of calculations now described only two regions of different resolution were considered: a coarse mesh region in front of and behind the shock with uniform knot spacing (region Fi, in fig. 2) and a region having a fine uniform knot spacing situated over the shock itself (region R,), the centre of the latter region moving with the shock. When the centre of the shock (i.e. the point at which u = 0.5) has moved a distance 2Ah (Ah is the distance between two adjacent knots in the higher resolution region), then the knot positions are redefined by increasing by one the number of knots in the region in front of the shock and reducing by one the number of knots in the region behind the shock. This redefinition of knots is performed in such a way that the uniform knot spacings in regions R, and R, are preserved. This rearrangement of knot positions corresponds to moving the region R, by a distance 2Ah through the region R, (note the difference in knot positions at time t = T compared with those at time t = 0 in fig. 2). Hence, starting with an initial knot distribution, denoted as distribution A, a set of B-splines were calculated and coefficients in the expansion computed from

g

a?)(O)

j

Mi’A)(x)MiA)(x) dx = j

0

0

u(x, 0)

M?)(x) dx

,

k = 1, 2, ... m

,

(15)

where the superscript (A) denotes that the B-splines were defined upon distribution A of knots, and a?) are the coefficients based on this distribution. Using eq. (8), the coefficients aiA) were advanced through time, and the value of u at the knot positioned 2Ah from the knot corresponding to the centre of the shock at t = 0 was calculated at each time step. When the value of u at this knot equalled or exceeded 0.5, the centre of the shock was judged to have moved over a distance 2Ah and the knot positions were redefined as described previously. Denoting this new arrangement of knots as distribution (B), new coefficients @j(T) were then calculated from, m

m

c i=1

ai’“)
.

(16)

It appeared reasonable to calculate the a{@(T) from the a{‘)(T) in a least squares sense since this is consistent with the Galerkin method. Hence, the a{*)(T) were calculated from,

g

@‘CT)j 0

MicB’(x>

M,E)(x) dx = &?)(7’)

J MfA)(x) MiE)(x) dx , 0

k = 1, 2, ... in .

(17)

315

A.M. Davies, Application of the Galerkin method to the solution of Burgers’ equation

An alternative method is to make the two expansions given in (16) agree at m specific x points, say xi, giving

(18) This method would appear appropriate if a collocation technique were used to solve Burgers’ equation. The xi could then be taken as the collocation points. In practice eq. (17) was solved for the new set of coefficients @I, and then eq. (8) with matrices C, A, E, and D containing integrals of the B-splines defined upon the new knot distribution B was used to advance the coefficients through time. Provided a constant knot separation is maintained in the coarse region R, and fine region R, (see fig. 2), it is not necessary to calculate any new integrals whenever a new knot distribution is used. The movement of the tine resolution region R, through the region R, only involves changing the positions of the integrals within the matrices C, A, E and D (see appendix) and does not require the evaluation of new integrals. If, however, the tine mesh region reaches the boundary at x = 1, then a few integrals involving B-splines based upon the knot at x = 1 have to be recalculated. In the examples considered in the present paper this situation never arose. Although no new integrals are required in forming the various matrices, it is necessary to compute a new inverse for the matrix C whenever a new knot distribution is required if the method described previously for the solution of eq. (8) is used. In any long-time integration an iterative solution of the set of coupled differential equations given in (8), avoiding the inversion of the matrix C, would probably be more economical in computer time. New coefficients @)(7’) would also be calculated iteratively from eq. (17), again avoiding the inversion of matrix C. In order to test the accuracy of the method, Burgers’ equation was integrated through time to t = 0.5 subject to initial conditions (14) corresponding to a discontinuity at x = 0.4 and the boundary conditions given by (4) for v = 0.00375 and 0.0075. Six initial knot distributions (namely distributions A to F, see table 3) were used to test the effect upon the accuracy of the solution of refining the knot resolution in the fine and coarse resolution regions and varying the extent of the region of finer resolution. For each of these distributions the region of finer resolution was maintained over the shock as the integration proceeded through time, using the technique described previously. The change in the accuracy of the solution arising from increasing the order of the spline was also examined for each distribution. Initially, the effect of increasing the resolution of the higher resolution region was examined Table 3. lnitial knot positions used in solving Burgers’ equation with a time-dependent knot distribution Knot distribution

Number of knot intervals k

Knot positions

A B c D E F

12 16 24 23 23 22

0.0 0.0 0.0 0.0 0.0 0.0

(0.1) (0.1) (0.05) (0.05) (0.05) (0.05)

0.3 0.3 0.3 0.3 0.35 0.3

(0.05) 0.5 (0.025)0.5 (0.025) 0.5 (0.025) 0.45 (0.025) 0.5 (0.025) 0.4

(0.1) (0.1) (0.05) (0.05) (0.05) (0.05)

1.0 1.0 1.0 1.0 1.0 1.0

316

A.M. Davies, Application of the Gale&in method to the sohrtiorl of Burgers’ equation Table 4. Numerical solutions of Burgers’ equation using a time-dependent Order of spline

h

2 2 2 2 2 2

12 16 24 23 23 22

4 4 4 4 4 4

12 16 24 23 23 22

Knot distribution

A B c D E F

Maximum

_-

knot distribution ___ _-.

for v = 0.00375

24

Width of shock

Speed of propagation U

1.2048 1.1007 1.0473 1.0475 1.0543 1.0755

17 15 10 10 12 13

0.4980 0.4980 0.4978 0.4977 0.4975 0.4975

1.1974 1.0963 1.0337 1.0331 1.0417 1.0703

17 15 10 10 10 12

0.4984 0.4986 0.4998 0.4997 0.4988 0.4984

(compare distributions B and A in table 3), and subsequently the effect of increasing the resolution of the coarser resolution region (compare distributions C and B) was considered. The change in the accuracy of the solution arising from reducing the extent of the finer resolution region in front of the shock (compare distributions C and D) and from reducing its extent behind the shock (compare distributions C and E) was also examined. In a final claculation the effect of having only the coarse resolution region in front of the shock was considered (distribution F). Comparing results for knot distributions A, B and C given in tables 4 and 5, it is evident that reducing the knot separation in the fine region (compare results A with B) improves the accuracy of the solution. However, a further refinement of the coarser region (compare results B with C) produced a further improvement in accuracy. Comparing the solution calculated with secondorder splines (knot distribution C, v = 0.0075, h = 24) with that obtained using a uniform knot distribution (V = 0.0075, h = 25, II = 2 in table I), it is evident that the magnitude of the maximum oscillation has been reduced significantly even though the number of spline functions has been reduced by one. A similar increase in accuracy occurred for v = 0.00375. Results obtained using Table 5. Numerical solutions of Burgers’ equation using a time-dependent

knot distribution for v = 0.0075

Order of spline

h

Knot distribution

Maximum u

Width of shock

Speed of propagation U

2 2 2 2 2 2

12 16 24 23 23 22

A B c D E F

1.2735 1.1041 1.0001 1.0001 1.0092 1.0090

22 24 28 28 27 29

0.4980 0.4979 0.4996 0.4996 0.4995 0.4995

4 4 4 4 4 4

12 16 24 23 23 22

A B c D E F

1.0560 1.0051 1.0001 1.0002 1.0024 1.0023

22 25 27 28 28 29

0.4992 0.4994 0.4999 0.4998 0.4999 0.4996

A.M. Davies, Application of the Galerkin method to the solution of Burgers’ equation

317

distribution C with fourth-order splines were slightly more accurate than those calculated with second-order splines. Results similar to those obtained with distribution C were calculated using distribution D in which the region of finer resolution was reduced in front of the shock. However, reducing the region of finer resolution by a similar amount (although in the rear of the shock) produced a deterioration in the accuracy of the solution (distribution E, tables 4 and 5). Limiting the area of finer resolution to the region behind the centre of the shock (distribution F) gave similar results for ZJ= 0.0075 to those obtained with distribution E, although for v = 0.00375 a deterioration in accuracy is evident. Solutions calculated with a finer resolution region which was not centred over the shock (i.e. distributions D, E and F) were slightly asymmetric with respect to the centre of the shock, and for this reason distribution C is preferable to distribution D even though the results presented in table 4 and 5 are of comparable accuracy. From the results it is evident that using a higher resolution region moving within one of coarser resolution reduced the magnitude of the spurious oscillation in the rear of the shock without requiring a fine knot spacing over the whole spatial domain. The method also yields accurate results for the width of the shock and its speed of propagation. In table 6 numerical values of the solution calculated using second-order and fourth-order splines with time-independent and time-dependent knots are compared with the analytical solution for v = 0.0075. The results illustrate the increased accuracy in the computed shock profile obtained with a time-dependent knot distribution, and using fourth-order rather than second-order splines. Table 6. Comparisons of the computed and analytical shock profiles X

---

0.550 0.555 0.560 0.565 0.570 0.575 0.580 0.585 0.590 0.595 0.600 0.605 0.610 0.615 0.620 0.625 0.630 0.635 0.640 0.645 0.65

Analytical solution

-_-___ 0.999 0.998 0.997 0.996 0.995 0.993 0.991 0.987 0.982 0.975 0.966 0.952 0.935 0.912 0.881 0.841 0.792 0.731 0.661 0.583 0.500

___

Numerical solution u with time-independent knots

Numerical solution u with time-dependent knots

n=2 h = 25

n=4 h = 25

n=2 h = 24

n=4 h = 24

0.994 0.992 0.991 0.994 0.998 1.002 1.006 1.009 1.013 1.016 1.020 0.973 0.926 0.878 0.832 0.784 0.737 0.690 0.643 0.575 0.509

0.983 0.980 0.979 0.982 0.986 0.991 0.996 0.999 0.999 0.993 0.983 0.965 0.940 0.908 0.870 0.825 0.773 0.716 0.651 0.582 0.508

1.000 0.999 0.999 0.998 0.997 0.997 0.993 0.990 0.986 0.983 0.979 0.954 0.928 0.903 0.877 0.851 0.780 0.709 0.639 0.568 0.497

0.999 0.999 0.998 0.997 0.995 0.993 0.991 0.985 0.979 0.971 0.962 0.951 0.937 0.916 0.886 0.846 0.792 0.728 0.656 0.578 0.499

318

A.M. Davies, Application

of the Galerkin method to the solution of Burgers’ equation

Obviously, in any prolonged time integration the shock would progress a considerable distance in space, and the saving in computer time and storage obtained using a higher resolution region which moved with the shock, rather than a fine knot spacing throughout the whole region, could be appreciable.

Concluding

remarks

It is evident from the results presented previously that the higher order B-splines yield solutions of higher accuracy than the second-order spline for the solution of Burgers’ equation considered here, especially for low values of u. A method for reducing the spurious oscillation in the rear of the shock, and for improving the accuracy of the solution, without using a large basis set has been presented. This method involves moving, as the shock advances in space, a finer resolution region through one of coarser resolution. The results presented here have only involved two regions of different resolution. However, the method of generating the B-splines, evaluating the necessary integrals and moving one region through another are sufficiently general to permit a number of regions of different resolution to be connected. The accuracy of the solution calculated using two regions suggest that a method involving a number of regions of finer resolution embedded within one another such that the region of finest resolution is central and follows the shock could yield particularly accurate results even for low values of V. The results presented here for Burgers’ equation suggest that in solving the shallow water equations with an expansion of B-splines, the steepening of the wave front in shallow water should be resolvable provided a region of higher resolution is used in the area where the wave front becomes particularly steep. However, in order to accurately compute the progress of this wave front, a moving fine resolution region would probably be necessary.

Acknowledgements

The author is indebted to Dr. N.S. Heaps for many useful comments and suggestions in connection with this work. The care taken by Mr. R. Smith in preparing the diagrams is much appreciated. The work described in this paper was funded by a Consortium consisting of the Natural Environment Research Council, the Ministry of Agriculture Fisheries and Food, and the Departments of Energy, Environment and Industry.

Appendix

The method used to determine new matrices C, A, E and n can be illustrated by considering an expansion of second-order B-splines defined upon the set of knots 0.0 (0.1) 0.3 (0.05) 0.5 (0.1) 1.0

319

A.M. Davies, Application of the Galerkin method to the solution of Burgers’equation t--Region

------+

I

Rl

fi

.2

.

3

1

I

MS

M6

MI

MS

Mg

MIO

Ml1

Ml2

Ml

B

B’

M2

B’

S

S

S’

S

T’

T’

T

T’

T’

U

U’

w

u

U’

U’

U

T’

T’

T

T’

T’

s

S’

S’

s

S’

S’

S

S’

S’

S

M4

MS

M4

Rl

R2

M2

M3

M3

Region -------4

-----It-----

MI

0s

f

t-----Region

Fig. 3. The matrix C at time t = 0.

at time t = 0 and subsequently defined upon a new set of knots 0.0 (0.1) 0.4 (0.05) 0.6 (0.1) 1.0 at a later time T when the shock has travelled a distance 2Ah (AX is the distance between adjacent knots in the region of finer resolution). It is convenient to consider the matrix C (neglecting for illustrative purposes knots outside the region 0 < x < 1) to demonstrate the change in the arrangement of integrals within the matrix when the knot positions change. Similar changes occur within the other matrices. Thus for the knot distribution given at time t= 0,C has the form in fig. 3, where

B = J:M,M,dx involves integrals of splines based on boundary knots, B'=JiM,Mjdxwithi= l,i=2, S = J: M,M,dx with Mi defined upon knots all of which are in the coarser resolution region R,, S’=J; MiMjdxwith Mj,Mi adjacent knots in region R 1, U = ld MiMi d.x with Mi defined upon knots all of which are in the fine resolution region R, , U’ = JA MiMjdx with Mi,Mi defined upon adjacent knots in region R,, T = JiM,M, dx with Mi defined upon one knot in region R, and a second knot in region R,, T'= J:M,Mjdx with Mf,Mj defined upon adjacent knots though with one knot in region R, and a second in region R, .

A.M. Davies, Application of the Galerkin method to the solution of Burgers’ equation

320

Region

‘

Region

rR,+-----+

Region

e---------rR2_

I

-RI+-----+

I

I

/ i = .z b d

1

M3

M4

M5

M6

M7

M8

M9

MI

M2

Ml

B

B’

M2

B’

S

S’

S’

s

As’

S’

S

T’

T’

T

T’

T’

U

u’

U’

u

I/’

U’

U

T’

T’

T

T’

T’

S

S’

S’

S

S’

S’

S

n? M3

M4 :

;

MS M6

g 'B

u"

MI

cz Ml3 T ;

J49

Fig. 4. The matrix Cat time t = T.

The matrix C is tridiagonal since the second-order B-splines are only non zero over two knot intervals. At time T, matrix C has the form in fig. 4, corresponding to the movement of the finer knot region along the x axis. Thus no new integrals have to be calculated, although it is necessary to rearrange the integrals within the various matrices. Hence, as the shock front moves through a distance 2Ah, the positions of the knots change, corresponding to a movement of the tine mesh region by a distance 2Ah. The integrals designated U and U’ move within the matrix C; if this tine resolution region reaches the boundary at x = 1, then it would be necessary to calculate new integrals. For the cases considered here the integration period was such that this never occurred.

References [l] A.M. Davies, A numerical investigation of errors arising in applying the Galerkin method to the solution of nonlinear partial differential equations, Comp. Meth. Appl. Mech. Eng. 11 (1977) 341-350. [ 21 M.J. Lighthill, Viscosity effects in sound waves of finite amplitude, in: G.K. Batchelor and R.M. Davies, Eds., Surveys in mechanics (Cambridge Univ. Press, 1956). [3] M.G. Cox, The numerical evaluation ofB-sphnes, J. Inst. Math. Appl. 10 (1972) 1344149.

A.M. Davies, Application of the Galerkin method to the solution of Burgers’equation

321

[4] P.J. Roache, On artificial viscosity, J. Comp. Phys. 10 (1972) 169-184. [Sj P.J. Roache, Computational fluid dynamics (Hermosa Publ., Albuquerque, NM, 1972). [6] P.J. Roache, Recent developments and problem areas in computational fluid dynamics, in: A. Donald and B. Eckman, eds., Lectures notes in physics 461 (Springer, 1974). [7] T.D. Taylor, E. Ndefo and B.S. Masson, A study of numerical methods for solving viscous and inviscid flow problems, J. Comp. Phys. 9 (1972) 99-119. [ 81 O.C. Zienkiewicz and P.N. Godbole. Viscous incompressible flow and special reference to non-Newtonian (plastic) fluids, in: Finite elements in fluid mechanics (Wiley, 1975) ch. 2. [9] G.E. Hill, Grid telescoping in numerical weather prediction, J. Appl. Meteorol. 7 (1968) 29-38. [ 101 N.A. Phillips and J. Shukla, On the strategy of combining coarse and fine grid meshes in numerical weather prediction, J. Appl. Meteorol. 12 (1973) 763-770. [ 111 F.J. Harrison and R.L. Elsberry, A method for incorporating nested finite grids in the solution of systems of geophysical equations, J. Atmos. Sci. 29 (1972) 1235-1245.