Inhomogeneous Heat Equation example 3

clear
clc
%Terms in Fourier series. 50 is enough. Any accuracy gained from increasing
%this further is probably not worth the extra computation time.
N=50;
%Diffusivity
diff= 2*10^-2;
%Set the number of terms in the x linspace to 100 for generating the gif
%and 10 or 20 for generating the mesh plot.
x=linspace(0,1,100);
y=zeros(1,length(x));
Y=linspace(0,5,5);
%Desired simulation time in seconds. This is ideally set to a multiple of
%10.
simtime=11;
%Desired frames per second. Set this to 1,2, 5, 10, 20, or 50. Set to 1, 2, or 5
%for the mesh plot, higher values for the gif. Framerates higher than 50 do
%not significantly contribute to the quality of the animation and will
%therefore add unnecessary computation time.
FPS=50;
%Number of frames
framenum=simtime*FPS;
%time axis
t=linspace(0,simtime,framenum+1);
%Solution matrix placeholder. Each row will represent T(x,t) for a value of t.
soln=zeros(length(t), length(y));
%Coefficients
forcing_coeffs=zeros(1,N);
phi_coeffs=zeros(1,N);
%Time when heat is turned off.
switchoff=5;
%Magnitude of forcing.
force=5;

for n=1:N
    forcing_coeffs(n) = 2*integral(@(z)force*sin(n*pi*z),0.45,0.55);
end

%Solution while forcin is ON.
for a=1:FPS*switchoff
    y = zeros(1,length(x));
    for n=1:N
        y = y +sin(n*pi*x)*forcing_coeffs(n)*integral(@(s)exp(diff*(n^2)*(pi^2)*(s-t(a))),0,t(a));
    end
    %Row 'a' of solution matrix is equal to the y array. Absolute
    %value to eliminate error causing y to be slightly less than 0 when it
    %should just be 0.
    soln(a,:) = abs(y);
end

%To find the 'initial' condition for the time period when the forcing is
%off, we need a functional representation of the temperature at the moment
%when the forcing is turned off. If you are on MATLAB then you will need
%the Symbolic Toolkit.  SciLab has limited symbolic abilities so I do not
%know if this will work for SciLab.
syms p TXSWITCHOFF(p)
TXSWITCHOFF(p) = 0;
for n=1:N
    TXSWITCHOFF(p) = TXSWITCHOFF(p) + sin(n*pi*p)*forcing_coeffs(n)*integral(@(s)exp(diff*(n^2)*(pi^2)*(s-switchoff)),0,switchoff);
end

%Now we use this to compute the Fourier coefficients for the 'initial' conditions.
phi=matlabFunction(TXSWITCHOFF(p));
for n=1:N
    phi_coeffs(n) = 2*integral(@(z)(phi(z).*sin(n*pi*z)),0,1);
end

%Solution for the temperature after the source has been turned off. We use
%the solution for the temperature at the instant of the switchoff as the
%initial condition, so we have to make the transformation T >> (T-switchoff).
for a=(FPS*switchoff+1):1:length(t)
    y=zeros(1,length(x));
    for n=1:N
        y=y+phi_coeffs(n)*exp(-1*diff*(n^2)*(pi^2)*(t(a)-switchoff))*sin(n*pi*x);
    end
    soln(a,:)=abs(y);
end
%Uncomment the next two lines if you want to store the solution array.
% filename='forcing_solution1.mat';
% save(filename, 'soln')

%Uncomment the following block to generate the frames to be used in the
%animation.

nImages = simtime*FPS;
fig=figure;
FF=zeros(1,length(Y));
for idx = 1:nImages
    clf
    hold on
    plot(x,soln(idx,:))
    plot(FF,Y,'white')
    hold off
    elapsed = (idx/framenum)*simtime;
    if idx < switchoff*FPS
        title(['Elapsed time: ' num2str(elapsed-0.02, '%.2f') ' seconds. Forcing is ON.'])
    else
        title(['Elapsed time: ' num2str(elapsed-0.02, '%.2f') ' seconds. Forcing is OFF.'])
    end
    frame = getframe(fig);
    im{idx} = frame2im(frame);
end
close;

%Uncomment the next two lines if you wish to save a backup of the images you
%that were generated by the last section.
% filename='ForcingSpike1FramesBackup.mat';
% save(filename,im);

%Uncomment the following block to save a gif animation of the solution.

% filename='Forcing_Spike1.gif';
% %Set this to 1 to include every frame, 2 to include every other frame, 4 to
% %include every fourth frame, etc.
% frameskip=1;
% %Animation length in seconds, up to simtime.
% animation_length=10;
% for idx = 1:nImages
%     map = im{idx};
%     A = rgb2gray(im{idx});
%     if idx == 1
%         imwrite(A,filename,'gif','LoopCount',Inf,'DelayTime',1.5);
%     elseif (idx>1) && (idx<=animation_length*FPS) && (mod(idx,frameskip)==0);
%         imwrite(A,filename,'gif','WriteMode','append','DelayTime',(1/FPS)*frameskip);
%     elseif idx==animation_length*FPS+1
%         imwrite(A,filename,'gif','WriteMode','append','DelayTime',1.5);
%     end
% end

%Uncomment the following code to produce a 3D mesh plot of the solution.
% mesh(x,t,soln)
% xlabel('x');
% ylabel('t');
% zlabel('T(x,t) in Celsius');