Contents

Sparse Grid Demos

Author(s): Paul R. Miles & Ralph C. Smith

Date Created: April 10, 2019

Exampled created using MATLAB2018a. Certain features may not be compatible with previous versions of MATLAB.

The sparsegrid toolbox is available for download on Andreas Klimke's ResearchGate page: https://www.researchgate.net/publication/330956911_Sparse_Grid_Interpolation_Toolbox_v511

% Setup workspace
close all; clear; clc;

% Add path to sparse grid toolbox
% Note: I assume the uncompressed folder has been placed in the current
% working directory.
addpath('spinterp_v5.1.1');

Example 1

$f(x,y) = \sin(x) + \cos(y)$ Domain: $[0,\pi]\times[0,\pi]$

% define function
f = @(x,y) sin(x) + cos(y);
% construct interpolant
spi = spvals(f, 2, [0 pi; 0 pi]);

Compute interpolated values - random points in domain

x1 = pi*rand(1, 5); x2 = pi*rand(1, 5);
y = spinterp(spi, x1, x2);
error = max(abs(y - f(x1,x2)));
% display comparison between interpolated and exact
fprintf('(%4s,%4s) -> %6s, %6s\n', 'x1', 'x2', 'y(intrp)', 'f(x1,x2)');
for ii = 1:5
    fprintf('(%4.2f,%4.2f) -> %+6.4f, %+6.4f\n', x1(ii), x2(ii), y(ii), f(x1(ii),x2(ii)));
end
fprintf('max(abs(y-f(x1,x2))) = %4.2e\n', error);
(  x1,  x2) -> y(intrp), f(x1,x2)
(0.34,2.73) -> -0.5869, -0.5873
(3.02,0.27) -> +1.0797, +1.0844
(0.01,1.26) -> +0.3228, +0.3242
(2.43,0.82) -> +1.3299, +1.3345
(2.57,2.51) -> -0.2644, -0.2662
max(abs(y-f(x1,x2))) = 4.69e-03

Visualize grid

figure(1);
h = plotgrid(spi.maxLevel,spi.d);
pcol = distribute_color_spectrum(length(h), jet);
for ii = 1:length(h)
    set(h(ii), 'Marker','o');
    set(h(ii), 'MarkerSize', 10);
    set(h(ii), 'LineWidth', 2);
    set(h(ii), 'Color', pcol(ii,:));
end
xlabel('X')
ylabel('Y')
set(gca, 'FontSize', 22);
% set_figure_dimensions(8,6)

Visualize the interpolant

figure(2)
subplot(1,2,1);
ezmesh(f,[0,pi]);
title('f(x,y)=sin(x) + cos(y)');
subplot(1,2,2);
ezmesh(@(x,y) spinterp(spi,x,y), [0 pi]);
title('Sparse grid interpolant');
set_figure_dimensions(10,6);

Compute integral

% using spquad
ns = 100; % perform integration ns times
tic
for ii = 1:ns
    Qsg = spquad(spi);
end
sgtime = toc;
% using MATLAB's integral2
tic
for ii = 1:ns
    Q2 = integral2(f, 0, pi, 0, pi);
end
i2time = toc;
% display comparison
fprintf('Qsg = %6.4f, Q2 = %6.4f\n', Qsg, Q2);
fprintf('Avg. Comp. Time (sec):\n\tQsg -> %4.3e\n\tQ2 -> %4.3e\n', sgtime/ns, i2time/ns);
Qsg = 6.2630, Q2 = 6.2832
Avg. Comp. Time (sec):
	Qsg -> 5.137e-04
	Q2 -> 7.105e-04

Example 2

peaks($x, y$) Domain: $[0,2]\times[0,2]$

% define function
f = @(x,y) peaks(x,y);
% define sparse grid options
options = spset('DimensionAdaptive', 'on', 'RelTol', 1e-4, 'GridType', 'Chebyshev');
spi = spvals(f, 2, [0,2; 0,2], options);

Visualize grid

