Example 6.13

Convection Equation with non constant coefficients

Contents

Initialization

close all;
clear N h u x j X U D k A


N  = 16;
h  = 0.001;

u  = @(x,t) sin(2*pi*x*exp(-2*t));
x  = [-1:0.01:1];

j  = [0:N]';
X  = cos(j*pi/N);
U  = u(X,0);

D  = zeros(N+1,N+1);

Creation of the matrix operator

for j = [1:N+1]
    for k = [1:N+1]
        if j==k
            switch j
                case 1
                    D(j,k) = (2*N^2+1)/6;
                case N+1
                    D(j,k) = -(2*N^2+1)/6;
                otherwise
                    D(j,k) = -X(j)/(2*(1-X(j)^2));
            end;
        else
            if or(j==1,j==N+1) cj = 2; else cj = 1; end;
            if or(k==1,k==N+1) ck = 2; else ck = 1; end;
            D(j,k) = cj*(-1)^(j+k)/(ck*(X(j)-X(k)));
        end;
    end;
end;

Initial condition

plot(x,u(x,0),'k');
axis([-1 1 -1.5 1.5]);
hold on;

Iterate and plot

A  = eye(N+1) - 2*h*diag(X)*D;

for t = [h:h:0.6]
    U = A*U;

    if t == 0.3
        plot(x,u(x,t),'k--',X,U,'ko');
    end
    if t == 0.6
        plot(x,u(x,t),'k--',X,U,'k*');
    end;
end;

Published with MATLAB® 7.5