Solution to Problem III using Excel in the Temperature Formulation and Simultaneous Solution of the Equations.

The problem to be solved is Problem III.

y(0) = 0, y(1) = 1

The finite difference formulation of this problem is

y₁ = 0, y_N = 1

When the transport coefficients at the mid-points are evaluated using the average of the points on either side, the equation becomes

The iterative scheme is then

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. A partial check of the result can be made by assuming that the transport coefficient a is constant; in that case the result should be the one derived for Problem II with b = 0, which it is:

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_i at each grid point. The second row is the transport coefficient, 1 + y_i. The third row is the position, x, and the cell A4 is the value of ∆x. The value of ∆x is placed in cell A4; first use 0.5. The first value of x, 0, is placed in the cell A3. Then the equation is written in cell B3 and copied to the right.

B3: = A3+$A$4, copied to the right

The $A$4 insures that when the cell formula is copied to the right, the next cell has B3+$A$4, then C2+$A$4; in other words, the A4 location is absolute and is not incremented during the copying. Once this is done, the third row should show the values 0, 0.5, 1.0. The second row is used for the value of the transport coefficient at the i-th node.

A2: = 1 + A1, copied to the right

The formulas for the differential equation are placed into the cells of the first row.

B1: = ((B2+C2)*C1 + (A2+B2)*A1) / ( C2 + 2*B2 + A2 ), 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.5811388	1.000000
2	1	1.5811388	2
3	0	0.5	1

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 third row, which shows the location of the grid points. They should be x = 0, 0.5, 1.0, and these are correct. Next we check the second row, which should be 1 + y, using y from the first row; this is also true. 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.3228757	0.5811388	0.802776	1
8	1	1.3228757	1.5811388	1.802776	2
9	0	0.25	0.5	0.75	1
10	0.25

Note that the value at ∆x = 0.5, y(0.5) = 0.5811388, is the same as 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. The answers obtained for y(0.5) are shown below (using a tolerance of 10^-8).

∆x	y(0.5)
0.5	0.58113883
0.25	0.58113883
0.125	0.581138876

The numbers are so close to each other than no further analysis is necessary.

Notice that this solution is also the same derived using Excel with the first iteration techniqueas well as 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. This version of the equations (in heat flux form) gives a more accurate answer for the same mesh size than version 1.