% This code solve the one-dimensional heat diffusion equation
% for the problem of a bar which is initially at T=Tinit and
% suddenly the temperatures at the left and right change to
% Tleft and Tright.
%
% Upon discretization in space by a finite difference method,
% the result is a system of ODE's of the form,
%
% u_t = Au + b
%
% The code calculates A and b. Then, these are used in
% various different time integration methods to evolve the
% temperature in time. See lecture notes for more info.
%
% Using this code, we will talk about the following issues:
%
% * Use of explicit methods for stiff systems
% * Implementation of implicit methods for linear systems
% * Sparse vs. full matrices
% * Behavior of trapezoidal when lambda dt -> -infinity.
clear all;
close all;
% Set non-dimensional thermal coefficient
k = 1.0; % this is really k/(rho*cp)
% Set length of bar
L = 1.0; % non-dimensional
% Set initial temperature
Tinit = 400;
% Set left and right temperatures for t>0
Tleft = 800;
Tright = 1000;
% Set up grid size
Nx = 1000;
h = L/Nx;
x = linspace(0,L,Nx+1);
% Set timestep size
dt = 1e-1;
% Calculate number of iterations (Nmax) needed to iterate to t=Tmax
Tmax = 0.5;
Nmax = floor(Tmax/dt);
% Initialize a sparse matrix to hold stiffness & identity matrix
A = spalloc(Nx-1,Nx-1,3*(Nx-1));
I = speye(Nx-1);
% These commands allocate full (non-sparse) matrices
%A = zeros(Nx-1,Nx-1);
%I = eye(Nx-1);
% Calculate stiffness matrix
for ii = 1:Nx-1,
if (ii > 1),
A(ii,ii-1) = k/h^2;
end
if (ii < Nx-1),
A(ii,ii+1) = k/h^2;
end
A(ii,ii) = -2*k/h^2;
end
% Set forcing vector
b = zeros(Nx-1,1);
b(1) = k*Tleft/h^2;
b(Nx-1) = k*Tright/h^2;
% Set initial vector
v = Tinit*ones(Nx-1,1);
% Iterate in time
for n = 1:Nmax,
% Backward Euler
%v = (I - dt*A)\(v + dt*b);
% Forward Euler
%v = v + dt*(A*v+b);
% Trapezoidal
rhs = v + 0.5*dt*A*v + dt*b;
G = I - 0.5*dt*A;
v = G\rhs;
T = [Tleft; v; Tright];
plot(x,T,'*');
axis([0,1,400,1200]);
drawnow;
end