Martin J. Gander Section de Mathématiques Université de Genève Rue du Conseil-Général 9, CP 64 1211 Genève 4, Suisse 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