Solution of Problem III in heat flux-temperature form with MATLAB using the Jacobi iteration (Version 2).

The problem to be solved is Problem III.

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

The finite difference formulation of this problem is

y1 = 0, yN = 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 MATLAB. The command file is shown below.

% Solve Problem III in Matlab using version2
tol = 1e-8;
delx = 0.5;
n = 1/delx + 1;
y1(1) = 0;
y1(n) = 1;
% set up initial guess
for i=2:n-1
y1(i) = delx*(i-1);
end
%begin iteration
change = 1.;
iter=0;
while change>tol
iter=iter+1;
sum = 0.;
y2(1) = 0.;
y2(n) = 1.;
for i=1:n
al(i) = 1 + y1(i);
end
for i=2:n-1
y2(i) = (al(i+1)+al(i))*y1(i+1) + (al(i)+al(i-1))*y1(i-1);
y2(i)=y2(i)/(al(i+1)+2*al(i)+al(i-1));
sum=sum+abs((y2(i)-y1(i)));
end
% transfer y2 -> y1 for next iteration
for i=1:n
y1(i) = y2(i);
end
change = sum
% test change < tol
end
y2
iter

The program initially sets the tolerance, grid size, and number of points. Then the first guess for y(x) is taken as a straight line. The iteration will continue until change is less than the desired tolerance. Change is defined as

The expression for y2(i) is simply the iterative equation, while sum is the cumulative sum of absolute values of the changes in y at each node. Finally, the new value, y2, is transferred to the old value, y1, and sum is transferred to change so that the 'while' change can be tested. Finally, the last two commands print the solution and number of iterations.

How do we check our work? First, we remove the ';' from the statements, use delx = 0.25, and do the calculations by hand to compare with the MATLAB results. This insures that the calculations follow the equations given above. Next, we solve using successively smaller ∆x and see if the solution changes. The results are (for tolerance = 10-8):

∆x y(0.5) No. of iterations
0.5 0.58113883 7
0.25 0.58113882 46
0.125 0.58113881 190
0.0625 0.58113878 739

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 and the other methods of MATLAB, version 1 and MATLAB, version 3. 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: Solving for a single variable in a step provides a method which is slow to converge. Version 3 is faster.