Graphics Programs Reference
In-Depth Information
distribution at time step
n
. (At the endpoints
j
=
0 and
j
=
J
, this equation
refers to temperatures outside the prescribed range for
x
, but at these points
we will ignore the equation above and apply the boundary conditions
instead.) We can interpret this equation as saying that the temperature at a
given location at the next time step is a weighted average of its temperature
and the temperatures of its neighbors at the current time step. In other
words, in time
t
, a given section of the wire of length
x
transfers to each
of its neighbors a portion
s
of its heat energy and keeps the remaining
portion 1
−
2
s
of its heat energy. Thus our numerical implementation of the
heat equation is a discretized version of the microscopic description of
diffusion we gave initially, that heat energy spreads due to random
interactions between nearby particles.
The following M-file, which we have named
heat.m
, iterates the
procedure described above:
function u = heat(k, x, t, init, bdry)
% solve the 1D heat equation on the rectangle described by
% vectors x and t with u(x, t(1)) = init and Dirichlet
% boundary conditions
% u(x(1), t) = bdry(1), u(x(end), t) = bdry(2).
J = length(x);
N = length(t);
dx = mean(diff(x));
dt = mean(diff(t));
s = k*dt/dxˆ2;
u = zeros(N,J);
u(1, :) = init;
for n = 1:N-1
u(n+1, 2:J-1) = s*(u(n, 3:J) + u(n, 1:J-2)) +...
(1 - 2*s)*u(n, 2:J-1);
u(n+1, 1) = bdry(1);
u(n+1, J) = bdry(2);
end
The function
heat
takes as inputs the value of
k
, vectors of
t
and
x
values, a
vector
init
of initial values (which is assumed to have the same length as
x
), and a vector
bdry
containing a pair of boundary values. Its output is a
matrix of
u
values. Notice that since indices of arrays in MATLAB must start
at 1, not 0, we have deviated slightly from our earlier notation by letting
n=1
Search WWH ::
Custom Search