#!/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);