% Solve Heat Eqn with Fourier Series and plot temperature
% distribution as a function of time. -A. Ruina
% Both bar ends are at T=0
% u is the temperature.
% Choice of 3 starting temp distributions.
% The only tricky Matlab is the 'erasemode', 'xor', drawnow stuff
% described in Pratap book.
%Set up basic constants
fprintf('Calculating ... ')
t_tot =.1; % total time of integration
ntimes=500; % temporal resolution
numbpts = 100; % spatial resolution
nterms = 100; % number of terms in Fourier series
L=1; %length of bar
beta=1; %thermal diffusivity
x=(0:numbpts)*L/numbpts; %value of x at discrete points
%Contribution to u comes from the Fourier terms
%The matrix A has one row for each term, each column is a
%value of x
A=zeros(numbpts+1,nterms); % initialization
for n=1:nterms
A(:,n)=[sin(n*pi*x/L)]';
end
%Set initial temperature distribution
%fx is the initial temp distribution
%Pick from one of the three below
%fx=x*0 +1; % Constant initial temperature
fx= exp(-100*(x-.5).^2); % Delta function, approximately
%fx=sin(1*pi*x/L); % sine wave
%Find Fourier coefficients
%Crude numerical integration to find Fourier coefficients
for n = 1:nterms;
a_n = 2 * fx*A/numbpts;
end
u0 = A*a_n'; %Add the rows of A to get initial temp dist.
h1 = plot(x,u0 ,'erasemode','xor'); %fancy way to plot
xlabel('x, position on bar');
ylabel('u, temperature')
title('Temperature vs position, evolving in time.')
%Plot temp at a sequence of times
for t=linspace(0,t_tot,ntimes);
%Add the rows of A, attenuated by exponential decay, to get u
u= A*[a_n.* exp(-beta*(1:nterms).^2 *pi^2*t/L^2)]';
set(h1,'ydata', u);drawnow; % fancy plot, continued
end
fprintf(' done.\n\n')