Code covered by the BSD License  

Highlights from
meshfree exponential integrator

meshfree exponential integrator

by

 

experimetal code for meshfree exponential integrators

main
function main
%====================================================================
%            Meshfree exponential Integrator 1D                 
%   Equation: (see function timeEvolution)
%     u_dt = F(u) = epsilon*(u_xx)-alpha*(u_x)+rho*u*(u-1/2)*(1-u) 
%   Reference: 
%        Meshfree exponential integrators
%        M. Caliari, A. Ostermann and S. Rainer
%        SIAM J. Sci. Comput (2013)
%===================================================================



rbf = @(r,e) max(1-e.*r,0).^6.*(35*(e.*r).^2+18*e.*r+3);
rbf_x = @(r,e,dx) -56*dx.*e.^2.*(5*e.*r+1).*max(1-e.*r,0).^5;
rbf_xx = @(r,e) 56*e.^2.*((35*(e.*r).^2-4*e.*r-1).*max(1-e.*r,0).^4);

% shape parameter of RBF
shapeParameter = 15;
% spatial tolerance per time step
tolerance = 1e-3;

borderPoints = [-1,1];
% final integration time
finalTime = 0.1;

% set of candidate intepolation points for first residual subsampling
candidatePoints = linspace(-1,0,10);
candidatePoints = candidatePoints(2:end-1);

% initial value
initfun = @(x) exp(-100*(x+0.5).^2);

[interpolant,points] = calculateRbfInterpolant(candidatePoints,initfun,tolerance);

x = linspace(-1,1,1000);

dt=1e-4;
t = 0;
while t<finalTime
  t = t+dt;
  funTimeEvolution = @(points,tol,eigenValues) timeEvolution(points,dt,...
    interpolant,tol,eigenValues);
  integrationPoints = [points;points+dt];
  
  integrationPoints = sort(unique_tol(integrationPoints,1e-3));
  
  integrationPoints = linspace(-0.9,0.9,200);
  
  % perform time step
  [newInterpolant,points,newDt]=calculateRbfInterpolant(integrationPoints,...
    @(p,eigenValues) funTimeEvolution(p,0.1*tolerance,eigenValues),...
    tolerance,dt);
  if newDt>0
    % find new interpolation points 
    [interpolant,points] = calculateRbfInterpolant(points+6*dt,...
      newInterpolant,tolerance);
    plot(x,interpolant(x));
    hold on;
    % plot interpolation points
    plot(points,zeros(size(points)),'k.');
    axis([-1,1,0,1]);
    s = sprintf('number of points: %d',length(points));
    title(s)
    hold off;
    pause(0.001);
  else
    disp('Timestep rejected');
    t = t-dt;
  end
  dt = abs(newDt);
  if t+dt > finalTime
    dt = finalTime-t;
  end
  dt
end

