|Using EXCEL Spreadsheets to Solve a 1D Wave Equation|
The goal of this tutorial is to describe how to use an EXCEL spreadsheet to compute numerically the solution to the following initial-boundary value problem for the one-dimensional wave equation:
Let us introduce a finite difference mesh xn = nΔx and tk = kΔt, and let the corresponding nodal values be denoted by
Derivatives will be approximated by central differences both in space and time,
By rearranging terms, one obtains
where r= (c Δt/Δx).
The resulting discrete equation involves three distinct time levels: known data from steps k - 1 and k are transferred to step k + 1. In other words, to compute the time evolution of the solution its values at both time tk and an earlier time tk-1 must be known.
This three-level scheme poses some challenges when imposing the initial condition. To calculate the value of un at the first time step k = 1, one needs information from step k = - 1 and k = 0. The initial condition u(x, 0) = f(x) provides information at step k = 0. Any information at an earlier time must be inferred. The initial velocity condition can be used to do that. Imagine a row of false mesh point at time t = - Δt = t-1, then the initial velocity condition can be approximated using central differences as
Assuming that the discrete wave equation also holds at t = 0, one can compute the value of the solution at time t = t1 as
and substitute in the expression for un-1 derived from the initial velocity condition. After re-arranging one obtains
Since un0=f(xn) and g(xn) are known at all nodes, we are now in a position to specify the first two rows of nodes. Let's see how to build an EXCEL spreadsheet that implements this calculation.
In line with the spreadsheet format used in other tutorials, we will be working with a spreadsheet where rows represent time levels, and columns correspond to spatial nodes. In particular, throughout this tutorial we will refer to this spreadsheet.
For this tutorial, consider a spatial interval [0, 1.57], and an increment Δx = 0.04. Note that the value of Δx is stored in cells C1.
As usual, each sample point is calculated as
for n = 1, ... N and x0 = 0.
Each time level is calculated as
for k = 1, ... M and t0 = 0, using a time step Δt = 0.01. Note that the value of Δt is stored in cell E1.
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, we chose
f(x) = exp(-40(x - π /4)2)
The values of f(xn) are entered in row 5.
The values of g(xn) from the initial velocity condition are entered in row 4. For this tutorial, we chose g(x) = 0. These values will be used in computing the first time step, as shown in Step 3.
The first time step is computed and stored in row 7. In each cell of row 7 that represents a node, an EXCEL formula is coded to compute the expression
which is obtained by re-arranging the terms in the expression for un1 derived in the introduction of this tutorial. Note that the value of r2 is stored in cell I1.
For example, u11 is computed as
and its value is stored in cell E7. The EXCEL formula coded in cell E7 must call the following cells:
The resulting formula is E6+$E$1*E5+0.5*$I$1*(D6-2*E6+F6). Note that Δt and r2 are constants, thus the cells containing their values are called by absolute referencing using the dollar signs in their cell address.
To compute the solution at the adjacent node, u21, it is sufficient to copy and paste the formula used in cell E7 into the cell to the immediate right, that is, cell F7. Copying and pasting the same formula along row 7 allows you to compute the solution at each node at time t = t1. Note that the value of the solution on the last node is given by the boundary condition of the problem.
To compute the solution at later times the recursion formula is implemented
This expression is obtained by re-arranging the terms in the equation for unk+1 derived in the introduction of this tutorial.
For example, u12 is computed as
and its value is stored in cell E8. The EXCEL formula coded in cell E8 must call the following cells:
The resulting formula is 2*E7-E6+$I$1*(D7-2*E7+F7). Again, since r2 is a constant, the cell containing its value is called by absolute referencing using the dollar signs in its cell address.
To compute the solution at the adjacent node, u22, simply copy and paste the formula used in cell E8 into the cell to the immediate right, that is, cell F8. Then, repeat the same procedure to fill in the cells in row 8 to compute the solution at each node at time t = t2, and similarly, into each row below it to find the time evolution of the initial wave.
The Edit|Fill|Fill Right and Edit|Fill|Fill Down features will allow you to fill in all necessary cells and calculate the solution quickly.
The wave propagation can be easily observed if one plots the solution computed in the spreadsheet at different time levels. Instead of creating several overlapping graph of u at different t, in EXCEL it is possible to generate an animation that runs through the values of the solution previously calculated in the spreadsheet. This animation is implemented in the spreadsheet described in this tutorial.
First of all, observe the plot shown in Fig 1. It depicts the solution u computed in the spreadsheet over the entire spacial domain at a certain time interval.
A slider has been added in cell X1. By clicking the bar on the slider in cell X1, the plot in Fig 1 changes. As the bar on the slider is moved from left to right, the plots show the time evolution of the wave signal.
Let's see how the animation works. The graph in Fig 1 plots the values stored in row 3, from cell D3 through AR3. To see how these values were generated, consider the formula coded in the cells along row 3. The formulas make use of the function OFFSET( ).
In its simplest version, the function OFFSET returns a reference to a cell that is offset by a certain number of rows and columns from a given cell. Specifically, OFFSET(reference, rows, cols) returns a reference as follows: reference is the cell address from which you want to base the offset; rows is the number of rows, up or down, to be added or removed to or from the cell address of reference; cols is the number of columns, to the left or right, to be added or removed to or from the cell address of reference. For example, if OFFSET(C2,1,2) is coded in a given cell, the value shown in that cell is the value stored in a cell one row below and two column to the right of cell C2, i.e., cell E3.
In the spreadsheet described in this tutorial, the OFFSET function is used in the following manner. The value shown, for example, in cell D3 is generated by OFFSET(D8,$M$1,0), that is, by offsetting the reference D8 by a number of rows equal to the value shown in cell M1, without changing column. When the bar on the slider is all the way to the left, the value in cell M1 is 0. Thus, the OFFSET function coded in cell D3 returns the value in cell D8 offset by 0 rows and 0 columns. In other words, it returns the value in D8, which is u02. Similarly, the value shown in cell E3 is simply the value read from cell E8, which is u12, and so forth. So, when the value in cell M1 is 0 the plot in Fig1 shows the solution computed at time t = t2.
Now observe what happens to cell M1 when the bar on the slider is moved slightly to the left. The value in cell M1 increases to 1, causing the values in row 3 to change. Now the cells in row 3 return the values stored in a row that is offset from row 8 by 1, i.e., row 9. In other words, now row 3 shows the solution computed at time t = t3. By moving the slider progressively to the left, the cells in row 3 point to subsequent rows and the solution computed at subsequent time steps is plotted in sequence in Fig 1.