Skip to the content of the web site.

Slides for this presentation are available here: 2.Diffusion.pptx. The presentation for this talk is available on youtube.

Submit your homework on UW Learn.

We will now look at the heat-conduction/diffusion equation: given an n-dimensional region R, the solution $u({\bf x}, t)$ defines, in the first case, the heat (or thermal energy) and, in the second, the density of a diffusing material, at the at a point x at a time t. Under ideal conditions, both of these phenomena are modeled by the partial differential equation:

$\frac{1}{\kappa}\frac{\partial}{\partial t}u = \nabla^2 u$

where κ > 0.

This is a partial differential equation in both space and time, and therefore we must solve in one higher dimension. For example, when solving Laplace's equation in one dimension, we only need a solution vector of size U( nx ); but if we now want to solve the heat-conduction/diffusion equation, we require a solution array of size U( nx, nt ) where U( :, ell ) will represent the solution at time $t = t_{\rm initial} + (\ell - 1) \Delta t$. This is shown in Figure 1.

Figure 1. Solving the heat-conduction/diffusion equation in one dimension.

In general, we will know what the initial temperature distribution or density of a diffusing material. That will give us our initial condition U( :, 1 ). What we want to do is to observe how that temperature distribution/density changes over time. Once we find that solution, if we now plot surf( U ), one axis will represent space while the other will represent time. If the boundary conditions are not changing, as time progresses, the solution will approach Laplace's equation.

Suppose, however, we have a problem in two dimensions. In that case, we need a third dimension to store how the solution changes over nt periods of time: U( nx, ny, nt ). Thus, the solution for a particular period of time may be given by U( :, :, ell ) and, as before, this is the approximation at time $t = t_{\rm initial} + (\ell - 1) \Delta t$. This is shown in Figure 2.

Figure 2. Solving the heat-conduction/diffusion equation in two dimension.

Once again, we will know the temperature or density at some initial point and then from that point, we will want to estimate the solution at some time in the future.

To visualize this, Laboratory 6 will demonstrate how we can make a movie in Matlab and so that surf( :, :, ell ) is one frame of the movie. That is, each frame shows the temperature/density and as the frames increase, you will see how the heat moves throughout the object or how the material diffuses through the object.

If we are solving a solution in three spatial dimensions, then we need a four dimensional array U( nx, ny, nz, nt ) where U( :, :, :, ell ) is the solution at time $t = t_{\rm initial} + (\ell - 1) \Delta t$.

It may seem odd to think of a four-dimensional array, but simply think of U( :, :, :, ell ) as a three-dimensional capture of the temperature of, say, a room and then U( :, :, :, ell + 1 ) is the temperature of this room one second (or whatever Δt is) into the future.

Comparison with Laplace

This is different from Laplace's equation: we now have a time component together with a space component. Suppose that we have a solution to Laplace's equation, say that $u({\bf x}, t)$ satisfies

$\nabla^2 u = 0$.

In this case, it follows that

$\frac{1}{\kappa}\frac{\partial}{\partial t}u = 0$

and therefore the solution will not change with time.