%=======================================================================
%        Residual Subsampling procedure
%=======================================================================
  function [interpolant,points,dt] = calculateRbfInterpolant(candidatePoints,...
      fun,tolerance,varargin)
    
    candidatePoints = sort(candidatePoints(:));
    numberCandidatePoints = length(candidatePoints);
    pointsCoarsable = [];
    
    % check for error at check points
    checkPoints = calculateCheckPoints(candidatePoints);
    allPoints = [candidatePoints;checkPoints];
    
    if nargout<=2
      rhs = fun(allPoints);
    else
      [rhs,dt,eigenValues] = fun(allPoints,[]);
    end
    
    checkPointsValues1 = rhs(numberCandidatePoints+1:end);
    rhs = rhs(1:numberCandidatePoints);
    
    % calculate inital interpolant
    DM = distanceMatrix(candidatePoints,candidatePoints);
    IM = rbf(DM,shapeParameter);
    coef = IM\rhs;
    interpolant = @(dsites) evaluateInterpolant(dsites,candidatePoints,coef);
    
    % evaluate at check points
    checkPointsValues2 = interpolant(checkPoints);
    
    % find points to be refined
    errorCheckPoints = abs(checkPointsValues2-checkPointsValues1);
    pointsToRefine = checkPoints(errorCheckPoints>tolerance);
    
    % if no points are refined find coarsable points
    if (isempty(pointsToRefine) & nargin<=3)
      testCoarsable = (errorCheckPoints(1:end-1)+errorCheckPoints(2:end))/2;
      pointsCoarsable = candidatePoints(testCoarsable<0.1*tolerance);
    end
    
    while (~isempty(pointsToRefine) | ~isempty(pointsCoarsable))
      if ~isempty(pointsToRefine)
        if ~isempty(pointsCoarsable)
          interpolant = prevInterpolant;
          candidatePoints = prevCandidatePoints;
          break;
        end
        candidatePoints = sort([candidatePoints;pointsToRefine]);
      else
        prevCandidatePoints = candidatePoints;
        prevInterpolant = interpolant;
        
        sortedPointsToCoarse = sortrows([candidatePoints,testCoarsable],2);
        candidatePoints = setdiff(candidatePoints,sortedPointsToCoarse(1));
      end
      
      % check for error at check points
      checkPoints = calculateCheckPoints(candidatePoints);
      allPoints = [candidatePoints;checkPoints];
      numberCandidatePoints = length(candidatePoints);
      
      if nargout<=2
        rhs = fun(allPoints);
      else
        [rhs,dt,eigenValues] = fun(allPoints,eigenValues);
      end
      
      checkPointsValues1 = rhs(numberCandidatePoints+1:end);
      rhs = rhs(1:numberCandidatePoints);
      
      % calculate inital interpolant
      DM = distanceMatrix(candidatePoints,candidatePoints);
      IM = rbf(DM,shapeParameter);
      coef = IM\rhs;
      interpolant = @(dsites) evaluateInterpolant(dsites,candidatePoints,coef);
      
      % evaluate at check points
      checkPointsValues2 = interpolant(checkPoints);
      
      
      % find points to be refined
      errorCheckPoints = abs(checkPointsValues2-checkPointsValues1);
      pointsToRefine = checkPoints(errorCheckPoints>tolerance);
      
      % if no points are refined find coarsable points
      if (isempty(pointsToRefine) & nargin<=3)
        testCoarsable = (errorCheckPoints(1:end-1)+errorCheckPoints(2:end))/2;
        pointsCoarsable = candidatePoints(testCoarsable<0.1*tolerance);
      end
    end
    points = candidatePoints;
  end

