Download (direct link):
A Stiff Problem
Consider the initial value problem
y = -100y + 100t + 1, y(0) = 1.
Find numerical approximations to the solution for 0 < t < 1 using the Euler, backward Euler, and Runge-Kutta methods. Compare the numerical results with the exact solution.
Since the differential equation is linear, it is easy to solve, and the solution of the initial value problem (16) is
y = ?(t) = e-100t + t.
Some values of the solution fi(t), correct to six decimal places, are given in the second column of Table 8.5.3, and a graph of the solution is shown in Figure 8.5.2. There is a thin layer (sometimes called a boundary layer) to the right of t = 0 in which the exponential term is significant and the solution varies rapidly. Once past this layer, however, $(t) = t and the graph of the solution is essentially a straight line. The width of the boundary layer is somewhat arbitrary, but it is certainly small. At t = 0.1, for example, exp(-100t) = 0.000045.
FIGURE 8.5.2 The solution of the initial value problem (16).
8.5 More on Errors; Stability
If we plan to approximate the solution (17) numerically, we might intuitively expect that a small step size will be needed only in the boundary layer. To make this expectation a bit more precise recall from Section 8.1 that the local truncation errors for the Euler and backward Euler methods are proportional to 0" (t). For this problem 0"(t) = 104 e—100t, which varies from a value of 104 at t = 0 to nearly zero for t > 0.2. Thus a very small step size is needed for accuracy near t = 0, but a much larger step size is adequate once t is a little larger.
On the other hand, the stability analysis in Eqs. (9) through (15) also applies to this problem. Since r = —100 for Eq. (16), it follows that for stability we need h < 0.02 for the Euler method, but there is no corresponding restriction for the backward Euler method.
Some results using the Euler method are shown in columns 3 and 4 of Table 8.5.3. The values obtained for h = 0.025 are worthless due to instability, while those for h = 0.01666... are reasonably accurate for t > 0.2. However, comparable accuracy for this range of t is obtained for h = 0.1 using the backward Euler method, as shown by the results in column 7 of the table.
The situation is not improved by using instead of the Euler method a more accurate one, such as Runge-Kutta. For this problem the Runge-Kutta method is unstable for h = 0.033 ... but stable for h = 0.025, as shown by the results in columns 5 and 6 in Table 8.5.3.
The results given in the table for t = 0.05 and for t = 0.1 show that in the boundary layer a smaller step size is needed to obtain an accurate approximation. You are invited to explore this matter further in Problem 3.
TABLE 8.5.3 Numerical Approximations to the Solution of the Initial Value Problem
y = —100y + 100t + 1, y(0) = 1
t Exact Euler 0.025 Euler 0.0166... Runge-Kutta 0.0333 ... Runge-Kutta 0.025 Backward Euler 0.1
0.0 1.000000 1 . 000000 1 . 000000 1.000000 1 . 000000 1 . 000000
0.05 0.056738 2.300000 -0.246296 0.470471
0.1 0.100045 5.162500 0.187792 10.6527 0.276796 0.190909
0.2 0.200000 25.8289 0.207707 111.559 0.231257 0.208264
0.4 0.400000 657.241 0.400059 1.24 x 104 0.400977 0.400068
0.6 0.600000 1.68 x 104 0.600000 1.38 x 106 0.600031 0.600001
0.8 0.800000 4.31 x 105 0.800000 1.54 x 108 0.800001 0.800000
1.0 1.000000 1 . 11 x 107 1 . 000000 1.71 x 1010 1 . 000000 1 . 000000
As a final example, consider the problem of determining two linearly independent solutions of the second order linear equation
f - 10n2y = 0 (18)
for t > 0. The generalization of numerical techniques for the first order equations to higher order equations or to systems of equations is discussed in Section 8.6, but that is not needed for the present discussion. Two linearly independent solutions of Eq. (18) are 0x(t) = coshV^??n t and 02(t) = sinhyT?n t. The first solution, 0x(t) = coshV^0n t, is generated by the initial conditions 0X (0) = 1,0'1 (0) = 0; the second solution, 02 (t) =
Chapter 8. Numerical Methods
sinh vT0nt, is generated by the initial conditions 02(0) = 0, 02(0) = V70n. While analytically we can tell the difference between cosh VT0n t and sinh VT0n t, for large t we have cosh VT0n t ~ e710^ 72andsinhVT0n t ~ e710^ /2; numerically these two functions look exactly the same if only a fixed number of digits are retained. For example, correct to eight significant figures, we find that for t = 1
sinhV70n = coshV70n = 10,315.894.
If the calculations are carried out on a machine that carries only eight digits, the two solutions 01 and 02 are identical at t = 1 and indeed for all t > 1. Thus, even though the solutions are linearly independent, their numerical tabulation would show that they are the same because we can retain only a finite number of digits. This phenomenon is called numerical dependence.
For the present problem we can partially circumvent this difficulty by computing, instead of sinhvT0n t and coshvT0n t, the linearly independent solutions 03(t) = e^/Tont and 04(t) = e~^nt corresponding to the initial conditions 03(0) = 1, 03(0) = VT0n and 04(0) = 1, 04 (0) = — VT0n, respectively. The solution 03 grows exponentially while 04 decays exponentially. Even so, we encounter difficulty in calculating 04 correctly on a large interval. The reason is that at each step of the calculation for 04 we introduce truncation and round-off errors. Thus at any point tn the data to be used in going to the next point are not precisely the values of 04(tn) and 04 (tn). The solution of the initial value problem with these data at tn involves not only e—710n t but also e710^t. Because the error in the data at tn is small, the latter function appears with