Applications - A Guide to MATLAB for Beginners and Experienced Users - page 189

This looks much better! As time increases, the temperature distribution
seems to approacha linear function of x . Indeed u ( x , t ) = 20 + x is
the limiting “steady state” for this problem; it satisfies the boundary
conditions and it yields 0 on bothsides of the partial differential
equation.
Generally speaking, it is best to understand some of the theory of partial
differential equations before attempting a numerical solution like we have
done here. However, for this particular case at least, the simple rule of
thumb of keeping the coefficients of the iteration positive yields realistic
results. A theoretical examination of the stability of this finite difference
scheme for the one-dimensional heat equation shows that indeed any value
of s between 0 and 0.5 will work, and it suggests that the best value of t to
use for a given x is the one that makes s = 0 . 25. (See Partial Differential
Equations: An Introduction, by Walter A. Strauss, John Wiley and Sons,
1992.) Notice that while we can get more accurate results in this case by
reducing x , if we reduce it by a factor of 10 we must reduce t by a factor of
100 to compensate, making the computation take 1000 times as long and use
1000 times the memory!

The Case of Variable Conductivity
Earlier we mentioned that the problem we solved numerically could also be
solved analytically. The value of the numerical method is that it can be
applied to similar partial differential equations for which an exact solution
is not possible or at least not known. For example, consider the
one-dimensional heat equation with a variable coefficient , representing an
inhomogeneous material with varying thermal conductivity k ( x ),

k ( x ) ∂ u
∂ x

2 u

∂ u

∂ x 2 + k ( x ) ∂ u ∂ t = ∂

= k ( x ) ∂

∂ x . ∂ x

For the first derivatives on the right-hand side, we use a symmetric finite
difference approximation, so that our discrete approximation to the partial
differential equations becomes

u j + 1 − u j − 1
2 x ,
where k j = k ( a + j x ). Then the time iteration for this method is

u n + 1

= k j u j + 1 − 2 u j + u j − 1
x 2

− u j
t

k j + 1 − k j − 1
2 x

j

+

= s j u j + 1 + u j − 1 + (1 − 2 s j ) u j + 0 . 25 ( s j + 1 − s j − 1 ) u j + 1 − u j − 1 , u n + 1

j