% 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')