Martin J. Gander
Section de Mathématiques
2-4 rue du Lièvre, CP 64
CH-1211 Genève
Martin.Gander(at)unige.ch   Fax: +41 22 379 11 76

Home Research Teaching Consulting  
   
Seminars
Genève:
- Analyse numérique
McGill:
- CSE
- Applied Mathematics
Conference
Ascona 2013

Available projects
Postdoctoral:
- upon request
Doctoral degree:
- upon request
Master thesis:
- upon request

:: Compact implementation of fast and oblivious convolution
Here is a compact implementation in Matlab of the fast and oblivious convolution algorithm by Lubich and Schädle, SIAM J. Sci. Comput., 24:161-182, 2002.
        
function v=Convolution(F,g,dt,N,B,mu0,nu,sigma)
% CONVOLUTION fast convolution algorithm by Lubich and Schaedle
%  v=Convolution(F,g,N,B,mu0,nu,sigma) computes approximations of the
%  convolution integrals int(f(j*dt-tau)*g(t),tau=0..j*dt),
%  j=1,2,...,length(g)-1, where f is defined by its Laplace transform
%  F, g is a vector, and dt is a time step, using the fast algorithm
%  of Lubich and Schaedle, with 2*N+1 quadrature points, base B and
%  paramters mu, nu and sigma. The result is the vector v of the same
%  length as g. 

Nt=length(g)-1;
[la0,w0]=NodesAndWeights(N,mu0/dt,nu,sigma);
Lg=ceil(log(1+Nt*(B-1))/log(B)-10*eps)-1; % number of Talbot contours
tau=zeros(Lg,1);                          % 10*eps to compensate for
for l=1:Lg                                % roundoff errors
  Tl=(2*B^l-1)*dt;
  [la(:,l),w(:,l)]=NodesAndWeights(N,mu0/Tl,nu,sigma);
  Bs(l)=2+sum(B.^(1:l-1));
end;
phi1=sum(w0.*feval(F,la0)./la0.*exp(la0*dt));
phi2=sum(w0.*feval(F,la0)./la0.^2.*exp(la0*dt));

v(1)=phi1*g(1)+phi2*(g(2)-g(1))/dt;
y=zeros(2*N+1,Lg);                        % to store ODE solutions
for n=2:Nt
  t=n*dt;                                 % do the last dt
  for l=1:Lg                              % advance all ODEs
    y(:,l)=y(:,l)+(exp(dt*la(:,l))-1)./(dt*la(:,l))...
           .*(dt*la(:,l).*y(:,l)+dt*g(n-1)+dt*(g(n)-g(n-1))/dt./la(:,l))...
           -dt*(g(n)-g(n-1))/dt./la(:,l);
    if l==1                               % always activate first
      ya(:,l)=y(:,l); tau(l)=tau(l)+dt;
    else
      if mod(n,B^(l-1))==1                % save for future use
        ys(:,l)=y(:,l); y(:,l-1)=zeros(2*N+1,1);
      elseif n>=Bs(l) & mod(n-Bs(l),B^(l-1))==0     % activate
        ya(:,l)=ys(:,l); tau(l)=tau(l)+B^(l-1)*dt;
      end;
    end;
  end;
  v(n)=phi1*g(n)+phi2*(g(n+1)-g(n))/dt;
  L=ceil(log(1+n*(B-1))/log(B)-10*eps)-1; % number of active Talbot contours
  for l=1:L                               % sum remaining contributions
    v(n)=v(n)+sum(w(:,l).*feval(F,la(:,l)).*exp((t-tau(l))*la(:,l)).*ya(:,l));
  end;
end;

function [la,w]=NodesAndWeights(N,mu,nu,sigma);
% NODESANDWEIGHTS computes nodes and weights of Talbot contour
                      
th=(-N:N)*pi/(N+1);              
for j=1:length(th)
  if th(j)~=0
    la(j,1)=sigma+mu*(th(j)*cot(th(j))+i*nu*th(j));
    w(j,1)=1/(2*i*(N+1))*mu*((cos(th(j))*sin(th(j))-th(j))/sin(th(j))^2+i*nu);
  else                                    % need to handle th=0 separately
    la(j,1)=sigma+mu;
    w(j,1)=mu*nu/(2*(N+1));
  end;
end;
       
        
The following commands run a simple example from the original paper:
       
F=inline('(sqrt(pi./s))','s');            
dt=1/100; N=12; B=2;
g=ones(2+sum(B.^(1:3)),1);
mu0=8; nu=0.6; sigma=0;
v=ConvolutionOld(F,g,dt,N,B,mu0,nu,sigma);
       
        
   
Home Research Teaching Consulting Personal