exact_solution.m
代码语言：javascript复制function ye = exact_solution(t,x,c)
% Function called: profile
yyy = x;
for j = 1:length(x), yyy(j) = profile(x(j) - c*t); end
ye = yyy;
kappa_scheme.m
代码语言：javascript复制kappa = 0; % Parameter of kappa-scheme
c = 1; % Velocity
sigma = 0.7; % CFL number
J = 31; % Number of grid cells
bc = 2; % Governs choice of numerical boundary conditions; see page 350
% Enter 1 for boundary conditions using upwind schemes
% or 2 for extrapolation using equation (9.36)
T = 0.4; % Final time
% ..............End of input...................................
h = 1.0/(J-1); % Size of grid cells
dt = sigma*h/c; % Time step
n = floor(T/dt); % Number of time steps
% Definition of grid numbering
% x=0 x=1
% grid |------|------|------|------|
% 1 2 3 J
x = [0:J-1]'*h; % Location of grid nodes
sig2 = sigma^2; sig4 = sigma/4;
ynew = zeros(J,1); yright = zeros(n 1,1);
t = 0; yright(1) = exact_solution(t,1,c); yold = exact_solution(t,x,c);
for i = 1:n, t = t dt;
a = 2*yold(1) - yold(2); % Extrapolated inflow value
ynew(1) = exact_solution(t,0,c);
if bc == 1
ynew(2) = yold(2) - sigma*(yold(2) - yold(1));
else
ynew(2) = yold(2) sig4*(-1 kappa (1-3*kappa)*sigma 2*kappa*sig2)*a ...
sig4*(5-3*kappa (9*kappa-1)*sigma-6*kappa*sig2)*yold(1)...
sig4*(-3 3*kappa -(1 9*kappa)*sigma 6*kappa*sig2)*yold(2) ...
sig4*(-1-kappa (1 3*kappa)*sigma-2*kappa*sig2)*yold(3);
end
for j = 3:J-1
ynew(j) = yold(j) sig4*(-1 kappa (1-3*kappa)*sigma 2*kappa*sig2)*...
yold(j-2) sig4*(5-3*kappa (9*kappa-1)*sigma-6*kappa*sig2)*yold(j-1)...
sig4*(-3 3*kappa -(1 9*kappa)*sigma 6*kappa*sig2)*yold(j) ...
sig4*(-1-kappa (1 3*kappa)*sigma-2*kappa*sig2)*yold(j 1);
end
if bc == 1 % Fully upwind scheme at outflow
ynew(J) = (1-1.5*sigma 0.5*sig2)*yold(J) sigma*(2-sigma)*yold(J-1)...
sigma*0.5*(sigma-1)*yold(J-2);
else % Extrapolation at outflow
ynew(J) = yold(J) sig4*(-1 kappa (1-3*kappa)*sigma 2*kappa*sig2)*...
yold(J-2) sig4*(5-3*kappa (9*kappa-1)*sigma-6*kappa*sig2)*yold(J-1)...
sig4*(-3 3*kappa -(1 9*kappa)*sigma 6*kappa*sig2)*yold(J) ...
sig4*(-1-kappa (1 3*kappa)*sigma-2*kappa*sig2)*(2*yold(J) - yold(J-1));
end
yright(i 1) = ynew(J); yold = ynew;
end
yexact = exact_solution(t,x,c); error = ynew - yexact;
norm1 = norm(error,1)/J % Compute error norms
norm2 = norm(error,2)/sqrt(J)
norminf = norm(error,inf)
figure(1), clf, hold on
xx = 0:0.005:1; plot (xx,exact_solution(t,xx,c),'-'); plot (x,ynew,'o');
s1 = ['kappa=',num2str(kappa),' bc =', num2str(bc),' sigma=',...
num2str(sigma),' cells=',num2str(J-1),' t=',num2str(t)];
title(s1,'fontsize',18);
Lax_Wendroff_scheme.m
代码语言：javascript复制c = 1; % Velocity
sigma = 0.7; % CFL number
J = 31; % Number of grid cells
bc = 1; % Governs choice of numerical boundary conditions; see page 350
% Enter 1 for boundary conditions using upwind schemes
% or 2 for extrapolation using equation (9.36)
T = 0.4; % Final time
% ..............End of input...................................
h = 1.0/(J-1); % Size of grid cells
dt = sigma*h/c; % Time step
n = floor(T/dt); % Number of time steps
% Definition of grid numbering
% x=0 x=1
% grid |------|------|---.........---|------|
% 1 2 3 J
x = [0:J-1]'*h; % Location of grid nodes
sig2 = sigma^2;
t = 0;
yold = exact_solution(t,x,c); ynew = zeros(J,1);
yright = zeros(n 1,1); yright(1) = exact_solution(t,1,c);
for i = 1:n, t = t dt; ynew(1) = exact_solution(t,0,c);
for j = 2:J-1
ynew(j) = 0.5*(sig2 sigma)*yold(j-1) (1-sig2)*yold(j) ...
0.5*(sig2-sigma)*yold(j 1);
end
if bc == 1, % First order upwind at outflow
ynew(J) = yold(J) sigma*(yold(J-1) - yold(J));
else % Extrapolation at outflow
ynew(J) = 0.5*(sig2 sigma)*yold(J-1) (1-sig2)*yold(J) ...
0.5*(sig2-sigma)*(2*yold(J) - yold(J-1));
end
yright(i 1) = ynew(J); yold = ynew;
end
yexact = exact_solution(t,x,c); error = ynew - yexact;
norm1 = norm(error,1)/J % Compute error norms
norm2 = norm(error,2)/sqrt(J)
norminf = norm(error,inf)
figure(1), clf, hold on
xx = 0:0.005:1; plot (xx,exact_solution(t,xx,c),'-'); plot (x,ynew,'o');
s1 = ['Lax-Wendroff',' bc =', num2str(bc),' sigma=',...
num2str(sigma),' cells=',num2str(J-1),' t=',num2str(t)];
title(s1,'fontsize',18);