%=======================================================================
%   Definition of the PDE
%   time march scheme (exponential Rosenbrock method order 3 (exprb32)
%=======================================================================
  function [u,dt,eigenValues] = timeEvolution(points,dt,interpolant,tol,...
      eigenValues)
    % Equation: u_dt = F(u) = epsilon*(u_xx)-alpha*(u_x)+rho*u*(u-1/2)*(1-u)
    
    epsilon = 1e-5;
    alpha = 6;
    rho = 70;
    
    center = points;
%     points = [center;borderPoints(:)];
    
    DM = distanceMatrix(points,points);
    IM = rbf(DM,shapeParameter);
    
    [R,dummy,p] = chol(sparse(IM),'vector');
    
    if dummy~=0
      error('singular IM matrix: some points probably coincide');
    end
    
    % inflow on left boundary, outflow on the right
    index = [ones(length(points),1)];
%     index(end-1)=0;
    un = interpolant(points);
    u = un;
    
        
    % Linearize at un
    EM = diffusion(DM,points);
    diffusionPart = @(u) evalOperator(u,R,p,EM);
    JnSymm = @(u) (epsilon*diffusionPart(u)+rho*(3*un.*(1-un)-1/2).*u);
    EM = advection(DM,points);
    advectionPart = @(u) evalOperator(u,R,p,EM);
    JnSSymm = @(u) -alpha*advectionPart(u);
    
    F = @(u) (epsilon*diffusionPart(u)-alpha*advectionPart(u)+...
              rho*u.*(u-1/2).*(1-u));
    
    g = @(u) F(u)-JnSymm(u)-JnSSymm(u);
    
    %time step
    Fn = F(un);
    
    [phi1,errest,info,eigenValues] = phileja(dt,JnSymm,JnSSymm,1,dt*Fn);
    
    Un2 = un+phi1;
    Dn2 = g(Un2)-g(un);
    
    [phi3,errest,info]=phileja(dt,JnSymm,JnSSymm,3,2*dt*Dn2,[],eigenValues);
    
    err = phi3;
    
    u = Un2+err;
    
    normerror = norm(err,inf);
    
    factor = min(2,max(0.6,abs(0.9*(tol/normerror)^(1/2))));
    
    dt = factor*dt;
    
    if normerror> tol
      fprintf('integration error too big %e',normerror);
      dt = -dt;
    end
    u(u<0) = 0;
    % inflow left boundary;
%     u(end-1) = 0;
    
%     u = u(1:end-2);
  end

  function EM = advection(DM,points)
    dx = DifferenceMatrix(points,points);
    EM = rbf_x(DM,shapeParameter,dx);
  end

  function EM = diffusion(DM,points)
    EM = rbf_xx(DM,shapeParameter);
  end

  function EM = differentialOperator(DM,points,epsilon,alpha)
    dx = DifferenceMatrix(points,points);
    EM_x = rbf_x(DM,shapeParameter,dx);
    EM_xx = rbf_xx(DM,shapeParameter);
    
    EM = epsilon*EM_xx-alpha*EM_x;
  end

  function Au = evalOperator(u,R,p,EM)
    coef = zeros(size(u));
    coef(p) = R\(R'\(u(p)));
    Au = EM*coef;
  end

  function [le,re,ri] = eigLNonlinear(diffusion,advection,reaction,N,varargin)
    
    %symmetic part
    
    numEig = 10;
    
    opts.disp = 0;
    opts.tol = 1;
    opts.maxit = 1;
    opts.issym = true;
    opts.isreal = true;
    
    [realeigs]=eigs(@(u) diffusion(u)+reaction(u),N,numEig,'lm',opts);
    [realeigs,index] = sort(realeigs);
    
    le = realeigs(1)
    re = 0;
    
    % skew symmetric part
    
    clear opts;
    
    opts.disp = 0;
    opts.tol = 1;
    opts.maxit = 1;
    opts.issym = true;
    opts.isreal = true;
    
    
    [imageigs] =eigs(@(u) advection(advection(u)),N,numEig,'lm',opts);
    [imageigs,index] = sort(imageigs);
    
    ri = sqrt(abs(imageigs(1)));
  end

%========================================================

  function checkPoints = calculateCheckPoints(points)
    points = [borderPoints(1);points(:);borderPoints(2)];
    checkPoints = (points(2:end)+points(1:end-1))/2;
  end

  function values = evaluateInterpolant(dsites,ctrs,coef)
    DM = distanceMatrix(ctrs,dsites);
    EM = rbf(DM,shapeParameter);
    values = EM*coef;
  end


  function DM = DifferenceMatrix(ctrs,dsites)
    [dr,cc] = ndgrid(dsites,ctrs);
    DM = dr-cc;
  end


  function DM = distanceMatrix(ctrs,dsites)
    DM = abs(DifferenceMatrix(ctrs,dsites));
  end

  function [result,pointer,idx]=unique_tol(v,tol)
    % result: set with unique points
    % pointer: vector that points from the old set to the new set
    % idx: indices of the old point contained in the new set
    DM = distanceMatrix(v,v);
    
    idx = 1:size(v,1);
    i=1;
    pointer = zeros(length(v),1);
    while ~isempty(idx(i+1:end))
      todelete = find(DM(idx(i+1:end),idx(i))<=tol)+i;
      pointer(idx(todelete))=i;
      idx = idx(setdiff(1:length(idx),todelete));
      i = i+1;
    end
    pointer(pointer==0) = 1:length(idx);
    result = v(idx,:);
  end

end

Contact us