#!/usr/bin/octave
% This program, called EULER.m, has been published by Abate and Whitt
% in "The Fourier series method for inverting transforms of
% probability distributions", Queueing Systems, 10, pp. 5-88, and in
% "Numerical Inversion of Laplace Transforms of Probability
% Distributions", ORSA Journal on Computing 7 (1995), pp 36-43. A&W
% wrote it originally in Ubasic. This is a port to Octave, a GNU version
% of Matlab, written by Nicky van Foreest, foreest@math.utwente.nl
% Date: 08/21/01
% I tried to stay as close as possible to the published code to
% facilitate the reading of this program. Doubtless there will be
% numerous improvements to make by exploiting specific functionality
% of Octave. I did not opt to include these in the present
% program (save the definition of the vector C---I was too lazy to type
% in the numbers; and the definition of the variable M). Easy of
% comparison with the published version had my highest priority.
M=11;
SU=zeros(1,M+2);
C=zeros(1,M+1);
C(1)=1;
for k=[2:M+1]
C(k)=(M-(k-2))*C(k-1)/(k-1);
endfor
T=input("TIME=");
A=19.1;
Ntr=15;
U=exp(A/2.)/T;
X=A/(2.*T);
H=pi/T;
% The function fnRf has to be defined before we can invoke it
function Rfs=fnRf(X,Y)
S=X+i*Y;
Rho=0.75;
Mean=1.;
Gs=1./sqrt(1.+2.*S);
Gse=(1.-Gs)/(Mean*S);
Fs=(1.-Gse)/(S*(1.-Rho*Gse));
Rfs=real(Fs);
endfunction
Sum=fnRf(X,0)/2.;
for N=[1:Ntr]
Y=N*H;
Sum=Sum+(-1)^N*fnRf(X,Y);
endfor
SU(1)=Sum;
for K=[1:M+1]
N=Ntr+K;
Y=N*H;
SU(K+1)=SU(K)+(-1)^N*fnRf(X,Y);
endfor
Avgsu=0.;
Avgsu1=0.;
for J=[1:M+1]
Avgsu=Avgsu+C(J)*SU(J);
Avgsu1=Avgsu1+C(J)*SU(J+1);
endfor
Fun=U*Avgsu/2048;
Fun1=U*Avgsu1/2048;
Errt=abs(Fun-Fun1)/2.;
printf("\nTime=%f\tFUNCTION=%10.9f\n",T,Fun1);
printf("Truncation Error Estimate=%e\n",Errt);