You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
353 lines
11 KiB
353 lines
11 KiB
function [ N, X , T , an , ae, c, nk, re] = expaper( ex )
|
|
|
|
%EXPAPER Numerical examples considered in the paper
|
|
% Input:
|
|
% EX = integer between 1 and 8, corresponding to Example 7.EX
|
|
% Output:
|
|
% N = vector containing the number of nodes and weights in the
|
|
% quadrature formulae
|
|
% X(1,:) = vector containing the exact (maximal) relative error as a
|
|
% function of N for RFEJER with one input argument
|
|
% X(2,:) = vector containing the exact (maximal) relative error as a
|
|
% function of N for RFEJER with two input arguments
|
|
% (the following outputs are not defined for EX == 2 and EX == 7)
|
|
% T(1,:) = vector containing the average time (of 100 measurements)
|
|
% to construct the N-point quadrature formula
|
|
% T(2,:) = vector containing the average time (of 100 measurements)
|
|
% to construct the k-point quadrature formulae for k=c:i:N
|
|
% T(3,:) = vector containing the average time (of 100 measurements)
|
|
% for the semi-automatic integrator
|
|
% [AN,AE,C,Y] = RFEJER(SGL,FUN [,ARRAY]) where SGL and FUN are
|
|
% respectively the sequence of poles and the integrand from
|
|
% Example 7.EX
|
|
% NK = length(Y(1,:))
|
|
% RE = exact relative error on the obtained approximation(s) AN
|
|
%
|
|
% Remarks:
|
|
% * For EX == 1 a choice should be made between OMEG == 1.1 and
|
|
% OMEG == 1.001 by removing/adding `%' on line 36
|
|
% * Possible messages of the form "warning: division by zero" can be
|
|
% avoided by using the alternatives from lines 325-353 instead on
|
|
% lines 38, 52, 76, and 110-113
|
|
|
|
switch ex
|
|
case 1, % Example 7.1
|
|
omeg = 1.1;
|
|
%omeg = 1.001; % remove/add `%' at the beginning of this line
|
|
fun = @(x)(pi.*x./omeg)./sin(pi.*x./omeg);
|
|
% fun = @(x) sinf(x,omeg); % use this instead of the previous line
|
|
% to avoid a warning message of the form
|
|
% "warning: division by zero"
|
|
sgl1=zeros(1,100);
|
|
sgl1(1:2:end) = omeg*[1:1:50];
|
|
sgl1(2:2:end) = -sgl1(1:2:end);
|
|
if (omeg == 1.1),
|
|
exact = 4.4677736463877657892361233985034;
|
|
else % omeg = 1.001
|
|
exact = 12.929256850002296208119060780724;
|
|
end
|
|
case 2, % Example 7.2
|
|
omeg = 1.1;
|
|
fun = @(x)(pi.*x./omeg)./sin(pi.*x./omeg);
|
|
% fun = @(x) sinf(x,omeg); % use this instead of the previous line
|
|
% to avoid warning messages of the form
|
|
% "warning: division by zero"
|
|
sgl1=zeros(5,100);
|
|
sgl1(1,1:2:end) = omeg*[1:1:50]; % optimal set
|
|
sgl1(1,2:2:end) = -sgl1(1,1:2:end);
|
|
|
|
del = 1e-10; % small perturbations
|
|
sgl1(2,:) = (1+del)*sgl1(1,:); % set 1
|
|
sgl1(3,:) = (1-del)*sgl1(1,:); % set 2
|
|
|
|
rho = omeg;
|
|
sgl1(4,1:2:end) = rho; % set 3
|
|
sgl1(4,2:2:end) = -rho;
|
|
|
|
rho = 1.3;
|
|
sgl1(5,1:2:end) = rho; % set 4
|
|
sgl1(5,2:2:end) = -rho;
|
|
|
|
sgl1(6,:) = inf; % set 5
|
|
exact = 4.4677736463877657892361233985034;
|
|
case 3, % Example 7.3
|
|
omeg = 0.1;
|
|
fun = @(x)(pi.*x./omeg)./sinh(pi.*x./omeg);
|
|
% fun = @(x) sihf(x,omeg); % use this instead of the previous line
|
|
% to avoid warning messages of the form
|
|
% "warning: division by zero"
|
|
sgl1=zeros(1,100);
|
|
sgl1(1:2:end) = -sqrt(-1)*omeg*[1:1:50];
|
|
sgl1(2:2:end) = -sgl1(1:2:end);
|
|
exact = 0.15707963267939592622470892489958;
|
|
case 4, % Example 7.4
|
|
omeg = 1.1;
|
|
fun = @(x) sin(1./(x-omeg));
|
|
sgl1 = omeg*ones(1,100);
|
|
exact = -1.19245706732219214078770039225;
|
|
case 5, % Example 7.5
|
|
omeg = 0.1;
|
|
fun = @(x) 1./(exp(pi.*x./omeg)+1);
|
|
sgl1 = zeros(1,100);
|
|
sgl1(1:2:end) = sqrt(-1)*omeg*[1:2:100];
|
|
sgl1(2:2:end) = -sgl1(1:2:end);
|
|
exact = 1;
|
|
case 6, % Example 7.6
|
|
fun = @(x) 1./sqrt((x+3).*(x+2));
|
|
sgl1 = -2.5*ones(1,100);
|
|
exact = 0.87116861981054736678487404465635;
|
|
case 7, % Example 7.7
|
|
fun = @(x) exp(x)./(25*x.^2+1);
|
|
sgl1 = [.2*sqrt(-1)*[1,-1],inf*ones(1,98)];
|
|
sgl2 = inf*ones(1,100);
|
|
exact = 0.579408643011693305029559997358;
|
|
otherwise % Example 7.8
|
|
omeg = [2, 1.1, 1.01, 1.001];
|
|
fun = @(x) [(pi.*x./omeg(1))./sin(pi.*x./omeg(1)), ...
|
|
(pi.*x./omeg(2))./sin(pi.*x./omeg(2)); ...
|
|
(pi.*x./omeg(3))./sin(pi.*x./omeg(3)), ...
|
|
(pi.*x./omeg(4))./sin(pi.*x./omeg(4))];
|
|
% fun = @(x) [sinf(x,omeg(1)), ... % use this instead of the
|
|
% sinf(x,omeg(2)); ... % previous lines to avoid
|
|
% sinf(x,omeg(3)), ... % warning messages of the form
|
|
% sinf(x,omeg(4))]; % "warning: division by zero"
|
|
sgl1 = zeros(1,102);
|
|
sgl1(1:6:end) = -omeg(4)*[1:1:17];
|
|
sgl1(3:6:end) = -omeg(3)*[1:1:17];
|
|
sgl1(5:6:end) = -omeg(2)*[1:1:17];
|
|
sgl1(2:2:end) = -sgl1(1:2:end);
|
|
exact = [2.33248723224655024110707565174 , ...
|
|
4.4677736463877657892361233985034; ...
|
|
8.43018458047084205897126420479 , ...
|
|
12.929256850002296208119060780724];
|
|
end
|
|
|
|
N = [2:1:41];
|
|
n = length(N);
|
|
T = zeros(3,n);
|
|
X = zeros(2,n);
|
|
|
|
if (ex < 7) && ~(ex == 2), % Examples 7.1, 7.3, 7.4, 7.5, and 7.6
|
|
|
|
% exact relative errors and execution times
|
|
|
|
for M = 1:1:n,
|
|
disp(sprintf('Execution time measurements countdown: %1.3f', n-M+1));
|
|
sgl = sgl1(1:N(M)-1);
|
|
for K=1:100,
|
|
tic;
|
|
[x,Aerr] = rfejer(sgl);
|
|
A=fun(x(1,:))*x(2,:)';
|
|
t1(K) = toc;
|
|
tic;
|
|
[an1,ae1,c1,y1] = rfejer(sgl,fun);
|
|
t2(K) = toc;
|
|
end
|
|
T(1,M) = mean(t1);
|
|
T(3,M) = mean(t2);
|
|
X(1,M) = abs(exact-A)/abs(exact);
|
|
X(2,M) = abs(exact-an1)/abs(exact);
|
|
|
|
if M == 1,
|
|
T(2,1) = T(1,1);
|
|
elseif (ex == 7), % c = 3, i = 1
|
|
T(2,M) = sum(T(1,2:M));
|
|
elseif (ex == 3) || (ex == 5), % c = 3, i = 2
|
|
T(2,M) = sum(T(1,2:2:M));
|
|
else
|
|
T(2,M) = sum(T(1,1:M));
|
|
end
|
|
end
|
|
index = find(X==0);
|
|
X(index)=eps/10;
|
|
|
|
% convergence / deterioration
|
|
|
|
[an,ae,c,y] = rfejer(sgl1,fun);
|
|
nk = length(y(1,:));
|
|
re = abs(exact-an)/abs(exact);
|
|
|
|
elseif (ex == 8), % Example 7.8
|
|
|
|
% exact relative errors and execution times
|
|
|
|
for M = 1:1:n,
|
|
disp(sprintf('Execution time measurements countdown: %1.3f', n-M+1));
|
|
sgl = sgl1(1:N(M)-1);
|
|
for K=1:100,
|
|
tic;
|
|
[x,Aerr] = rfejer(sgl);
|
|
F = fun(x(1,:));
|
|
m = length(x(1,:));
|
|
for j = 1:2,
|
|
A(:,j) = F(:,(j-1)*m+1:j*m)*x(2,:)';
|
|
end
|
|
t1(K) = toc;
|
|
tic;
|
|
[an1,ae1,c1,y1] = rfejer(sgl,fun,'Array',true);
|
|
t2(K) = toc;
|
|
end
|
|
T(1,M) = mean(t1);
|
|
T(2,M) = sum(T(1,1:M));
|
|
T(3,M) = mean(t2);
|
|
Xa = abs(exact-A)./abs(exact);
|
|
X(1,M) = max(Xa(:));
|
|
Xa = abs(exact-an1)./abs(exact);
|
|
X(2,M) = max(Xa(:));
|
|
end
|
|
index = find(X==0);
|
|
X(index)=eps/10;
|
|
|
|
% convergence / deterioration
|
|
|
|
[an,ae,c,y] = rfejer(sgl1,fun,'Array',true);
|
|
nk = length(y(1,:));
|
|
re = abs(exact-an)./abs(exact);
|
|
|
|
elseif (ex == 2), % Example 7.2
|
|
|
|
X = inf*ones(5,n);
|
|
T = NaN; an = NaN; ae = NaN; c = NaN; nk = NaN; re = NaN;
|
|
|
|
% n-point rational Fejer
|
|
|
|
for M = 1:1:n,
|
|
% optimal set
|
|
sgl = sgl1(1,1:N(M)-1);
|
|
[x,Aerr] = rfejer(sgl);
|
|
A = fun(x(1,:))*x(2,:)';
|
|
X(1,M) = abs(exact-A)/abs(exact);
|
|
|
|
% set 1
|
|
sgl = sgl1(2,1:N(M)-1);
|
|
[x,Aerr] = rfejer(sgl);
|
|
A = fun(x(1,:))*x(2,:)';
|
|
X(2,M) = abs(exact-A)/abs(exact);
|
|
|
|
% set 2
|
|
sgl = sgl1(3,1:N(M)-1);
|
|
[x,Aerr] = rfejer(sgl);
|
|
A = fun(x(1,:))*x(2,:)';
|
|
X(3,M) = abs(exact-A)/abs(exact);
|
|
|
|
% set 3
|
|
sgl = sgl1(4,1:N(M)-1);
|
|
[x,Aerr] = rfejer(sgl);
|
|
A = fun(x(1,:))*x(2,:)';
|
|
X(4,M) = abs(exact-A)/abs(exact);
|
|
|
|
% set 4
|
|
sgl = sgl1(5,1:N(M)-1);
|
|
[x,Aerr] = rfejer(sgl);
|
|
A = fun(x(1,:))*x(2,:)';
|
|
X(5,M) = abs(exact-A)/abs(exact);
|
|
|
|
% set 5
|
|
sgl = sgl1(6,1:N(M)-1);
|
|
[x,Aerr] = rfejer(sgl);
|
|
A = fun(x(1,:))*x(2,:)';
|
|
X(6,M) = abs(exact-A)/abs(exact);
|
|
|
|
end
|
|
index = find(X==0);
|
|
X(index)=eps/10;
|
|
X1 = X(1,1:19); X2 = X(2,1:19); X3 = X(3,1:19);
|
|
N1 = N(1:19);
|
|
X4 = X(4,:); X5 = X(5,:); X6 = X(6,:);
|
|
figure;
|
|
semilogy(N1,X1,'-',N1,X2,'--',N1,X3,'-.',N,X4,'-+',N,X5,'--o',N,X6,'-o');
|
|
xlabel('n');
|
|
ylabel('relative error');
|
|
hleg1 = legend('opt','Delta = 1e-10','Delta = - 1e-10','rho = omega',...
|
|
'rho = 1.3', 'infinity');
|
|
else % Example 7.7
|
|
X = inf*ones(2,n);
|
|
T = NaN; an = NaN; ae = NaN; c = NaN; nk = NaN; re = NaN;
|
|
|
|
% n-point rational Fejer
|
|
|
|
for M = 1:1:n,
|
|
% rational
|
|
sgl = sgl1(1:N(M)-1);
|
|
[x,Aerr] = rfejer(sgl);
|
|
A = fun(x(1,:))*x(2,:)';
|
|
X(1,M) = abs(exact-A)/abs(exact);
|
|
|
|
% polynomial
|
|
sgl = sgl2(1:N(M)-1);
|
|
[x,Aerr] = rfejer(sgl);
|
|
A = fun(x(1,:))*x(2,:)';
|
|
X(2,M) = abs(exact-A)/abs(exact);
|
|
end
|
|
index = find(X==0);
|
|
X(index)=eps/10;
|
|
X1 = X(1,2:end); N1=N(2:end);
|
|
X2 = X(2,:);
|
|
figure;
|
|
semilogy(N1,X1,'-',N,X2,'--');
|
|
xlabel('n');
|
|
ylabel('relative error');
|
|
hleg1 = legend('rational','polynomial');
|
|
end
|
|
|
|
if ~(ex==2) && ~(ex==7),
|
|
|
|
if ~(ex==3) && ~(ex==5),
|
|
p = [1:1:n];
|
|
else
|
|
p = [2:2:n];
|
|
end
|
|
|
|
figure;
|
|
semilogy(N(p),X(1,p),'-',N(p),X(2,p),'--');
|
|
xlabel('n');
|
|
ylabel('Relative error');
|
|
title('Relative error against number of points');
|
|
|
|
figure;
|
|
loglog(T(3,p),X(2,p));
|
|
xlabel('t');
|
|
ylabel('Relative error');
|
|
title('Relative error against execution time');
|
|
|
|
T1 = T(1,p); T2 = T(2,p); T3 = T(3,p); N1=N(p);
|
|
figure;
|
|
semilogy(N1,T3,'-',N1,T2,'--',N1,T1,'-.');
|
|
xlabel('n');
|
|
ylabel('t');
|
|
hleg1 = legend('t_{n}^{(sa)}','t_{c:i:n}','t_n');
|
|
title('Execution time against number of points');
|
|
|
|
end
|
|
|
|
end
|
|
|
|
%--------------------------------------------------------------------------
|
|
|
|
function [ f ] = sinf( x , w )
|
|
%SINF @(x)(pi.*x./omeg)./sin(pi.*x./omeg)
|
|
|
|
warning('off'); % division by zero is allowed in the next line
|
|
f = (pi.*x./w)./sin(pi.*x./w);
|
|
warning('on');
|
|
index = find(x==0);
|
|
if ~isempty(index),
|
|
f(index)=1;
|
|
end
|
|
|
|
end
|
|
|
|
%--------------------------------------------------------------------------
|
|
|
|
function [ f ] = sihf( x , w )
|
|
%SINF @(x)(pi.*x./omeg)./sinh(pi.*x./omeg)
|
|
|
|
warning('off'); % division by zero is allowed in the next line
|
|
f = (pi.*x./w)./sinh(pi.*x./w);
|
|
warning('on');
|
|
index = find(x==0);
|
|
if ~isempty(index),
|
|
f(index)=1;
|
|
end
|
|
|
|
end
|
|
|