Fall 2013

Math 471

Problem Set 5 (10/16, 18, 21)
Due on Fri, Nov 11) Let us consider the two-dimensional boundary value problem of a metal plate on
D = {(x, y); 0 < x, y < 1}. The plate temperature φ(x, y) obeys −∆φ = 0 on Dwith boundary conditions φ(x, 1) = 1, φ(x, 0) = φ(0, y) = φ(1, y) = 0. This means
there are no heat sources inside the plate, and one side of the plate is kept at ahigh temperature while the other three sides are kept at a low temperature.
We will obtain φ(x, y) inside the plate. Let w be the numerical solution, i.e.,yyxxwij ≈ φ(xi , yj ). By the finite-difference scheme, we have (D+D−+ D+D−)wij = 01
1 1 1
with mesh size h = n+1 , for h = 4 , 8 , 16 . This yields a linear system Aw = f .
The mesh points are given by (xi , yj ) = (ih, jh), i, j = 0, 1, . . . , n, n + 1. The
finite-difference equations can be written in component form as1(4wij − wi+1,j − wi−1,j − wi,j+1 − wi,j−1 ) = fij .
h2Thus for example by Jacobi’s method, we have
1(k)
(k)
(k)
(k)
(k+1)
− wi+1,j − wi−1,j − wi,j+1 − wi,j−1 ) = fij ,(4wij2h(k)where wij is the numerical solution at step k.
(a) Solve the problem by Jacobi’s method and the Gauss-Seidel method. Do notform the full matrix A (because it is sparse and that would be inefficient). If
you like, you can use the Matlab template below. Submit a copy of the code.
(b) For each value of h, plot the computed temperature wij at the final step(including the boundary values) using a contour plot and a mesh plot (typehelp contour and help mesh for instructions).
(c) Present the following results in a table. column 1: h, column 2: number of
iterations needed to reach the stopping criterion.
(d) What is the value of the temperature at the corners of the plate in the limith → 0? Explain your answer.
A few tips: (1) Since Matlab doesn’t accept zero indices, take i=1:n+2, j=1:n+2.
(2) In the case of the two-point boundary value problem in one dimension, weput the boundary values in f . However in a two-dimensional problem, it is moreconvenient to keep the boundary values in w. Therefore the boundary values andinterior values of w are set at the initial step and the interior values are updated atevery new step. The boundary values can be stored in elements of w with indices
i=1,n+2, j=1,n+2. (3) The interior values are set to zero at the initial step. (4)You can use the stopping criterion krk k2 /kr0 k2 ≤ 10−4 , where rk = f − Awk is theresidual at step k.

function bvp2dclear; clf;
tol = ...; % set tolerance for stopping criterion
for icase=1:3
n = 2∧ (icase+1)-1; h = 1/(n+1); % set mesh sizex = 0:h:1; y = 0:h:1; % create x and y arrays for plots% initialize solution and residual arrays
w new = zeros(n+2,n+2);
w old = zeros(n+2,n+2);
res = zeros(n+2,n+2);% set nonzero boundary values
for j = ...; w new(...,...) = ...; w old(...,...) = ...; end% initialize control variables
k = 0; ratio = 1;
% start iteration
while ratio > tolk = k+1;% compute residual vector
for i = ...; for j = ...;
res(i,j) = ...;
end; end% compute ratio of residual norms
rn(k) = norm(res,’fro’);ratio = rn(k)/rn(1);% compute numerical solution
for i = ...; for j = ...;
w new(i,j) = ...;
end; end;
w old = w new; % reset numerical solution for next stepend % end while% store results for output
table(icase,1) = h; table(icase,2) = k;
% draw contour plotsubplot(2,3,icase)contour(x,y,w new); axis square
string = sprintf(’h=1/%d’,n+1); title(string)% draw surface plot
subplot(2,3,3+icase)mesh(x,y,w new)string = sprintf(’h=1/%d’,n+1); title(string)endtable