以单元格为中心的网格对二维定常奇异摄动问题的数值解。

function ye = exact_solution(x,y,D)
z = (1/sqrt(2-x))*(exp(-y*y/(4*D*(2-x))) + exp(-(2-y)*(2-y)/(4*D*(2-x))));
ye = z;

function ye = fL(y, D, hom)
% Neumann condition at outflow
z =  2^(-2.5)*(y*y*0.25/D - 1)* exp(-y*y/(8*D));
z = z + 2^(-2.5)*((2-y)*(2-y)*0.25/D - 1)*exp(-(2-y)*(2-y)/(8*D));
ye = hom*z;

D = 0.001;  % Diffusion coefficient ( = 1/[Peclet number] )
nx = 8;    % Horizontal number of cells
nyc = 4;  % Vertical number of cells outside refinement zone
nyf = 32;  % Vertical number of cells inside refinement zone
central = 0;  % Enter 1 for central scheme or something else for upwind scheme
hom = 0;  % Neumann condition at outflow boundary: 1 for homogeneous,
    % or something else for inhomogeneous (exact) condition
th = 8;    % Governs thickness of refinement zone = del = th*sqrt(D)
    % .....................End of input......................

del = th*sqrt(D);  % Thickness of refinement zone

% Implementation of upwind scheme by artificial diffusion in x-direction
dx = 1/nx;
if central == 1
  Da = D;    % Artificially increased diffusion coefficient
else
  Da = D + dx/2;
end

%     Definition of volume boundaries and grid points
 
%     x = 0    
% grid |---o---|---o---|---o---|---o---|
%      1   1   2   2   3           J

%      |<---- del----->|                      
%     y=0    K1 cells  |        K2 cells      y=1
% grid |-o-|-o-|-o-|-o-|---o---|---o---|---o---|
%      1 1 2 2 3 3         

J = nx; K1 = nyf; K2 = nyc;
x = [1:J]'*dx - dx/2;    % Horizontal coordinates of cell centers

K = K1 + K2;
by = zeros(K+1,1);    % Vertical coordinates of cell boundaries
dy = del/K1;      % Vertical mesh size in refinement zone
for k = 1:K1+1,  by(k) = (k-1)*dy; end
dy = (1-del)/K2;    % Vertical mesh size in unrefined zone
for k = K1+2:K+1,  by(k) = del + (k-K1-1)*dy; end

n=length(by); y=zeros(n-1,1);  % Vertical coordinates of cell centers
dy=zeros(n-1,1);    % Vertical mesh sizes
for k = 1:n-1
  y(k) = 0.5*(by(k)+by(k+1)); dy(k) = by(k+1)-by(k);
end

rhs = zeros(J*K,1);    % Right-hand side of differential equation
for k = 1:K      % Lexicographical ordering of unknowns:
  for j = 1:J      %   j: x-direction  k: y-direction
    rhs(j + (k-1)*J) = q(x(j),y(k),D)*dx*dy(k);
  end
end
  % Modification of rhs due to boundary conditions
  % x = 0: outflow: homogeneous or inhomogeneous Neumann, depending on 
  %    parameter hom
for k = 1:K
  m = 1 + (k-1)*J; rhs(m) = rhs(m) + dy(k)*(Da - 0.5*dx)*fL(y(k),D,hom);
end

  % x = 1: Dirichlet
for k = 1:K
  m = k*J; rhs(m) = rhs(m) + dy(k)*(1 + 2*Da/dx)*exact_solution(1,y(k),D);
end

  % y = 1: homogeneous Neumann, y = 0: Dirichlet
for j = 1:J,  rhs(j) = rhs(j) + 2*dx*D*exact_solution(x(j),0,D)/dy(1); end

  % Matrix construction by flux evaluation
n = J*K; z = zeros(n,1); d = [-J  -1  0  1  J];
A = spdiags([ z  z  z  z  z], d',n,n);     % Declaration of sparsity pattern

for k = 1:K           % Loop over vertical faces
  m1 = (k-1)*J; m = 1 + m1;     A(m,m) = A(m,m) + dy(k);
  a = - dy(k)/2 + dy(k)*Da/dx;  b = - dy(k)/2 - dy(k)*Da/dx; 
  for j = 2:J
    m = j + m1;  
    A(m-1,m-1) = A(m-1,m-1) + a;    A(m-1,m) = A(m-1,m) + b;
    A(m,m-1)   = A(m,m-1)   - a;      A(m,m) = A(m,m)   - b;
  end
  m = J+m1; A(m,m) =  A(m,m) + 2*Da*dy(k)/dx;
end

for j=1:J        % Loop over horizontal faces;
  m = j; A(m,m) = A(m,m) + 2*D*dx/dy(1);
  for k = 2:K
    m = j + (k-1)*J;           a = 2*dx*D/(dy(k) + dy(k-1));
    A(m-J,m-J) = A(m-J,m-J) + a;    A(m-J,m) = A(m-J,m) - a;
    A(m,m-J)   = A(m,m-J)   - a;    A(m,m)   = A(m,m)   + a;
  end
end

numerical_solution = A\rhs;

exsol = zeros (K*J,1);      % Exact solution
for k = 1:K
  for j = 1:J,   exsol(j + (k-1)*J) = exact_solution(x(j),y(k),D);  end
end

% Text output
error = numerical_solution - exsol; n = K*J;
norm1 = norm(error,1)/n
norm2 = norm(error,2)/sqrt(n)
[norminf,m] = max(abs(error))

  % Solution profile at x = x(1)
hh = zeros(K,1); for k = 1:K,   hh(k) = numerical_solution(1 + (k-1)*J);  end
figure (1), clf, subplot(2,2,1), hold on, plot(hh,y,'o')
yy = 0:0.02:1;  exact = zeros (length(yy));
for kk = 1:length(yy),  exact(kk) = exact_solution(x(1), yy(kk), D);  end
plot(exact,yy)

  % Solution profile at y = y(kkk)
vv = zeros(J,1);  kkk = 1;
for j = 1:J,  vv(j) = numerical_solution(j + (kkk-1)*J); end
subplot(2,2,4), plot(x,vv,'o')
for kk = 1:length(yy),  exact(kk) = exact_solution( yy(kk), y(kkk), D); end
hold on, plot(yy,exact)

  % Solution profile at x = x(J)
hh = zeros(K,1); for k = 1:K,  hh(k) = numerical_solution(J + (k-1)*J); end
subplot(2,2,2),  plot(hh,y,'o'),  hold on
for kk = 1:length(yy),  exact(kk) = exact_solution(x(J), yy(kk), D);  end
plot(exact,yy)

  % Contour plot
consol = zeros(K,J);
for k = 1:K,for j = 1:J, consol(k,j) = numerical_solution((k-1)*J + j);end,end
subplot(2,2,3), contour(x,y,consol)

numsol = zeros(K,J);    % waterfall plot
for k = 1:K,for j = 1:J, numsol(k,j) = numerical_solution(j + (k-1)*J);end,end
figure(2),clf, waterfall(x,y,numsol), view(200,70), colormap([0 0 0]), grid off