Maximum norm contractivity and accuracy of the Crank-Nicolson

In document Numerical Treatment of Linear Parabolic Problems (Pldal 66-69)

2.3 Discrete analogs of the qualitative properties - reliable discrete models

2.4.3 Maximum norm contractivity and accuracy of the Crank-Nicolson

Applying the Crank-Nicolson discretization to the heat conduction problem, our aim is to get reliable numerical model. This means that we have to guarantee the maximum norm contractivity, too. However, as we have already seen, this requires a bound for the discretization step-sizes. In this section we consider this condition and analyze its effect to the numerical accuracy.

For one-step finite difference methods which are based on some rational approximation of the exponential function, Spijker has shown in [129] that there is an order barrier:

only methods with first order of accuracy can be contractive in the maximum norm for all q > 0. (Such a method is the backward Euler method with the stability function rBE(qK0) = (I−qK0)−1.) For higher order methods it is necessary to restrict the choice ofq(from above) in order to preserve the maximum norm contractivity. As we have seen in the previous section, in Theorem 2.4.8, for the second order Crank-Nicolson scheme, when rCN(qK0) = (I+q2K0)(Iq2K0)−1 is applied to the one-dimensional heat equation, then the sharp restriction is q≤1.5. Therefore, the use of the Crank-Nicolson scheme requires the choice of a very small time stepτ in the case of a small space discretization parameter hif we want to preserve the maximum norm contractivity. One of the main problems when using the Crank-Nicolson scheme is that it preserves the maximum norm contractivity only for q 1.5. This means that if we have a fine mesh for the space variable, we must choose the time step τ 1.5h2 in order to have maximum norm contractivity. For large values ofN (whereN denotes the number of the partition in space, i.e., h= 1/N) it would mean very small (even useless)τ and it requires considerable computational efforts.

Moreover, this also results in an essential loss of accuracy: the computational error, as the result of the large number of iteration steps, might cause an essential loss in the accuracy, i.e., the scheme may loose its second order accuracy.

We illustrate this problem on the numerical solution of the initial-boundary value problem

∂u

∂t −∂2u

∂x2 = 0, t >0, x(0,1), u(0, t) =u(1, t) = 0, t≥0,

u(x,0) =u0(x), x∈[0,1].

(2.4.24)

We demonstrate the behaviour of the numerical solution on two examples, namely, one with a smooth initial function and one with a non-smooth initial function. In the examples we compare the errors at the same fixed time level T = 1.

Example 2.4.9 The first model problem is (2.4.24) with the smooth initial function u0(x) = sin(πx). The exact solution is u(x, t) = sin(πx) exp(−π2t).

Example 2.4.10 The second model problem is (2.4.24) with the non-smooth initial function

u0(x) =

½ 1 if x∈[0.25,0.75]

0 otherwise .

q 100 50 40 20 10 8 CN-error 6.35(−5) 4.14(−5) 3.10(−5) 8.79(−6) 1.53(−6) 6.04(−7) BE-error 6.90(−3) 1.60(−3) 1.00(−3) 2.81(−4) 9.83(−5) 7.29(−5)

q 7 6 5 4 2 1.5

CN-error 2.22(−7) 1.147(−7) 4.05(−7) 6.40(−7) 9.54(−7) 9.89(−7) BE-error 6.27(−5) 4.86(−5) 4.07(−5) 3.15(−5) 1.50(−5) 1.20(−5)

q 1 0.4 0.1 0.05 0.01 0.005

CN-error 1.03(−6) 1.06E(−6) 1.06(−6) 1.06(−6) 1.06(−6) 1.06E(−6) BE-error 7.75E(−6) 3.67(−6) 1.70(−6) 1.38(−6) 1.12(−6) 1.09(−6) Table 2.4.2: The maximum norm error for h= 0.05 for the CN and BE methods.

q 4000 2000 1500 1000 500

CN-error 3.16(−5) 9.70(−6) 5.12(−6) 2.54(−6) 6.35(−7) BE-error 9.91(−4) 2.76(−4) 1.58(−4) 9.59(−5) 3.91(−5)

q 400 200 100 40 20

CN-error 4.03(−7) 9.30(−8) 1.54(−8) 6.35(−9) 9.46(−9) BE-error 3.00(−5) 1.38(−5) 6.59(−6) 2.58(−6) 1.28E(−6)

q 10 4 2 1.5 1

CN-error 1.02(−8) 1.05(−8) 1.05(−8) 1.05(−8) 1.05(−8) BE-error 6.43(−7) 2.63(−7) 1.37(−7) 1.05(−7) 7.35(−8)

Table 2.4.3: The maximum norm error forh= 0.005 for the CN and BE methods.

The exact solution is u(x, t) = 2

π X

m=1

1 m

µ

cos

4 cos3mπ 4

sin(mπx) exp(−m2π2t).

