Recruiting Solver to Iteratively Solve Finite Difference Equations

Problem

You've set up a finite difference solution to an elliptic boundary value problem, but instead of directly solving a system of algebraic equations, you want to use an iterative approach.

Solution

Set up your finite difference equations in Excel and use Solver to iteratively find a solution.

Discussion

In Recipe 12.1, I explained how you can use Excel's matrix functions to directly solve systems of finite difference equations. For small systems, direct matrix operations are viable; however, for large systems the use of matrix functions in Excel is cumbersome. Other standard direct methods include Gauss elimination and Fourier's method, among others. These methods could be programmed in Excel; however, there's another class of methods known as iterative methods that are more attractive.

Iterative solutions to finite differences equations abound. Various standard methods such as Jacobi's method , Gauss-Seidel iteration (also called Liebmann's method ), and the alternating direction implicit (ADI) method are commonly used in lieu of direct methods. Iterative methods have an advantage over direct methods in that they typically require less memory and are relatively easy to program. Moreover, when using Excel to help solve your problem, you already have access to a powerful iterative tool, Solver. This means that you can use Solver to find a solution to your finite difference equations without any programming beyond setting up a spreadsheet.

By way of example, consider the steady state temperature distribution in a square plate. The left, right, and bottom edges of the plate are maintained at a constant 100°C, whereas the top of the plate is maintained at a constant 0°C. The temperature distribution in the interior of the plate is governed by the Laplace equation shown earlier, where u represents temperature. The Laplace equation is repeated here for convenience:

A finite difference mesh for this problem might look something like that shown in Figure 12-2.

Figure 12-2. Finite difference grid for temperature distribution in square plate

Each little circle shown in the grid in Figure 12-2 represents a computational node where the following finite difference equation, from the previous recipe, applies:

In this example, I'm assuming the grid spacing, h, is uniform in both x and y directions and equal to 1 cm.

If you rearrange this finite difference equation, solving for u(x, y), you get the following:

You can see that u (the temperature) at each node is simply the average of the temperatures of adjacent nodes. For nodes adjacent to the plate boundary, the specified boundary conditions are included in the average.

This example results in 49 finite difference equations with 49 unknown temperatures. A temperature distribution is sought that satisfies the above equation at each node.

Another way to look at this problem is to reconsider this form of the finite difference equation:

This is the same equation shown earlier, prior to solving for u(x, y). This form of the equation represents the residual at each nodethat is, the temperature minus the estimated temperature at each node. This residual must ideally be 0. If you take all these residuals, one for each node, and square them and then sum the squares, you end up with a least-squares minimization problem. You can use Solver to minimize the sum of squared residuals to arrive at a solution that satisfies the governing equations. Figure 12-3 shows the spreadsheet I set up to solve this problem in Excel.

Figure 12-3. Finite difference spreadsheet setup

The table under the heading Finite Difference Solution represents the temperature distribution satisfying the governing equations and boundary conditions. At the moment, the table is filled with initial guesses except at the boundaries, where the given boundary conditions are applied. Each of the shaded cells in this table corresponds to a computational node in the finite difference grid illustrated in Figure 12-2. The nonshaded cells bounding the shaded cells correspond to boundary nodes with specified boundary conditions.

The lower table in Figure 12-3, under the heading Residuals-squared, contains a similar set of shaded cells. In this case, the shaded cells contain the finite difference equation for each node, squared. The cell formulas are of the following form: =(-4*S6+S5+S7+R6+T6)^2. This is the cell formula for cell S18. All the other cells contain similar formulas with relative cell references. This form is essentially the square of the residual for each node. Note that these formulas refer to the values for each node contained in the upper table.

Cell R28 contains the formula =SUM(S18:Y24). This is simply the sum of all squared residuals contained in the lower table.

Now you can use Solver to minimize the value in cell R28 by changing all the values contained in the shaded cells in the upper table. Figure 12-4 shows the Solver model I used for this example.

Figure 12-4. Solver model for finite difference solution

You can see that this model aims to minimize the value in cell R28, the sum of squared residuals, by changing all the values contained in cells S6 to Y12. If Solver is successful, cells S6 to Y12 in the upper table in Figure 12-3 will contain a temperature distribution that satisfies the governing equations and boundary conditions.

After you press the Solve button, Solver does indeed find a solution, as shown in Figure 12-5.

Figure 12-5. Solution to finite difference problem

As you can see, the sum of squared residuals is a very small number and the squared residuals for each cell are all 0 (to 14 decimal places). The upper table now contains the solution to the temperature distribution problem.

Just for fun, I plotted these results using Excel's 3-D Surface chart feature (see Chapter 4). Figure 12-6 shows the resulting plot.

This highlights another advantage of using Excel: you can quickly visualize your results without having to program your own visualizations or use a separate visualization tool.

Figure 12-6. Surface plot of temperature distribution