Solution of Problem IV with MATLAB using simultaneous solution for the conduction terms and successive substitution for the heat generation terms (Version 1).

The problem to be solved is

The finite difference formulation of this problem is

Henceforth the primes on the T will be dropped for convenience. Rearrange the differential equation to the following form.

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. For a' = 2, b' = 0.5 the equation is

These equations are applied here using MATLAB. The command file is shown below.

% Solution to Problem IV using successive substitution for heat generation.
dx = 1;
% Solve the problem for 5 different x.
for kk=1:5
dx = dx/2;
n = 1 + 1/dx
% Need to iterate.
tol = 1e-12;
% set up initial guess
for i=1:n
y1(i) = dx*(i-1);
end
%begin iteration
change = 1.;
iter=0;
while change>tol
iter=iter+1;
sum = 0.;
% Calculate the matrices
for i=1:n
heat(i) = 2*dx*dx*exp(0.5*y1(i));
end
for i=2:n-1
f(i) = -heat(i);
for j=1:n
aa(i,j) = 0.;
end
aa(i,i-1) = 1;
aa(i,i) = -2;
aa(i,i+1) = 1;
end
% First and last row of matrices are special.
f(1) = 0;
f(n) = 1;
aa(1,1) = 1;
aa(1,2) = 0;
aa(n,n-1) = 0;
aa(n,n) = 1;
% Solver for this iteration.
y2 = aa\f';
% Calculate the quantities to test for convergence.
sum=0.;
for i=1:n
sum = sum + abs(y1(i) - y2(i));
y1(i) = y2(i);
end
change = sum
% Test to see if change < tol.
end
% Save some results to plot
nm = (n-1)/2 + 1
x=dx*(nm-1)
ys(kk) = y1(nm);
xs(kk) = dx*dx;
iters(kk)=iter
end

% Plot the value of y(0.5) versus x^2.
plot(xs,ys,'-*')
xlabel('x^2')
ylabel('y(0.5)')
clear aa
clear f
clear y1
clear y2
ys
xs
iters

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

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 few commands print the solution and number of iterations.

How do we check our work? First, we set heat(i) = 0, which is the same as taking the no heat generation. In that case the solution should be a straight line, and it is. If we need to troubleshoot, 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. When we have a solution, we can take the values and substitute them into the equation they are supposed to satisfy. If they do, then the equations are satisfied. You have to use judgment to decide how many of the equations you need to test; it depends on how the MATLAB program is written. Next, we solve using successively smaller ∆x and see if the solution changes. The results are (for tolerance = 10^-12):

∆x	y(0.5)	No. of iterations
0.5	0.890152	18
0.25	0.880706	17
0.125	0.878543	17
0.0625	0.878013	17
0.03125	0.877882	17

The numbers are close to each other. Furthermore the increments between the answers are:

(0.5, 0.25): 0.890152 ­ 0.880706 = 0.0945

(0.25, 0.125): 0.880706 ­ 0.878543 = 0.00216

(0.125, 0.0625): 0.878543 ­ 0.878013 = 0.000533

(0.0625, 0.03125): 0.878013 - 0.877882 = 0.000131

The error should go as ∆x'²; 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.000533 / 4 = 0.000134, which is close enough to say the convergence is followed. Extrapolation to zero gives T'(0.5) = 0.877836. The figure also shows convergence which is proportional to ∆x².

Notice that this solution is also the same derived using Excel. 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.

This case is a particularly hard one to solve when b' is large. In fact, for the case ∆x = 0.5 it can be shown that a solution (real) does not exist for b' larger than 0.926112 (proof). You can try this value yourself and see what happens. This is an example that underscores the necessity to carefully check your program. If you had not done that, the logical assumption would be that you had made a mistake. However, if the program is carefully checked then the results can be depended upon. In this case, no solution is possible, but to know that requires other analysis (link).

How do we know the accuracy of our result? Firstly, we have run the problem with five different ∆x and find similar results. Secondly, we calculated the truncation error, found the error was proportional to ∆x², and found that the results were also proportional to ∆x² for the smallest ∆x. These two results confirm that we have solved some problem correctly. Thirdly, we carefully checked every line of the code and made sure that we solved the equation we intened to solve. Fourthly, we extrapolated the results in ∆x² to zero to obtain the answer, 0.877836. Thus, the last calculation, with ∆x = 0.03125, was accurate to within 0.000046, or with 0.005% accuracy.

Take Home Message: When solving nonlinear finite difference equations, it is sometimes possible to evaluate the nonlinear terms using the past iterate, thereby making the programming easier. Usually, more iterations are then necessary.