figure(3);
h = plotgrid(spi.maxLevel,spi.d);
tmpcol = jet;
pcol = distribute_color_spectrum(2*length(h), tmpcol(1:end,:));
for ii = 1:length(h)
    set(h(ii), 'Marker','o');
    set(h(ii), 'MarkerSize', 6);
    set(h(ii), 'LineWidth', 2);
    set(h(ii), 'Color', pcol(2*ii,:));
end
xlabel('X')
ylabel('Y')
cbar=colorbar;
caxis([1, length(h)])
set(cbar, 'YTick', 1:1:length(h))
set(gca, 'FontSize', 22);
set_figure_dimensions(8,6)

Compute integral

% using spquad
ns = 100; % perform integration ns times
tic
for ii = 1:ns
    Qsg = spquad(spi);
end
sgtime = toc;
% using MATLAB's integral2
tic
for ii = 1:ns
    Q2 = integral2(f, 0, 2, 0, 2);
end
i2time = toc;
% display comparison
fprintf('Qsg = %6.4f, Q2 = %6.4f\n', Qsg, Q2);
fprintf('Avg. Comp. Time (sec):\n\tQsg -> %4.3e\n\tQ2 -> %4.3e\n',sgtime/ns, i2time/ns);
Qsg = 9.9553, Q2 = 9.9553
Avg. Comp. Time (sec):
	Qsg -> 6.463e-04
	Q2 -> 9.383e-04

Example 3

$f = \sin(\pi x)\cdot\sin(\pi y)\cdot2z$ Domain: $[0,1]\times[0,1]\times[0,1]$

% define function
f = @(x,y,z) sin(pi*x).*sin(pi*y).*2.*z;
% construct interpolant
options = spset('DimensionAdaptive', 'on', 'RelTol', 1e-4, 'GridType', 'Chebyshev');
spi = spvals(f, 3, [0,1; 0,1; 0,1], options);

Visualize grid

figure(4);
h = plotgrid(spi.maxLevel,spi.d);
tmpcol = jet;
pcol = distribute_color_spectrum(2*length(h), tmpcol(1:end,:));
for ii = 1:length(h)
    set(h(ii), 'Marker','o');
    set(h(ii), 'MarkerSize', 6);
    set(h(ii), 'LineWidth', 2);
    set(h(ii), 'Color', pcol(2*ii,:));
end
xlabel('X')
ylabel('Y')
zlabel('Z')
cbar=colorbar;
caxis([1, length(h)])
set(cbar, 'YTick', 1:1:length(h))
set(gca, 'FontSize', 22);
set_figure_dimensions(8,6)

Compute integral

% using spquad
ns = 100; % perform integration ns times
tic
for ii = 1:ns
    Qsg = spquad(spi);
end
sgtime = toc;

% using MATLAB's integral3
tic
for ii = 1:ns
    Q3 = integral3(f, 0, 1, 0, 1, 0, 1);
end
i3time = toc;
% display comparison
fprintf('Qsg = %6.4f, Q3 = %6.4f\n', Qsg, Q3);
fprintf('Avg. Comp. Time (sec):\n\tQsg -> %4.3e\n\tQ3 -> %4.3e\n',sgtime/ns, i3time/ns);
Qsg = 0.4053, Q3 = 0.4053
Avg. Comp. Time (sec):
	Qsg -> 7.263e-04
	Q3 -> 3.017e-02

Support Functions

function [psty] = distribute_color_spectrum(nds, cmap)
% evenly distribution color palette across nds data sets
if nargin < 2 || isempty(cmap)
    cmap = jet; % default color map
end
cmp = colormap(cmap);
[mc,~] = size(cmp);
nc = linspace(1,mc,nds)';
psty = zeros(nds,3);
for ii = 1:nds
    psty(ii,:) = interp1((1:1:mc)', cmp, nc(ii));
end
end

function set_figure_dimensions(width, height)
% adjust figure dimensions
set(gcf,'units','inches');
pos = get(gca, 'pos');
set(gcf, 'Position', [pos(1) pos(2) width height]);
end