Example 5.2

Convection Equation

Contents

Setup

close all;
L = 1;
dx = 0.01;
dt = 0.01;
c = 1;

X = (0:dx:L)';
N = length(X);
B = gallery('tridiag',N-1,-1,0,1);
B(end,end) = 2;
B(end,end-1) = -2;

u = exp(-200*(X-0.25).^2);
u(1) = 0;

System of ODEs

f = @(u) [0; -c/2/dx*B*u(2:end)];

Explicit Euler

t_final = .4;
time = 0:dt:t_final;
pt = [0 .12 .25 .375];
pn = length(pt);
pc = 1;
rt = zeros(1,pn);
S = zeros(N,pn);

for t = time
    % plot storage
    if ( t >= pt(pc) )
        S(:,pc) = u;
        rt(pc) = t;
        pc = pc + 1;
        if (pc > pn)
            break
        end
    end
    % time advancement
    u = u + dt*f(u);

end

Plot unstable figure

figure(1)
plot(X,S)
axis([0 .75 -4 4])
ax = gca;
set(ax,'XTick',[0 .25 .5 .75])
set(ax,'YTick',[-4 -2 0 2 4])
grid on
xlabel('x','FontSize',14)
ylabel('u(x)','FontSize',14)
legend('t = 0.00','t = 0.12','t = 0.25','t = 0.38',...
    'Location','NorthWest')

Fourth-order Runge-Kutta

% reset initial condition
u = exp(-200*(X-0.25).^2);
u(1) = 0;

t_final = .8;
time = 0:dt:t_final;
pt = [0 .25 .50 .75];
pn = length(pt);
pc = 1;
rt = zeros(1,pn);
S = zeros(N,pn);

for t = time
    % plot storage
    if ( t >= pt(pc) )
        S(:,pc) = u;
        rt(pc) = t;
        pc = pc + 1;
        if (pc > pn)
            break
        end
    end
    % time advancement
    u1 = dt*f(u);
    u2 = dt*f(u+u1/2);
    u3 = dt*f(u+u2/2);
    u4 = dt*f(u+u3);
    u = u + (u1 + 2*u2 + 2*u3 + u4)/6;

end

Plot stable figure

figure(2)
plot(X,S)
axis([0 1 -.1 1.5])
ax = gca;
set(ax,'XTick',[0 .25 .5 .75 1])
set(ax,'YTick',[0 .25 .50 .75 1 1.25 1.5])
grid on
xlabel('x','FontSize',14)
ylabel('u(x)','FontSize',14)
legend('t = 0.00','t = 0.25','t = 0.50','t = 0.75',...
    'Location','NorthWest')

Published with MATLAB® 7.5