Code covered by the BSD License

# meshfree exponential integrator

### Stefan Rainer (view profile)

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);
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,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```