First we consider the maximum norm error of the Crank-Nicolson (CN) and the backward Euler (BE) methods, applied to the numerical solution of Example 2.4.9, on different meshes. Table 2.4.2 and 2.4.3 show that by refining the mesh under the maximum norm contractivity condition, the Crank-Nicolson scheme loses its higher accuracy with respect to the backward Euler method and in limit they result in the same accuracy.

The following Tables 2.4.4 - 2.4.6 serve to demonstrate the behavior of the maximum norm error of the Crank-Nicolson scheme with further different discretization parameters.

These results show that the optimal accuracy of the Crank-Nicolson scheme is attained at some value qopt = qopt(h) which is greater than 1.5. Moreover, by decreasing h (that

q 10 5 2 1.5 1

error 2.93(−5) 5.93(−6) 2.62(−6) 3.23(−6) 3.92(−6)

q 0.5 0.1 0.05 0.01 0.005

error 4.25(−6) 4.35(−6) 4.36(−6) 4.36(−6) 4.36(−6)

Table 2.4.4: The maximum norm error of the Crank-Nicolson scheme for h= 0.1.

q 1000 500 250 160 80 40 20 error 0.245 0.014 9.54(−5) 3.14(−5) 9.48(−6) 2.29(−6) 3.84(−7)

q 16 4.8 1.6 1.5 1 0.8 0.16

error 1.52(−7) 1.57(−7) 2.59(−7) 2.59(−7) 2.61(−7) 1.85(−7) 1.86(−7) Table 2.4.5: The maximum norm error of the Crank-Nicolson scheme for h= 0.025.

q 1000 500 100 50 40 30

error 3.16E(−5) 9.67(−6) 3.72(−7) 6.16(−8) 2.43(−8) 4.61(−9)

q 20 10 5 1.5 1 0.5

error 2.54(−8) 3.79(−8) 4.10(−8) 2.96(−8) 2.97(−8) 2.97(−8) Table 2.4.6: The maximum norm error of the Crank-Nicolson scheme for h= 0.01.

is, by refining the mesh), the values qopt are increasing. Table 2.4.7 shows the loss in the accuracy. The fifth column in this table shows how much more CPU-time is used to obtain the less accurate result (with q = 1.5). The accuracy with the choice q = 1.5 is attained with some qbig >1.5, too. The approximate values of these parameters and the corresponding CPU ratios are included in the last two columns.

Let us compare the above numerical methods on the non-smooth problem, i.e., on the Example 2.4.10 by using the Crank-Nicolson and the backward Euler methods.

Table 2.4.8 summarizes the errors in maximum norm for the Crank-Nicolson and back-ward Euler methods. The behaviour of the Crank-Nicolson scheme is similar as for the smooth initial function. However, the smoothing property of the backward Euler method is considerable. We remark that the same conclusions can be made for the other choices of h.

The Crank-Nicolson scheme has a local errorO(τ2+h2). Therefore, the optimal accuracy (i.e., second order), is achieved when τ 1/N, (i.e., τ h) and not at τ 1/N2, which is required by the contractivity condition. (The latter implies that for a fixedh the order of the error is defined only by τ.)

h qopt error error forq = 1.5 CPU ratio qbig CPU ratio

0.1 2 2.62(−6) 3.23(−6) 1.266 4.6 3.07

0.05 4 6.40(−7) 9.87(−7) 2.71 8.7 5.80

0.025 16 1.52(−7) 2.59(−7) 10.81 18 12.00

0.01 30 4.61(−9) 4.19(−8) 20.5 45 30.00

0.005 40 6.35(−9) 1.05(−8) 27.6 94 62.67

0.004 62.5 2.57(−9) 6.71(−9) 42.9 112.5 75.00 Table 2.4.7: Comparison of the accuracy and consumed CPU time.

q 4000 2000 400 200 100 CN-error 0.4765 0.4438 0.2397 7.86(−2) 2.30(−3) BE-error 8.99(−4) 2.51(−4) 2.76(−5) 1.29(−5) 6.35(−6)

q 75 60 50 45 40

CN-error 8.18(−5) 1.48(−6) 3.68(−7) 3.69(−7) 3.71(−7) BE-error 4.83(−6) 3.88(−6) 3.30(−6) 2.99(−6) 2.70(−6)

q 20 10 5 1.5 1

CN-error 3.73(−7) 3.74(−7) 3.74(−7) 3.74(−7) 3.74(−7) BE-error 1.53(−6) 9.48(−7) 6.61(−7) 4.60(−7) 4.31(−7)

Table 2.4.8: The maximum norm error for h = 0.005 for the CN and BE methods for a non-smooth initial function.

2.4.4 Maximum norm contractivity for the modified Crank-Nicolson

In document Numerical Treatment of Linear Parabolic Problems (Pldal 66-69)