PHSTUD8%DKNKURZ1.BITNET@cunyvm.cuny.edu (05/23/89)
Date: 22 May 1989, 19:49:37 MEZ
From: Nikolaj Schutzbach (0049-(0)7531)56464, PHSTUD8 at DKNKURZ1
Hinterhauser Str. 13, D-7750 Konstanz
P720, Universitaet Konstanz, D-7750 Konstanz
To: info-pascal-request at vim.brl
I'm searching desperatly for a function which integrates
(any) function from zero or a point larger than that to infinity.
Every language is welcome, as long as I get source code, so I may
translate it to one of my favourite languages.
Thanks in advance
Niko.ts@chyde.uwasa.fi (Timo Salmi LASK) (05/23/89)
In article <19715@adm.BRL.MIL> PHSTUD8%DKNKURZ1.BITNET@cunyvm.cuny.edu writes: >I'm searching desperatly for a function which integrates >(any) function from zero or a point larger than that to infinity. >Every language is welcome, as long as I get source code, so I may >translate it to one of my favourite languages. Try the following for an interim solution until you get sources for what you are looking for from somewhere. Ftp anonymous pc/ts/tsnum12.arc from 128.214.12.3. This package includes interation of general functions between any two values. The values must be numerical, though. The option of true infinite is not included. - Also (as is evident) Turbo Pascal Numerical Toolbox might be of some use to you if you do not already have it. ................................................................... Prof. Timo Salmi School of Business Studies, University of Vaasa, SF-65101, Finland Internet: ts@chyde.uwasa.fi Funet: vakk::salmi Bitnet: salmi@finfun
dmurdoch@watstat.waterloo.edu (Duncan Murdoch) (05/24/89)
In article <19715@adm.BRL.MIL> PHSTUD8%DKNKURZ1.BITNET@cunyvm.cuny.edu writes: >I'm searching desperatly for a function which integrates >(any) function from zero or a point larger than that to infinity. Numerical Recipes, by Press et al, suggests using the identity b 1/a 1 1 I f(x) dx = I - f(-) dt a 1/b t^2 t which works provided a and b have the same sign, to turn an upper limit of infinity into a lower limit of zero. If a=0, you'll have to break the integral into two parts, i.e. 0 to c, c to infinity, and sum the results. To integrate over a finite range, I've had very good luck with the routine QUANC8, from Forsythe, Malcolm and Moler, Computer Methods for Mathematical Computations. Here's a copy, translated to Turbo Pascal v 5. Duncan Murdoch ---- cut here ----- FMM.PAS ------------- unit fmm; interface type float=double; integrand = function(var x;info:pointer):float; procedure quanc8(fun:integrand;info:pointer;a,b,abserr,relerr:float; var result,errest:float; var nofun:integer;var flag:float); { Subroutine from Forsythe, Malcolm and Moler (1977) Computer Methods for Mathematical Computations Chapter 5 to estimate the integral of fun(x) from a to b to a user provided tolerance. (Translated from Fortran to Ratfor and from Ratfor to Turbo Pascal by djm.) An automatic adaptive routine based on the 8-panel Newton-Cotes rule is used. Input: fun: The name of the integrand function subprogram fun(x) info: A pointer which will be passed to fun a: The lower limit of integration. b: " upper " " " . (b may be less than a) relerr: A relative error tolerance (non-negative) abserr: An absolute " " " Output: result: An approximation to the integral hopefully satisfying the least stringent of the two error tolerances errest: An estimate of the magnitude of the actual error nofun: The number of function evaluations flag: A reliability indicator. If flag is zero, then result probably satisfies the error tolerance. If flag is XXX.YYY, then XXX = the number of intervals which have not converged and 0.YYY = the fraction of the interval (a,b) left to do when the limit on nofun was approached } implementation function max(x,y:float):float; begin if x>=y then max := x else max := y; end; procedure quanc8; var w0,w1,w2,w3,w4,area,x0,f0,stone,step,cor11,temp, qprev,qnow,qdiff,qleft,esterr,tolerr:float; qright: array [1..31] of float; f,x: array[1..16] of float; fsave,xsave: array[1..8,1..30] of float; levmin,levmax,levout,nomax,nofin,lev,i,j:integer; nim:longint; label break; begin { Stage 1 General Initialization } { Set constants } levmin := 1; levmax := 30; levout := 6; nomax := 5000; nofin := nomax - 8*(levmax - levout + 1 shl (levout+1)); { Trouble when nofun reaches nofin } w0 := 3956.0E0/14175.0E0; w1 := 23552.0E0/14175.0E0; w2 := -3712.0E0/14175.0E0; w3 := 41984.0E0/14175.0E0; w4 := -18160.0E0/14175.0E0; { Initialize running sums to zero. } flag := 0.0E0; result := 0.0E0; cor11 := 0.0E0; errest := 0.0E0; area := 0.0E0; nofun := 0; if (a = b) then exit; { Stage 2 Initialization for first interval } lev := 0; nim := 1; x0 := a; x[16] := b; qprev := 0.0; f0 := fun(x0,info); stone := (b-a)/16.0E0; x[ 8] := ( x0 + x[16])/2.0E0; x[ 4] := ( x0 + x[ 8])/2.0E0; x[12] := (x[ 8] + x[16])/2.0E0; x[ 2] := ( x0 + x[ 4])/2.0E0; x[ 6] := (x[ 4] + x[ 8])/2.0E0; x[10] := (x[ 8] + x[12])/2.0E0; x[14] := (x[12] + x[16])/2.0E0; for j:=1 to 8 do f[2*j] := fun(x[2*j],info); nofun := 9; repeat { Stage 3 Central Calculation } { Requires qprev,x0,x2,...,x16,f0,f2,...,f16. } { Calculates x1,x3,...,x15, f1,f3,...,f15,qleft,qright,qnow } { qdiff,area. } x[1] := (x0 + x[2])/2.0E0; f[1] := fun(x[1],info); for j:=2 to 8 do {3 to 15 by 2} begin x[2*j-1] := (x[2*j-2] + x[2*j])/2.0E0; f[2*j-1] := fun(x[2*j-1],info); end; nofun := nofun + 8; step := (x[16] - x0)/16.0E0; qleft := (w0*(f0+f[8]) + w1*(f[1]+f[7]) + w2*(f[2]+f[6]) + w3*(f[3]+f[5]) + w4*f[4])*step; qright[lev+1] := (w0*(f[8]+f[16]) + w1*(f[9]+f[15]) + w2*(f[10]+f[14]) + w3*(f[11]+f[13]) + w4*f[12])*step; qnow := qleft + qright[lev+1]; qdiff := qnow - qprev; area := area + qdiff; { Stage 4 Interval Convergence test } esterr := abs(qdiff)/1023.0E0; tolerr := max(abserr,relerr*abs(area))*(step/stone); if (lev < levmin) or ((lev < levmax) and (nofun <= nofin) and (esterr > tolerr)) then begin { Stage 5 No convergence } { Locate next interval. } nim := 2*nim; lev := lev + 1; { Store right hand elements for future use } for i:=1 to 8 do begin fsave[i,lev] := f[i+8]; xsave[i,lev] := x[i+8]; end; { Assemble left hand elements for immediate use } qprev := qleft; for i:=1 to 8 do begin j := -i; f[2*j+18] := f[j+9]; x[2*j+18] := x[j+9]; end; { Go on to next iteration } end else begin if (nofun > nofin) then begin { Stage 6 Trouble Section } { Number of function evaluations is about to exceed } { limit } nofin := 2*nofin; levmax := levout; flag := flag + (b-x0)/(b-a); end else if (lev >= levmax) then begin { Current level is levmax } flag := flag + 1.0; end; { Stage 7 Interval converged } { Add contributions into running sums } result := result + qnow; errest := errest + esterr; cor11 := cor11 + qdiff/1023.0E0; { Locate next interval } while (nim <> 2*(nim div 2)) do begin nim := nim div 2; lev := lev - 1; end; nim := nim + 1; if (lev <= 0) then goto break; { Completely done } { Assemble elements required for the next interval } qprev := qright[lev]; x0 := x[16]; f0 := f[16]; for i:=1 to 8 do begin f[2*i] := fsave[i,lev]; x[2*i] := xsave[i,lev]; end; end; until false; { Stage 8 Finalize and return } break: result := result + cor11; { Make sure errest not less than roundoff level } if (errest = 0.0E0) then exit; temp := abs(result) + errest; while (temp = abs(result)) do begin errest := 2.0E0*errest; temp := abs(result) + errest; end; end; end. --------- end of FMM.PAS -----------
nmm@cl.cam.ac.uk (Nick Maclaren) (05/25/89)
dmurdoch@watstat.waterloo.edu (Duncan Murdoch) writes: >Numerical Recipes, by Press et al, suggests using the identity > b 1/a 1 1 > I f(x) dx = I - f(-) dt > a 1/b t^2 t >which works provided a and b have the same sign, to turn an upper limit of >infinity into a lower limit of zero. If a=0, you'll have to break the integral >into two parts, i.e. 0 to c, c to infinity, and sum the results. This sort of thing is the reason that 'cookbooks for the naive' are so dangerous, and it is much better to get access to a good (and hence probably expensive) numerical library or consult a numerical analyst. Yes, this transformation may work and may even give the right answer. Alternatively, it may just APPEAR to work, and give an answer that is almost unrelated to the correct one. For example, it is extremely common to get functions of the general nature of f(x)exp(-ax)sin(bx) or f(x)sin(ax)/(x^n). If you apply this transformation to such functions, you end up with a function that is very nasty indeed (as far as numerical integration is concerned). Many free numerical integration routines will happily return an answer, but it may be nonsense. This area is a minefield, and should be treated as such. It is critical to work out at least the gross properties of a function before trying to integrate it, especially over an infinite range. Nick Maclaren University of Cambridge Computer Laboratory nmm @ uk.ac.cam.cl