Contents

clear all;
close all;

c=1;
delx=.01;
a=-3;
b=5;
x=a:delx:b;
J=(b-a)/delx; % set number of spatial steps. The nr. of dividing points is J+1
delt=.01;
r=c*delt/delx; % c/rho.
nsteps=100;
tfinal=nsteps*delt;


u=f(x); % initialize the vector u

Forward difference

figure

for n=1:nsteps
    v=u; % v stores the solution at the old time level

    for j=1:J
        u(j)=(1+r)*v(j)-r*v(j+1);
    end;
    %u(J)=v(J)+delt*v(1); %boundary condition

    %u(J)=1;
    if(n==1)
        plot(x,u,'g-+');
    end;
    hold on;
    if(n==nsteps)
        plot(x,u,'r-s')
        hold on
        plot(x,sol(c,x,delt*nsteps)); % plot the solution with the method of characteristics
        %xlim([0 2])
        %ylim([-1 1])
    end;
end;

%say1=['tfinal=' num2str(tfinal)];
title('The solution using a forward difference approx.')
xlabel('x','FontSize',14);
ylabel('u_{fwd}(x,t)','FontSize',14);
legend('u(x,t_1)', 'u(x,tfinal)','u_{exact}(x,tfinal)')
hold off;

Upwind scheme

u=f(x);

figure

for n=1:nsteps
    v=u; % v stores the solution at the old time level

    for j=2:J+1
        u(j)=(1-r)*v(j)+r*v(j-1);
    end;
    %u(J+1)=1;

    if(n==1)
       plot(x,u,'g-+');
    end;
    hold on;
    if(n==nsteps)
        plot(x,u,'r-s')
        hold on
        plot(x,sol(c,x,delt*nsteps)); % plot the solution with the method of characteristics
    end;
end;

title('The solution using an upwind difference approx.')
xlabel('x','FontSize',14);
ylabel('u_{uwd}(x,t)','FontSize',14);
legend('u(x,t_1)', 'u(x,tfinal)','u_{exact}(x,tfinal)')
hold off;

The solution with the method of characteristics

figure
plot(x,f(x),'r'); % plot the initial solution
hold on
plot(x,sol(c,x,1)); % plot the solution with the method of characteristics
hold off;
title('Initial solution and Char.solution u(x(t),t)');
xlabel('x','FontSize',14);
ylabel('u(x,t)','FontSize',14);
legend('f','u_{char}(x,tfinal)')

Published with MATLAB® R2015b