\ Exponential Integral, fE1(x), and Ivantsov function, fIv(x).
\ --
\ Reference for Exponential Integral:
\ [1] Walter Gautschi and William F. Cahill, "Exponential Integral and
\ Related Functions," in Handbook of Mathematical Functions,
\ Miltion Abrmowitz and Stegun, eds., Dover Publications
\ New York, 1965, pp. 227-254.
\ --
\ References for Ivantsov function:
\ [2] GP Ivantsov, Doklady Akademmii Nauk SSSR [58] (1947) 567
\ [3] W. Kurz, D.J. Fisher, "Fundamental of Solidification, 4th Ed,"
\ Trans Tech Publications, 1998, pp. 242-244.
\ --
\ Notes:
\ Two versions for the Exponential Integral, E1(x), are given. One is based
\ on a Taylor series expansion and converges faster as x<1, while the other
\ is based on a continued fraction is converges faster as x>1.
\ These are then combined into one Forth word. Clearly, other methods
\ such as approximate polynomials, etc., would be faster and more efficient;
\ however, that was not the point of this particular programming exercise.
\ The code should work with both the zenfloat and ansfloat 4tH packages
\ as well as an ansforth system. David Johnson 6/8/2009
\ --
\ ** For an ansforth system, the following are needed: **
\ include easy.4th ( which defines words: FLOAT, S>F, S>FLOAT, ARRAY )
\ include taylor.4th ( from the 4tH library )
true constant ZENFLOAT?
ZENFLOAT? [IF] \ Use the ZENFLOAT package
include lib/zenfloat.4th
include lib/zentaylr.4th
include lib/zenfexp.4th
include lib/zenfln.4th
include lib/zenans.4th
[THEN]
ZENFLOAT? [NOT] [IF] \ Use the ANSFLOAT package
include lib/ansfloat.4th
include lib/ansfpio.4th
include lib/taylor.4th
include lib/flnflog.4th
include lib/fexp.4th
fclear 9 set-precision
[THEN]
\ --Exponential Integral--
\ E1(x) = Integral from x to infinity of (exp(-t)/t)dt, or
\ E1(x) = -Ei(-x) (ref[1], p. 228) with x being real and x>0.
\ (i) As a Taylor series expansion of N terms.
\ E1(x) = -g - ln(x)-series{((-1)^N)(x^N))/(N*N!)}
\ with g=0.5772156649... (Euler's constant)
: E1Low ( f1 N -- f2 )
>r fdup fln s" 577215665e-9" s>float f+ fswap
fnegate fdup fdup \ setup for Taylor series
1 r> 2 Do \ initial of value N! is 1
i * dup i * swap >r +taylor r> \ N! and i*N! for series
loop
drop fdrop fdrop \ clean up Taylor series calculation
f+ fnegate \ finish the approximation
;
\ (ii) As an Nth continued fraction.
\ E1(x) = exp(-x)(1/(x+),1/(1+), 1/(x+), 2/(1+), 2/(x+), ...)
FLOAT array px \ variable use for portability
: (E1high) ( r1 N -- r2 )
>r fdup px f!
r@ 1+ s>f f+
0 r> do
I 0= if leave then \ for consistent behavior in 4tH and ansforth
I s>f fswap f/ 1 s>f f+
I s>f fswap f/ px f@ f+
-1 +loop
1 s>f fswap f/
;
: E1High ( f1 N -- f2) >r fdup r> (E1high) fswap fexp f/ ;
: fE1 ( f1 -- f2 ) \ Expontial Integral, E1(x)
fdup f0= abort" fE1: divide by zero"
fdup f0< abort" fE1: negative float "
fdup s" 67e-2" s>float f< if 32 E1Low else 40 E1High then ;
\ --Invantsov Function--
\ fIv(x)= x*Exp(x)*E1(x)
: fIv ( f1 -- f2)
fdup f0= abort" fIv: divide by zero"
fdup f0< abort" fIv: negative float"
fdup fdup s" 67e-2" s>float f<
if 32 E1Low fswap fdup fexp f* f*
else 40 (E1High) f* then
;
\ ----------Testing ---------------------
\ Test the Invantsov function: (note behavior relative to erf(x)).
true [if] \ ref [3] limited precision,see next Test.
cr ." x fIv(x) Compared to:"
s" 1e-6" s>float fdup cr f. space fIv f. ." -----"
s" 1e-2" s>float fdup cr f. space fIv f. ." 0.04079"
s" 1e-1" s>float fdup cr f. space fIv f. ." 0.2015"
1 s>f fdup cr f. space fIv f. ." 0.5963"
10 s>f fdup cr f. space fIv f. ." 0.9156"
100 s>f fdup cr f. space fIv f. ." ------"
[then]
\ Test the exponential Integral, E1(x).
true [if] \ refs: [1] & mathworld.wolfram.com
cr cr ." x fE1(x) Compared to:"
s" 1e-7" s>float fdup cr f. space fE1 f. ." 15.540880086"
s" 1e-6" s>float fdup cr f. space fE1 f. ." 13.238295893"
s" 1e-5" s>float fdup cr f. space fE1 f. ." 10.935719800"
s" 1e-4" s>float fdup cr f. space fE1 f. ." 8.633224705"
s" 1e-3" s>float fdup cr f. space fE1 f. ." 6.331539364"
s" 1e-2" s>float fdup cr f. space fE1 f. ." 4.037929577"
s" 1e-1" s>float fdup cr f. space fE1 f. ." 1.822923958"
s" 65e-2" s>float fdup cr f. space fE1 f. ." 0.411516976" \ ref [1]
s" 66e-2" s>float fdup cr f. space fE1 f. ." 0.403586275" \ ref [1]
s" 67e-2" s>float fdup cr f. space fE1 f. ." 0.395852563" \ ref [1]
s" 1e0" s>float fdup cr f. space fE1 f. ." 0.219383934" \ ref [1]
s" 1e1" s>float fdup cr f. space fE1 f. ." 4.156968930E-6"
cr
[then]
\ ---- cut here -----
\ EXAMPLE (Ivantsov function): Snowflakes and Dendrites.
\ --
\ From [4], Kenneth Libbrecht, "The enigmatic snowflake"
\ PysicsWorld [21] (2008) p. 19 (physicsworld.com):
\ "In 1947 the Russian mathematician GP Invantov discovered a family of
\ dynamically stable solutions to the diffusion equations that shed
\ considerable light on the growth of dendritic structures." ....
\ "Interestingly, ice forms nearly the same dendritic structures whether it
\ is grown from water vapour in air or from freezing liquid water. In the
\ latter case, the growth is mainly limited by the diffusion of latent heat
\ generated at the solid-liquid interface. In a snow crystal, on the other
\ hand, growth is mainly limited by the diffusion of water-vapour molecules
\ through the surrounding air"
\ --
\ From [5], J.S. Langer, "Dendrites, Viscous Fingers, and the Theory
\ of Pattern Formation" SCIENCE [243] (1989) p. 1150:
\ A description of the Invantsov solution is given as (and the paradox) --
\ Delta (dimensionless undercooling) = (Tmelting-T)/(L/C)
\ and Peclet_No = Vel*Radius/[2*thermalDiffusivity] with
\ the Ivantsov solution as Delta = Iv(Peclet_No)
\ where:
\ L = Latent heat of freezing
\ C = Heat capacity of liquid
\ a = Thermal diffusivity
\ Tm = equilibrium melting temperature and T=temperature of water
\ Vel = ice crystal growth velocity
\ Radius = ice crystal tip radius
\ --
\ Experimental Results for H2O from [6] Y. Teraoka, A. Saito, S. Okawa,
\ "Ice crystal growth in supercooled solution"
\ Int. J refrigeration [25] (2002) p. 218-225. -- we find:
\ T=-3 deg C, Radius=1E-5 m, Vel=0.0002 m/s
\ --
\ Looking up data for water we also find near freezing:
\ L = 334 J/g
\ C = 4 J/[g*K]
\ a = 135E-9 [m^2]/s
\ --
\ Thus Delta = 3/(334/4) = 0.036
\ and Peclet_No=(0.0002)(1E-5)/(2*135e-9) = 0.0074
\ and from the Ivantsov solution: Delta=Iv(0.0074) which is:
true [IF]
cr cr ." Growth of Ice Crystals:"
cr ." Peclet Number = 0.0074"
cr ." Iv(0.0074) = " s" 74e-4" s>float fIv f. ." dimensionless undercooling"
cr ." as compared to 0.036 dimensionless undercooling experimentally"
cr ." with ice crystal growth vel=0.0002 m/s,"
cr ." ice crystal tip Radius=1x10-5 m, and"
cr ." water temperature of -3 deg. Celsius where"
cr ." dimensionless undercooling=(Tmelting-T)/(latent_heat/heat_capacity)"
cr ." and Peclet Number= Vel*Radius/(2*thermal_Diffusivity)"
cr cr
[then]
\ ---cut here---
|
