Example 5.1

Inhomogeneous heat equation

Contents

Setup

close all;
N = 21;
L = 1;
alpha = 1;
dx = L/(N-1);
X = linspace(0,1,N)';

% initial condition
T = sin(pi*X);

Spatial derivative operator

A = gallery('tridiag',N-2,1,-2,1);

inhomogeneous term

f = @(t) (pi^2-1)*exp(-t)*sin(pi*X(2:N-1));

time advancement

g = @(t,dt,T) T(2:N-1)+alpha*dt/(dx^2)*A*T(2:N-1) + dt*f(t);

Stable run

dt = 0.001;           % time step
t_final = 2.0;        % final time
time = 0:dt:t_final;  % time array
pt = [0.0 0.5 1.0 1.5 2.0];   % desired plot times
pn = length(pt);          % number of desired plots
pc = 1;               % plot counter
rt = zeros(1,pn);
S = zeros(N,pn);      % solution storage


for t = time
    % plot storage
    if ( t >= pt(pc) )
        S(:,pc) = T;
        rt(pc) = t;
        pc = pc + 1;
        if (pc > pn)
            break
        end

    end

    % time advancement
    T(2:N-1) = g(t,dt,T);
end

Plot stable run

figure(1)
plot(X,S,'LineWidth',1)
xlabel('x','FontSize',14)
ylabel('T(x)','FontSize',14)
legend('t = 0.0','t = 0.5','t = 1.0','t = 1.5','t = 2.0')
axis([0 1 0 1.1])
ax = gca;
set(ax,'XTick',[0 .25 .5 .75 1])
set(ax,'YTick',[0 .25 .5 .75 1])
grid on

Unstable run

T = sin(pi*X);            % reset initial condition
dt = 0.0015;              % time step
t_final = .3;             % final time
time = 0:dt:t_final;      % time array
pt = [0.0 0.15 0.153 0.158 0.166];   % desired plot times
pn = length(pt);  % number of desired plots
pc = 1;                   % plot counter
rt = zeros(1,pn);
S = zeros(N,pn);          % solution storage


for t = time
    % plot storage
    if ( t >= pt(pc) )
        S(:,pc) = T;
        rt(pc) = t;
        pc = pc + 1;
        if (pc > pn)
            break
        end
    end
    % time advancement
    T(2:N-1) = g(t,dt,T);
end

Plot unstable run

figure(2)
plot(X,S,'LineWidth',1)
xlabel('x','FontSize',14)
ylabel('T(x)','FontSize',14)
legend('t = 0.0000','t = 0.1500','t = 0.1530','t = 0.1590','t = 0.1665')
axis([0 1 0 1.3])
ax = gca;
set(ax,'XTick',[0 .25 .5 .75 1])
set(ax,'YTick',[0 .25 .5 .75 1 1.25])
grid on

Published with MATLAB® 7.5