% Matlab script: dif1d_main.m
%
% 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, uses one of Matlab's
% ODE integrators, either ode45 (which is based on a Runge-Kutta
% method and is not designed for stiff problems) or ode15s (which
% is based on an implicit method and is designed for stiff problems).
%
clear all;
close all;
sflag = input('Use stiff integrator? (1=yes, [default=no]): ');
% 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 = input(['Enter number of divisions in x-direction: [default=' ...
'51]']);
if (isempty(Nx)),
Nx = 51;
end
h = L/Nx;
x = linspace(0,L,Nx+1);
% Calculate number of iterations (Nmax) needed to iterate to t=Tmax
Tmax = 0.5;
% Initialize a sparse matrix to hold stiffness & identity matrix
A = spalloc(Nx-1,Nx-1,3*(Nx-1));
I = speye(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
v0 = Tinit*ones(1,Nx-1);
if (sflag == 1),
% Call ODE15s
options = odeset('Jacobian',A);
[t,v] = ode15s(@dif1d_fun,[0 Tmax],v0,options,A,b);
else
% Call ODE45
[t,v] = ode45(@dif1d_fun,[0 Tmax],v0,[],A,b);
end
% Get midpoint value of T and plot vs. time
Tmid = v(:,floor(Nx/2));
plot(t,Tmid);
xlabel('t');
ylabel('T at midpoint');