UBC Home Page
UBC Home Page - UBC Home Page UBC Home Page - - -
-

 

 
Using EXCEL Spreadsheets to Solve a 1D Heat Equation
 

The goal of this tutorial is to create an EXCEL spreadsheet that calculates the numerical solution to the following initial-boundary value problem for the one-dimensional heat equation:

The basic idea of the numerical approach to solving differential equations is to replace the derivatives in the heat equation by difference quotients and consider the relationships between u at (x,t) and its neighbours a distance Δx apart and at a time Δt later. In particular, in this tutorial the following expressions will be used:

By rewriting the heat equation in its discretized form using the expressions above and rearranging terms, one obtains

Hence, given the values of u at three adjacent points x-Δx, x, and x+Δx at a time t, one can calculate an approximated value of u at x at a later time t+Δt.

Using EXCEL spreadsheets allows you to perform these calculations repeatedly and effortlessly. The following are the necessary steps to set up a spreadsheet to calculate the solution to the initial-boundary value problem shown above.

  1. How to Tabulate the Spatial Interval
  2. How to Tabulate the Time Interval
  3. How to Include the Initial Condition
  4. How to Include the Boundary Conditions
  5. How to Compute the Numerical Solution
  6. How to Plot the Solution
  7. How to Implement Derivative Boundary Conditions



Step 1: How to Tabulate the Spatial Interval

The first step is to tabulate the sample points in the spatial interval. Subdivide the interval [0, 1] into N+1 equally-spaced sample points a distance Δx apart.

In a blank EXCEL spreadsheet, enter the sample points along a row, as shown below. For this tutorial, let Δx = 0.05 and N = 19. Note that the value of Δx is shown in cell C2.

Each sample point is calculated as

for n = 1, ... N and x0 = 0.

For example, in the spreadsheet below the first sample point in the interval is x1 = x0 + 0.05. The value of x1 is shown in cell E5: it is computed by addying the fixed space increment (cell C2) to the value contained in the cell to the immediate left of E5 (cell D5). The expression in the formula box in the image below shows the specific EXCEL formula used to compute the value in cell E5.



Step 2: How to Tabulate the Time Interval

The second step is to tabulate the time interval. Subdivide the interval [0,T] into M+1 equal time levels Δt long.

Enter the time levels in a column in the EXCEL spreadsheet, as shown below. For this tutorial, let Δt = 0.004 and M = 20. Note that the value of Δt is shown in cell F2.

Each time level is calculated as

for k = 1, ... M and t0 = 0.

For example, in the spreadsheet below the second time level is t2=t1 + 0.004. The value of t2 is shown in cell B9: it is computed by adding the fixed time increment (cell F2) to the value contained in the cell immediately above B9 (cell B8). The expression in the formula box in the image below shows the specific EXCEL formula used to compute the value in cell B9.



Step 3: How to Include the Initial Condition

The next step is to assign the initial condition u(x,0) = f(x). Note that the specific form of f(x) depends on the characteristic of the problem under consideration. For this tutorial, choose

f(x) = 2x for 0 < x < 0.5 and f(x) = 2-2x for 0.5 < x < 1.

For clarity, in the spreadsheet below space-time points have been highlighted in blue and the corresponding values of the solution at t = 0 in green. In other words, the values of the initial condition are shown in row 7, in cells D7 through X7. For example, the value of u(0.1, 0), that is, the value of u at x2 at t0, is shown in cell F7.



Step 4: How to Include the Boundary Conditions

The initial-boundary value problem discussed in this tutorial has two boundary conditions:

u(0, t) = 0 and u(1, t) = 0.

In the spreadsheet shown below, column D, from cells D7 through D27, contains the values corresponding to the first boundary condition u(0, t) = 0, that is, it shows the constant value of u at x0. Column X, from cells X7 through X27, instead, contains the second boundary condition u(1, t) = 0, that is, it shows the constant value of u at xN. Note that at t0, u(1, 0) = 0 (cell X7) in agreement with the initial condition.



Step 5: How to Compute the Numerical Solution

To calculate the value of u at each space-time sample point (xn , tk ), the following algorithm is used

for n = 1, 2, ... N. Thus, unk+1, the approximated value of u at point xn at a time tk+1, can be calculated given the values of u at three adjacent points, xn-1, xn, xn+1, at a time tk.

