Solution to Problem III using Excel in the Temperature Formulation and Jacobi iteration.
The problem to be solved is Problem III.
The finite difference formulation of this problem is
Rearrange the differential equation to the following form.
or
The iterative scheme can be represented as
where the superscript k represents the value of the quantity at the k-th iteration. Thus, one takes values at the k-th iteration and calculates the values with the superscript k+1. These values are then used in the equation again, and again. The process repeats until the values no longer change, to within a specified tolerance specified by you.
These equations are applied here using Microsoft Excel Spreadsheet, but the commands are generic ones that most spreadsheets allow. Consider the grid shown in the Figure.
The first row is reserved for the value of y at each grid point. The left-most cell is the value of y at the first point, and the right-most cell is the value of y at the last point. Since these values are known we simply insert the known values in the cells. For the points in between, we insert the equation given above.
It is convenient if we define some quantities, first, though, to simplify the algebra. Let the second row represent the value of x at the node i, which is used for convenience. The value of ∆x is placed in cell A3; first use 0.5. The first value of x, 0, is placed in the cell A2. Then the equation is written in cell B2 and copied to the right.
B2: = A2+$A$3, copied to the right
The $A$3 insures that when the cell formula is copied to the right, the next cell has B2+$A$3, then C2+$A$3; in other words, the A3 location is absolute and is not incremented during the copying. Once this is done, the second row should show the values 0, 0.5, 1.0. The second row is not used for this problem.
The formulas for the differential equation are placed into the cells of the first row.
B1: = 0.5*(A1 + C1) + (C1 A1)^2 / ((1 + B1) * 8), copied to the right.
When you are placing these formulas in your spreadsheet, the program, will give you a warning since you are defining a circular reference, that is some cells are defined in terms of other cells, which are in tern defined in terms of the the first cells. To begin the calculation in Microsoft Excel you have to choose Tools, then Preferences, then Calculation, and check the Iteration Box. That tells the computer that you will be performing iterations. (You may wish to prepare the spreadsheet with 'Iteration' turned off and turn it on when preparation is complete. This will prevent getting weird answers during the preparation when all the cells are not complete. This is especially true for non-linear problems.) There you can also set the maximum number of iterations to 500, and the convergence parameter to 10-6. After you choose Calculate Now the spreadsheet should look like this.
| A | B | C | |
| 1 | 0 | 0.5791562 | 1.000000 |
| 2 | 0 | 0.5 | 1 |
| 3 | 0.5 |
How do we check our work? First, we make sure that the numbers we put in are correct. We check cell A4, which has the value of ∆x; 0.5 is correct. Then we check the second row, which shows the location of the grid points. They should be x = 0, 0.5, 1.0, and these are correct. Lastly we check the cell B1. We can use the numbers in A1 and C1, since the cell formula does, too. We get
Since these calculations agree, we conclude we have set the spreadsheet up correctly.
Next do the same thing with a ∆x = 0.25. The spreadsheet results are shown in the Table.
| A | B | C | D | E | |
| 7 | 0 | 0.3221458 | 0.580559972 | 0.8024805 | 1 |
| 8 | 0 | 0.25 | 0.5 | 0.75 | 1 |
| 9 | 0.25 |
Note that the value at ∆x = 0.5, y(0.5) = 0.5806, is very close to the one derived using ∆x = 0.5. To be absolutely sure that the spreadsheet is prepared correctly, however, this spreadsheet is checked in the same way. Finally, use ∆x = 0.125 and then 0.0625. The answers obtained for y(0.5) are shown below (using a tolerance of 10-8).
| ∆x | y(0.5) |
| 0.5 | 0.579156 |
| 0.25 | 0.58056 |
| 0.125 | 0.580986 |
| 0.0625 | 0.5811 |
The numbers are close to each other. Furthermore the increments between the answers are:
(0.5, 0.25): 0.579156 0.580560 = 0.00140
(0.25, 0.125): 0.580560 0.580986 = 0.000426
(0.125, 0.0625): 0.580986 0.581100 = 0.000114
The error should go as ∆x2; with a ∆x half as big, the error should be one-fourth as big. If we are in the region (small enough ∆x) where this applies, the differences should decrease by a factor of 4 each time. In this case 0.00140 / 4 = 0.00035, which isn't quite 0.00043. However, at the next step 0.000426 / 4 = 0.000107, which is close enough to 0.000114 to say the convergence is followed. This shows that the formulas put into the spreadsheet are correct second-order expressions. Richardson extraplolation gives y(0.5) = 0.581138 as the answer.
Notice that this solution is also the same derived using MATLAB. That is as it should be, since we are solving the same equations. The only reason for a difference is that we might not have set the tolerance to a small enough value. Needless to say, the more accurate we want the solution, the smaller is the tolerance, and the more iterations we need.
Take Home Message: A spreadsheet program provides an easy way to solve the finite difference equations. Version 2 of the equations (in heat flux form) gives a more accurate answer for the same mesh size.