However, we may have an initial state which does not satisfy Laplace's equation. In this case, what is the behaviour of the solution? We will look at watching various solutions as they change through time. In this laboratory, we will focus on a one-dimensional space where the solution to Laplace's equation are straight lines (if $u(a) = u_a$ and $u(b) = u_b$ and the function $u(x)$ satisfies Laplace's equation, it follows that the solution is the straight line connecting the two points $(a, u_a)$ and $(b, u_b)$.

In dealing with one or two dimensions, we may write $u(x, t)$ and $u(x, y, t)$, respectively, as opposed to using the vector notation.

Recall that if u is a solution to Laplace's equation, then $\nabla^2 u = 0$ and therefore any solution to Laplace's equation is a steady-state solution to heat-conduction/diffusion equation. If, however, $\nabla^2 u > 0$ at a point, this indicates in one dimension that the function is concave up and therefore the rate of change of u at that point with respect to time will be positive: the solution will move towards a solution of Laplace's equation as time increases. This is shown in Figure 3.

Figure 3. A solution to a function u which is locally concave up at a point x.

In two dimensions, $\nabla^2 u > 0$ implies that the concavities do not cancel. Figure 4 shows one point on a surface where the surface is locally concave up along the x axis and locally concave down along the y-axis, but the concavity is greater along the x-axis which is positive. Such a move is likely to decrease the concavity along the x-axis and increase the negative concavity along the y axis until the concavities are equal and opposite; again, a solution to Laplace's equation.

Figure 4. A solution to a function u where the positive concavity along the x is greater than the negative concavity along the y axis.

Figure 5 shows a the entire surface in a region where the surface is locally concave up along one axis moreso than it is concave down along the other axis.

Figure 5. A function with with locally opposite and different magnitudes of concavity.

Solving the Heat-conduction/Diffusion Equation in One Dimension

In one dimension, we begin with an initial condition: the temperature or density distribution at time $t = t_{\rm initial}$, that is, $u({\bf x}, t_{\rm initial})$ is given. We will always begin with time $t = 0$ and therefore the initial condition will be specified by giving the function $u({\bf x}, 0)$.

We will begin with a one-dimensional space: as with the time-independent boundary-value problem, we will divide space into n_x points where the boundaries are two values $x = a$ and $x = b$ where $h = \frac{b - a}{n_x - 1}$ and we define $x_k = a + (k - 1)h$; however, in this case, we must also take into account the time variable, and thus, we will approximate our solution at intervals of $\Delta t$ where $t_\ell = \ell \Delta t$. We could, for example, solve up to time $t_{\rm final}$ with $n_t$ time intervals and thus set $\Delta t = \frac{t_{\rm final} - t_{\rm initial}}{n_t - 1}$

Thus, we will approximate the solution $u(x_k, t_\ell) \approx u_{x_k, t_\ell}$ where $u(x_k, 0) = u_{k, t_{\rm initial}}$ are our initial conditions and $u(a, t_\ell) = u(x_0, \ell) = u_a$ and $u(b, t_\ell) = u(x_{n_x}, \ell) = u_b$ are our boundary conditions. The unknown points are $u_{k, \ell} \approx u(x_k, t_\ell)$ for $k = 1, 2, \ldots, n_x - 1$ and $\ell = 1, 2, \ldots, $.

This is shown in Figure 6 which shows the boundary values, the initial values, and the unknown values.

Figure 6. The points at which solutions to the heat-conduction/diffusion equation will be approximated.

Note: while mathematically, it is convenient to discuss the range $0, 1, \ldots, n$, in Matlab, the matrices will be indexed from $1, 2, \ldots, n + 1$. It is therefore necessary to compensate appropriately.

Method

Recall that for boundary-value problems, we substituted divided difference formula for each of the components of the equation.

In one dimension, the heat-conduction/diffusion equation simplifies to

$\frac{1}{\kappa}\frac{\partial}{\partial t}u = \frac{\partial^2}{\partial x^2} u$

and thus, because we need to approximate the solution at a point in the future, we could use the approximations

$\frac{\partial^2}{\partial x^2}u(x, t) \approx \frac{u(x + h, t) - 2u(x, t) + u(x - h, t)}{h^2} + O(h^2)$

and

$\frac{\partial}{\partial t}u(x, t) \approx \frac{u(x, t + \Delta t) - u(x, t)}{\Delta t} + O(\Delta t)$.

This would yield

$\frac{1}{\kappa} \frac{u(x, t + \Delta t) - u(x, t)}{\Delta t} = \frac{u(x + h, t) - 2u(x, t) + u(x - h, t)}{h^2}$.

Assuming we have divided the plane up appropriate

$\frac{1}{\kappa} \frac{u_{k, \ell + 1} - u_{k, \ell}}{\Delta t} = \frac{u_{k + 1, \ell} - 2u_{k, \ell} + u_{k - 1, \ell}}{h^2}$

which can be rewritten as

$u_{k, \ell + 1} - u_{k, \ell} = \lambda \left( u_{k + 1, \ell} - 2u_{k, \ell} + u_{k - 1, \ell} \right )$

where $\lambda = \frac{\kappa \Delta t}{h^2}$. This can be shown diagrammatically as sampling the points shown in Figure 7.

Figure 7. The four points sampled in the formula substituting the divided difference approximations of the partial derivatives.

Notice that if $\ell = 0$ then we have

$u_{k, 1} - u_{k, 0} = \lambda \left( u_{k + 1, 0} - 2u_{k, 0} + u_{k - 1, 0} \right )$

This would give an equation in three knowns and only one unknown which we can therefore solve for the one unknown. When $k = 1$, we get the sample shown in Figure 8.

Figure 8. Approximating the points at time $t = t_2$.

We could do this for every x value at this time and thus determine an approximation for each of the $u_{k, 1}$ for $k = 1, 2, 3, ..., n_x - 1$. Once we have all those values, we can use those values to approximate the values for $u_{k, 2}$, then $u_{k, 3}$, and so on.

The formula is

$u_{k, \ell + 1} = u_{k, \ell} + \lambda \left( u_{k + 1, \ell} - 2u_{k, \ell} + u_{k - 1, \ell} \right )$.