For example, the value of u at x1 at time t1 can be approximated as u(x1, t1) ~ u11, and calculated using the expression above with n = 1 and k = 0,

where

Hence, given the value of u at x0, x1, x2 at time t0, one can find an approximated value of u at x1 one time step later. The spreadsheet below shows how to perform this calculation:

Note that cell O2 is called by absolute referencing ($O$2 instead of O2). This is because the formula used in cell E8 will be copied and pasted in other cells, as explained later on. It is worth noticing that the calculation performed in cell E8 uses the numerical values stored in three cells in the row just above it, specifically in cell E7, the cell immediately above E8, and cells D7 and F7, the cells to the left and right of E7, respectively. This way of storing information in rows and columns allows you to organize your computation in a structured manner. As you fill in the spreadsheet to calculate u over the entire domain, information will always be retrieved from cells one row above the row you are working on, in the same column and in the columns to the immediate left and right.

One could perform a similar calculation to find u21, the (approximated) value of u at a sample point a distance Δx to the right of x1 and at the same time level as u11. The result of the computation would be stored in the same row (or time level) as u11 (row 8), but one column to the right (column F), to be consistent with the way sample space-time points are tabulated in the spreadsheet. Furthermore, the computation of u21 would involve the values of u at x1, x2, and x3 at time t0, which are stored in cells E7, F7, and G7. In other words, the computation of u21 would involve information in three cells located to the immediate right with respect to the cells used to compute u11. Therefore, instead of typing a new formula for u21 in cell F8, one could simply Copy and Paste the formula from cell E8 into cell F8, EXCEL will then automatically update the formula. In fact, when a formula is copied from one cell to another, EXCEL does not create an exact copy of the formula, instead it changes cell addresses relative to the row and column they are moved to. When the formula from cell E8 is pasted one cell to the right (into F8), all column labels in the formula, with the exception of cell O2, will be shifted to the right, and the correct values of u will be accessed. The cell address for O2 will not change because the column and row labels are preceded by the "$" sign (this is called "absolute referencing" of cells), allowing for the constant value of K to be read always from cell O2.

The same method can be used to calculate values of u at later times. Instead of typing a new formula in subsequent rows, simply copy and paste the formula from the cell corresponding to the previous time level, i.e. the cell in the row above.

Keep in mind that in order to compute the time evolution of u at a particular point in the space interval, its earlier value at that point and at adjacent points must be known. Therefore, before calculating how u changes with time at a given point, it is convenient to find earlier values of u at all points in the spatial interval. In other words, when using the spreadsheet to calculate u, first fill in all cells in one row, and then move to the row below, so that the necessary values of u are always available. The Edit|Fill|Fill Right and Edit|Fill|Fill Down features will allow you to fill in all cells and calculate the solution quickly.



Step 6: How to Plot the Solution

The time evolution of u can be easily observed if one plots the solution computed in the spreadsheet at different time levels, and compare it with the initial condition. Here is a sequence of steps to make a plot of the solution at time t = 0 and at time t = 0.08 (corresponding to M = 20).

Your plot may look like the image below. By increasing M and implementing the computation on more time level (i.e., by filling more rows in the spreadsheet), you can investigate the evolution of u at even later time.



Step 7: How to Implement Derivative Boundary Conditions

Suppose the boundary conditions to the initial-boundary value problem discussed at the beginning of this tutorial are changed to

Unlike the previous case, in which the boundary conditions provided the value of the solution on both boundaries of the domain, now the value of the solution is known only at x = 0. Thus, u(1, t) must be calculated using the same difference approximation for the heat equation used in Step 5, but now n = N, that is,

By doing so, one finds that the value of u at xN depends on the value of u at xN+1. Recall that xN+1 = xN + Δx. Since xN = 1, then, xN+1 is outside the domain and must be defined. In practical terms, this means adding a column in the spreadsheet used to calculate the solution. Alternatively, one could seek a new expression for uNk that does not depend on the value of u at xN+1, avoiding therefore the additional column in the spreadsheet. In both cases, it will be essential to use the information provided by the derivative boundary condition, as described below.

Consider a central difference approximation of the derivative boundary condition,

It follows that

Since xN+1 = xN + Δx and xN-1 = xN - Δx, it follows that uN+1k = uN-1k. This information can be used to calculate uNk+1 in two ways: