Step-by-step look at a numerical (in)stability

I warned you last time that simulating time-dependent heat conduction might be tricky, due to numerical instabilities. What did I mean by that? Let's find out today by following the operations of a computer which is implementing the equations we derived last time, step by step.

The situation is this: a metal rod of length L = 2 m is initial at temperture T = 0. Then, at time t = 0 sec we attach to the left end a source of heat at a fixed temperature TL = 100 K, and to the right end of the rod a source of heat at a fixed temperature of TR = 10 K. The coefficient of thermal diffusivity of this metal is k = 10-4 m^2/s. What happens next?

Well, that depends on exactly how we choose to simulate the physical situation. Suppose we try this combination of parameters:

• break the rod into sections which are dx = 0.1 m wide
• use a timestep of dt = 100 sec

Let's focus on the left-hand end of the rod, connected to the source at fixed temperature TL = 100 K.

```

d2T        T(i+1) + T(i-1) - 2T(i)                          d2T
Time           Temperature                      ---- =    -------------------------                     k  -----  dt
dx2                dx*dx                                    dx2

100    0    0    0    0    0    0       X   10000    0     0     0     0     0       X    100    0    0    0    0    0
0     ____ ____ ____ ____ ____ ____ ____     ____ _____ _____ _____ _____ _____ _____     ____ ____ ____ ____ ____ ____ ____

x=0   0.1  0.2  0.3  0.4  0.5  0.6     x=0   0.1   0.2   0.3   0.4   0.5   0.6      x=0   0.1  0.2  0.3  0.4  0.5  0.6

100     ____ ____ ____ ____ ____ ____ ____     ____ _____ _____ _____ _____ _____ _____     ____ ____ ____ ____ ____ ____ ____

200     ____ ____ ____ ____ ____ ____ ____     ____ _____ _____ _____ _____ _____ _____     ____ ____ ____ ____ ____ ____ ____

300     ____ ____ ____ ____ ____ ____ ____     ____ _____ _____ _____ _____ _____ _____     ____ ____ ____ ____ ____ ____ ____

400     ____ ____ ____ ____ ____ ____ ____     ____ _____ _____ _____ _____ _____ _____     ____ ____ ____ ____ ____ ____ ____

```

That doesn't look very good, does it? Let's watch the evolution of this simulation graphically:

Hey! The temperature at some locations has exceeded our initial range of 0 - 100 Kelvin. In order to see all the values in our numerical array, we need to expand the range. All the following graphs will use the same, much larger range.

If we let this simulation run for one hour, we end up with baloney; exceedingly bogus baloney, at that.

Look at your worksheet. Can you explain what's going on in this simulation to cause the instability?

```

```

How can we avoid that kind of behavior? The answer is to decrease the size of the "steps" in our simulation. We can decrease the steps in time, or the steps in position, or both.

For example, let's cut in half the size of the steps in time:

• break the rod into sections which are dx = 0.1 m wide
• use a timestep of dt = 50 sec

Again, let's focus on the left-hand end of the rod, connected to the source at fixed temperature TL = 100 K.

```

d2T        T(i+1) + T(i-1) - 2T(i)                          d2T
Time           Temperature                      ---- =    -------------------------                     k  -----  dt
dx2                dx*dx                                    dx2

100    0    0    0    0    0    0       X   10000    0     0     0     0     0       X     50    0    0    0    0    0
0     ____ ____ ____ ____ ____ ____ ____     ____ _____ _____ _____ _____ _____ _____     ____ ____ ____ ____ ____ ____ ____

x=0   0.1  0.2  0.3  0.4  0.5  0.6     x=0   0.1   0.2   0.3   0.4   0.5   0.6      x=0   0.1  0.2  0.3  0.4  0.5  0.6

100                                     X                                            X
50     ____ ____ ____ ____ ____ ____ ____     ____ _____ _____ _____ _____ _____ _____     ____ ____ ____ ____ ____ ____ ____

100                                     X                                            X
100     ____ ____ ____ ____ ____ ____ ____     ____ _____ _____ _____ _____ _____ _____     ____ ____ ____ ____ ____ ____ ____

100                                     X                                            X
150     ____ ____ ____ ____ ____ ____ ____     ____ _____ _____ _____ _____ _____ _____     ____ ____ ____ ____ ____ ____ ____

100                                     X                                            X
200     ____ ____ ____ ____ ____ ____ ____     ____ _____ _____ _____ _____ _____ _____     ____ ____ ____ ____ ____ ____ ____

100                                     X                                            X
250     ____ ____ ____ ____ ____ ____ ____     ____ _____ _____ _____ _____ _____ _____     ____ ____ ____ ____ ____ ____ ____

100                                     X                                            X
300     ____ ____ ____ ____ ____ ____ ____     ____ _____ _____ _____ _____ _____ _____     ____ ____ ____ ____ ____ ____ ____

```

That looks much better, doesn't it? Let's watch the evolution of this simulation graphically:

Notice that there's no need to change the scale in this case: all the simulated temperatures fall within the initial range of 0 - 100 Kelvin.

If we let this simulation run for one hour, we end up with a rather nice result:

Look at your worksheet. Can you explain what's going on in this simulation to AVOID the instability?

```

```

It turns out that one condition for stability in this time-dependent solution is pretty simple: the combination of coefficient of diffusivity k, time-step dt and position-step dx must satisfy

This is necessary, but not sufficient: it's possible that under some circumstances, even if the condition above is true, the numerical solution to a problem will still oscillate endlessly (but at least the amplitude of the oscillations won't grow). In order to be guaranteed that the solution will converge to a proper solution,

Using even smaller fractions in this condition may be beneficial, up to a point -- when the fraction is (1/6), truncation errors tend to be minimized.

Note that these conditions hold only for the algorithm we've examined: a simple, Euler-like method which uses the current second derivative in position to compute the current first derivative in time, and assumes that the first derivative is constant over the entire timestep. There are more complex algorithms which are less likely to become unstable, which means that one might use larger time-steps and/or position-steps to reach the goal more quickly.