Skip to content

Added complete elliptic functions for 1st and 2nd kind #135

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 10 commits into from
Dec 19, 2018
Merged

Added complete elliptic functions for 1st and 2nd kind #135

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
1e96e2d
Added Complete Elliptic Functions of 1st and 2nd kind
Sumeghtech Dec 16, 2018
3d1af22
Updated docs with elliptic1 and elliptic2
Sumegh-git Dec 16, 2018
c770b93
Added elliptic1 and elliptic2.
Sumegh-git Dec 17, 2018
0a26c34
Updated ellip.jl
Sumegh-git Dec 17, 2018
053da6b
Added wiki and dlmf link for ellip.jl
Sumegh-git Dec 17, 2018
9cefeb4
Updated src/SpecialFunctions.jl
Sumegh-git Dec 17, 2018
d8659f9
Renamed , Added tests from Wolframalpha
Sumeghtech Dec 17, 2018
88e0e77
Updated runtests.jl
Sumegh-git Dec 17, 2018
930a084
Updated runtests.jl (1)
Sumegh-git Dec 18, 2018
a642382
Final tests.
Sumegh-git Dec 18, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions docs/src/index.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,8 @@ libraries.
| [`invdigamma(x)`](@ref SpecialFunctions.invdigamma) | [invdigamma function](http://bariskurt.com/calculating-the-inverse-of-digamma-function/) (i.e. inverse of `digamma` function at `x` using fixed-point iteration algorithm) |
| [`trigamma(x)`](@ref SpecialFunctions.trigamma) | [trigamma function](https://en.wikipedia.org/wiki/Trigamma_function) (i.e the logarithmic second derivative of `gamma` at `x`) |
| [`polygamma(m,x)`](@ref SpecialFunctions.polygamma) | [polygamma function](https://en.wikipedia.org/wiki/Polygamma_function) (i.e the (m+1)-th derivative of the `lgamma` function at `x`) |
| [`ellipk(x)`](@ref SpecialFunctions.ellipk) | [complete elliptic integral of 1st kind](https://en.wikipedia.org/wiki/Elliptic_integral) (i.e evaluates complete elliptic integral of 1st kind at `x`) |
| [`ellipe(x)`](@ref SpecialFunctions.ellipe) | [complete elliptic integral of 2nd kind](https://en.wikipedia.org/wiki/Elliptic_integral) (i.e evaluates complete elliptic integral of 2nd kind at `x`) |
| [`eta(x)`](@ref SpecialFunctions.eta) | [Dirichlet eta function](https://en.wikipedia.org/wiki/Dirichlet_eta_function) at `x` |
| [`zeta(x)`](@ref SpecialFunctions.zeta) | [Riemann zeta function](https://en.wikipedia.org/wiki/Riemann_zeta_function) at `x` |
| [`airyai(z)`](@ref SpecialFunctions.airyai) | [Airy Ai function](https://en.wikipedia.org/wiki/Airy_function) at `z` |
Expand Down
2 changes: 2 additions & 0 deletions docs/src/special.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,8 @@ SpecialFunctions.besseli
SpecialFunctions.besselix
SpecialFunctions.besselk
SpecialFunctions.besselkx
SpecialFunctions.ellipk
SpecialFunctions.ellipe
SpecialFunctions.eta
SpecialFunctions.zeta
SpecialFunctions.gamma
Expand Down
5 changes: 4 additions & 1 deletion src/SpecialFunctions.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,8 @@ export
bessely1,
besselyx,
dawson,
ellipk,
ellipe,
erf,
erfc,
erfcinv,
Expand All @@ -57,12 +59,13 @@ export

include("bessel.jl")
include("erf.jl")
include("ellip.jl")
include("sincosint.jl")
include("gamma.jl")
include("deprecated.jl")

for f in (:digamma, :erf, :erfc, :erfcinv, :erfcx, :erfi, :erfinv,
:eta, :gamma, :invdigamma, :lfactorial, :lgamma, :trigamma)
:eta, :gamma, :invdigamma, :lfactorial, :lgamma, :trigamma, :ellipk, :ellipe)
@eval $(f)(::Missing) = missing
end
for f in (:beta, :lbeta)
Expand Down
202 changes: 202 additions & 0 deletions src/ellip.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,202 @@
using Base.Math: @horner

#Using piecewise approximation polynomial as given in
#'Fast Computation of Complete Elliptic Integrals and Jacobian Elliptic Functions'
# Fukushima, Toshio. (2014). F09-FastEI. Celest Mech Dyn Astr
# DOI 10.1007/s10569-009-9228-z
#Link : https://pdfs.semanticscholar.org/8112/c1f56e833476b61fc54d41e194c962fbe647.pdf
#
#For m<0 , followed Fukushima, Toshio. (2014). Precise, compact, and fast computation of complete elliptic integrals by piecewise minimax rational function approximation. Journal of Computational and Applied Mathematics. 282. 10.13140/2.1.1946.6245.
# Link: https://www.researchgate.net/profile/Toshio_Fukushima/publication/267330394_Precise_compact_and_fast_computation_of_complete_elliptic_integrals_by_piecewise_minimax_rational_function_approximation/links/544b81a40cf2d6347f43074f/Precise-compact-and-fast-computation-of-complete-elliptic-integrals-by-piecewise-minimax-rational-function-approximation.pdf?origin=publication_detail
#Also suggested in this paper that we should consider domain only from (-inf,1].



"""
ellipk(x)

DLMF : https://dlmf.nist.gov/19.2#E4 , https://dlmf.nist.gov/19.2#E8
Wiki : https://en.wikipedia.org/wiki/Elliptic_integral
Computes Complete Elliptic Integral of 1st kind at `x`-> K(x)--- given by: ``\\int_{0}^{\\pi/2} \\frac{1}{\\sqrt{1 - x(\\sin \\theta )^{2}}} d\\theta``
"""
function ellipk(a::Float64)
flag = false
if a<0.0
x=a/(a-1) #dealing with negative args
flag=true
elseif a>=0.0
x=a
flag=false
end

if x==0.0
return pi/2

elseif x==1.0
return Inf

elseif x>1.0
throw(DomainError(x,"`x` must lie between -Inf and 1 ---- Domain: (-Inf,1.0]"))

elseif 0.0 <= x < 0.1 #Table 2 from paper
t= x-0.05
t= @horner(t, 1.591003453790792180 , 0.416000743991786912 , 0.245791514264103415 , 0.179481482914906162 , 0.144556057087555150 , 0.123200993312427711 , 0.108938811574293531 , 0.098853409871592910 , 0.091439629201749751 , 0.085842591595413900 , 0.081541118718303215)

elseif 0.1 <= x <0.2 #Table 3
t=x-0.15
t= @horner(t , 1.635256732264579992 , 0.471190626148732291 , 0.309728410831499587 , 0.252208311773135699 , 0.226725623219684650 , 0.215774446729585976 , 0.213108771877348910 , 0.216029124605188282 , 0.223255831633057896 , 0.234180501294209925 , 0.248557682972264071 , 0.266363809892617521)

elseif 0.2 <= x < 0.3 #Table 4
t=x-0.25
t= @horner(t , 1.685750354812596043 , 0.541731848613280329 , 0.401524438390690257 , 0.369642473420889090 , 0.376060715354583645 , 0.405235887085125919 , 0.453294381753999079 , 0.520518947651184205 , 0.609426039204995055 , 0.724263522282908870 , 0.871013847709812357 , 1.057652872753547036)

elseif 0.3 <= x < 0.4 #Table 5
t= x-0.35
t= @horner(t , 1.744350597225613243 , 0.634864275371935304 , 0.539842564164445538 , 0.571892705193787391 , 0.670295136265406100 , 0.832586590010977199 , 1.073857448247933265 , 1.422091460675497751 , 1.920387183402304829 , 2.632552548331654201 , 3.652109747319039160 , 5.115867135558865806 , 7.224080007363877411)

elseif 0.4 <= x < 0.5 #Table 6
t=x-0.45
t= @horner(t, 1.813883936816982644 , 0.763163245700557246 , 0.761928605321595831 , 0.951074653668427927 , 1.315180671703161215 , 1.928560693477410941 , 2.937509342531378755 , 4.594894405442878062 , 7.330071221881720772 , 11.87151259742530180 , 19.45851374822937738 , 32.20638657246426863 , 53.73749198700554656 , 90.27388602940998849)

elseif 0.5 <= x < 0.6 #Table 7
t= x-0.55
t= @horner(t , 1.898924910271553526 , 0.950521794618244435 , 1.151077589959015808 , 1.750239106986300540 , 2.952676812636875180 , 5.285800396121450889 , 9.832485716659979747 , 18.78714868327559562 , 36.61468615273698145 , 72.45292395127771801 , 145.1079577347069102 , 293.4786396308497026 , 598.3851815055010179 , 1228.420013075863451 , 2536.529755382764488)
elseif 0.6 <= x < 0.7 #Table 8
t=x-0.65
t =@horner(t , 2.007598398424376302 , 1.248457231212347337 , 1.926234657076479729 , 3.751289640087587680 , 8.119944554932045802 , 18.66572130873555361 , 44.60392484291437063 , 109.5092054309498377 , 274.2779548232413480,697.5598008606326163 , 1795.716014500247129 , 4668.381716790389910 , 12235.76246813664335 , 32290.17809718320818 , 85713.07608195964685 , 228672.1890493117096 , 612757.2711915852774)

elseif 0.7 <= x <0.8 #Table 9
t=x-0.75
t= @horner(t, 2.156515647499643235 , 1.791805641849463243 , 3.826751287465713147 , 10.38672468363797208 , 31.40331405468070290 , 100.9237039498695416 , 337.3268282632272897 , 1158.707930567827917 , 4060.990742193632092 , 14454.00184034344795 , 52076.66107599404803 , 189493.6591462156887 , 695184.5762413896145 , 2567994.048255284686 , 9541921.966748386322 , 35634927.44218076174 , 133669298.4612040871 , 503352186.6866284541 , 1901975729.538660119 , 7208915015.330103756)

elseif 0.8 <= x <0.85 #Table 10
t = x-0.825
t =@horner(t , 2.318122621712510589 , 2.616920150291232841 , 7.897935075731355823 , 30.50239715446672327 , 131.4869365523528456 , 602.9847637356491617 , 2877.024617809972641 , 14110.51991915180325 , 70621.44088156540229 , 358977.2665825309926 , 1847238.263723971684 , 9600515.416049214109 , 50307677.08502366879 , 265444188.6527127967 , 1408862325.028702687 , 7515687935.373774627)

elseif 0.85 <= x <0.9 #Table 11
t= x-0.875
t =@horner(t, 2.473596173751343912 , 3.727624244118099310 , 15.60739303554930496 , 84.12850842805887747 , 506.9818197040613935 , 3252.277058145123644 , 21713.24241957434256 , 149037.0451890932766 , 1043999.331089990839 , 7427974.817042038995 , 53503839.67558661151 , 389249886.9948708474 , 2855288351.100810619 , 21090077038.76684053 , 156699833947.7902014 , 1170222242422.439893 , 8777948323668.937971 , 66101242752484.95041 , 499488053713388.7989 , 37859743397240299.20)

elseif x>=0.9
td= 1-x
td1=td-0.05
qd = @horner(td, 0.0, (1.0/16.0), (1.0/32.0), (21.0/1024.0), (31.0/2048.0), (6257.0/524288.0), (10293.0/1048576.0), (279025.0/33554432.0), (483127.0/67108864.0), (435506703.0/68719476736.0), (776957575.0/137438953472.0) , (22417045555.0/4398046511104.0) , (40784671953.0/8796093022208.0) , (9569130097211.0/2251799813685248.0) , (17652604545791.0/4503599627370496.0))

kdm = @horner(td1 , 1.591003453790792180 , 0.416000743991786912 , 0.245791514264103415 , 0.179481482914906162 , 0.144556057087555150 , 0.123200993312427711 , 0.108938811574293531 , 0.098853409871592910 , 0.091439629201749751 , 0.085842591595413900 , 0.081541118718303215)
km = -Base.log(qd) * (kdm/pi)
t = km
end
if flag == false
return t
elseif flag == true
ans = t / sqrt(1.0-a)
return ans
end
end



#Using piecewise approximation polynomial as given in
#'Fast Computation of Complete Elliptic Integrals and Jacobian Elliptic Functions'
# Fukushima, Toshio. (2014). F09-FastEI. Celest Mech Dyn Astr
# DOI 10.1007/s10569-009-9228-z
#Link : https://pdfs.semanticscholar.org/8112/c1f56e833476b61fc54d41e194c962fbe647.pdf
#
#For m<0 , followed Fukushima, Toshio. (2014). Precise, compact, and fast computation of complete elliptic integrals by piecewise minimax rational function approximation. Journal of Computational and Applied Mathematics. 282. 10.13140/2.1.1946.6245.
# Link: https://www.researchgate.net/profile/Toshio_Fukushima/publication/267330394_Precise_compact_and_fast_computation_of_complete_elliptic_integrals_by_piecewise_minimax_rational_function_approximation/links/544b81a40cf2d6347f43074f/Precise-compact-and-fast-computation-of-complete-elliptic-integrals-by-piecewise-minimax-rational-function-approximation.pdf?origin=publication_detail
#Also suggested in this paper that we should consider domain only from (-inf,1].


"""
ellipe(x)

DLMF : https://dlmf.nist.gov/19.2#E5 , https://dlmf.nist.gov/19.2#E8
Wiki : https://en.wikipedia.org/wiki/Elliptic_integral
Computes the Complete Elliptic Integral of 2nd kind at `x` ->E(x)---gven by: ``\\int_{0}^{\\pi/2} \\sqrt{1 - x(\\sin \\theta )^{2}} d\\theta``
"""
function ellipe(a::Float64)
flag=false
if a<0.0
x=a/(a-1) #for dealing with negative args
flag=true
elseif a>=0.0
x=a
flag=false

end

if x==0.0
return pi/2
elseif x==1.0
return 1.0

elseif x>1.0
throw(DomainError(x,"`x` must lie between -Inf and 1 ---- Domain : (-inf,1.0]"))

elseif 0.0 <= x < 0.1 #Table 2 from paper
t= x-0.05
t= @horner(t , +1.550973351780472328 ,-0.400301020103198524 ,-0.078498619442941939 ,-0.034318853117591992 ,-0.019718043317365499 ,-0.013059507731993309 ,-0.009442372874146547 ,-0.007246728512402157 ,-0.005807424012956090 ,-0.004809187786009338)

elseif 0.1 <= x < 0.2 #Table 3
t=x-0.15
t= @horner(t , +1.510121832092819728 , -0.417116333905867549 , -0.090123820404774569 , -0.043729944019084312, -0.027965493064761785 , -0.020644781177568105 , -0.016650786739707238 , -0.014261960828842520 , -0.012759847429264803 , -0.011799303775587354 , -0.011197445703074968)

elseif 0.2 <= x < 0.3 #Table 4
t=x-0.25
t= @horner(t , 1.467462209339427155 , -0.436576290946337775 , -0.105155557666942554 , -0.057371843593241730 , -0.041391627727340220 , -0.034527728505280841 , -0.031495443512532783 , -0.030527000890325277 , -0.030916984019238900 , -0.032371395314758122 , -0.034789960386404158)

elseif 0.3 <= x < 0.4 #Table 5
t= x-0.35
t= @horner(t, +1.422691133490879171, -0.459513519621048674, -0.125250539822061878, -0.078138545094409477, -0.064714278472050002, -0.062084339131730311, -0.065197032815572477,-0.072793895362578779, -0.084959075171781003, -0.102539850131045997, -0.127053585157696036, -0.160791120691274606)

elseif 0.4 <= x < 0.5 #Table 6
t=x-0.45
t= @horner(t ,+1.375401971871116291 ,-0.487202183273184837 ,-0.153311701348540228 ,-0.111849444917027833 ,-0.108840952523135768 ,-0.122954223120269076 ,-0.152217163962035047 ,-0.200495323642697339 ,-0.276174333067751758 ,-0.393513114304375851 ,-0.575754406027879147 ,-0.860523235727239756 ,-1.308833205758540162)

elseif 0.5 <= x < 0.6 #Table 7
t= x-0.55
t= @horner(t ,+1.325024497958230082 ,-0.521727647557566767 ,-0.194906430482126213 ,-0.171623726822011264 ,-0.202754652926419141 ,-0.278798953118534762 ,-0.420698457281005762 ,-0.675948400853106021 ,-1.136343121839229244 ,-1.976721143954398261 ,-3.531696773095722506 ,-6.446753640156048150 ,-11.97703130208884026)
elseif 0.6 <= x < 0.7 #Table 8
t=x-0.65
t= @horner(t, +1.270707479650149744, -0.566839168287866583, -0.262160793432492598, -0.292244173533077419, -0.440397840850423189, -0.774947641381397458, -1.498870837987561088, -3.089708310445186667, -6.667595903381001064, -14.89436036517319078, -34.18120574251449024, -80.15895841905397306, -191.3489480762984920, -463.5938853480342030, -1137.380822169360061)

elseif 0.7 <= x < 0.8 #Table 9
t=x-0.75
t= @horner(t, +1.211056027568459525, -0.630306413287455807, -0.387166409520669145, -0.592278235311934603, -1.237555584513049844, -3.032056661745247199, -8.181688221573590762, -23.55507217389693250, -71.04099935893064956, -221.8796853192349888, -712.1364793277635425, -2336.125331440396407, -7801.945954775964673, -26448.19586059191933, -90799.48341621365251, -315126.0406449163424, -1104011.344311591159)

elseif 0.8 <= x < 0.85 #Table 10
t = x-0.825
t= @horner(t, +1.161307152196282836, -0.701100284555289548, -0.580551474465437362, -1.243693061077786614, -3.679383613496634879, -12.81590924337895775, -49.25672530759985272, -202.1818735434090269, -869.8602699308701437, -3877.005847313289571, -17761.70710170939814, -83182.69029154232061, -396650.4505013548170, -1920033.413682634405)

elseif 0.85 <= x < 0.9 #Table 11
t= x-0.875
t = @horner(t, +1.124617325119752213, -0.770845056360909542, -0.844794053644911362, -2.490097309450394453, -10.23971741154384360, -49.74900546551479866, -267.0986675195705196, -1532.665883825229947, -9222.313478526091951, -57502.51612140314030, -368596.1167416106063, -2415611.088701091428, -16120097.81581656797, -109209938.5203089915, -749380758.1942496220, -5198725846.725541393, -36409256888.12139973)

elseif x>=0.9
td= 1-x
td1=td-0.05

kdm = @horner(td1 , 1.591003453790792180 , 0.416000743991786912 , 0.245791514264103415 , 0.179481482914906162 , 0.144556057087555150 , 0.123200993312427711,0.108938811574293531 , 0.098853409871592910 ,0.091439629201749751 ,0.085842591595413900 ,0.081541118718303215)
edm = @horner(td1 ,+1.550973351780472328 ,-0.400301020103198524 ,-0.078498619442941939 ,-0.034318853117591992,-0.019718043317365499 ,-0.013059507731993309 ,-0.009442372874146547 ,-0.007246728512402157 ,-0.005807424012956090 ,-0.004809187786009338)
hdm = kdm -edm
km = ellipk(Float64(x))
#em = km + (pi/2 - km*edm)/kdm
em = (pi/2 + km*hdm)/kdm #to avoid precision loss near 1
t= em
end
if flag == false
return t
elseif flag == true
return t * sqrt(1.0-a)

end
end

for f in (:ellipk,:ellipe)
@eval begin
($f)(x::Float16) = Float16(($f)(Float64(x)))
($f)(x::Float32) = Float32(($f)(Float64(x)))
($f)(x::Real) = ($f)(float(x))
($f)(x::AbstractFloat) = throw(MethodError($f,(x,"")))
end
end
34 changes: 33 additions & 1 deletion test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,39 @@ relerrc(z, x) = max(relerr(real(z),real(x)), relerr(imag(z),imag(x)))
@test_throws MethodError erf(big(1.0)*im)
@test_throws MethodError erfi(big(1.0))
end

@testset "elliptic integrals" begin
#Computed using Wolframalpha EllipticK and EllipticE functions.
@test ellipk(0) ≈ 1.570796326794896619231322 rtol=2*eps()
@test ellipk(0.92) ≈ 2.683551406315229344 rtol=2*eps()
@test ellipk(0.5) ≈ 1.854074677301371918 rtol=2*eps()
@test ellipk(0.01) ≈ 1.57474556151735595 rtol=2*eps()
@test ellipk(0.45) ≈ 1.81388393681698264 rtol=2*eps()
@test ellipk(-0.5) ≈ 1.41573720842595619 rtol=2*eps()
@test ellipk(0.75) ≈ 2.15651564749964323 rtol=2*eps()
@test ellipk(0.17) ≈ 1.6448064907988806 rtol=2*eps()
@test ellipk(0.25) ≈ 1.685750354812596 rtol=2*eps()
@test ellipk(0.69) ≈ 2.0608816467301313 rtol=2*eps()
@test ellipk(0.84) ≈ 2.3592635547450067 rtol=2*eps()
@test ellipe(0.15) ≈ 1.5101218320928197 rtol=2*eps()
@test ellipe(0.21) ≈ 1.4847605813318776 rtol=2*eps()
@test ellipe(0.42) ≈ 1.3898829914929717 rtol=2*eps()
@test ellipe(0.66) ≈ 1.2650125751607508 rtol=2*eps()
@test ellipe(0.76) ≈ 1.2047136418292115 rtol=2*eps()
@test ellipe(0.865) ≈ 1.1322436887003925 rtol=2*eps()
@test ellipe(0) ≈ 1.570796326794896619231322 rtol=2*eps()
@test ellipe(0.8) ≈ 1.17848992432783852 rtol=2*eps()
@test ellipe(0.5) ≈ 1.3506438810476755 rtol=2*eps()
@test ellipe(0.01) ≈ 1.5668619420216682 rtol=2*eps()
@test ellipe(0.99) ≈ 1.0159935450252239 rtol=2*eps()
@test ellipe(-0.1) ≈ 1.6093590249375295 rtol=2*eps()
@test ellipe(0.3) ≈ 1.4453630644126652 rtol=2*eps()
@test ellipe(1.0) ≈ 1.00
@test ellipk(1.0)==Inf
@test_throws MethodError ellipk(BigFloat(0.5))
@test_throws MethodError ellipe(BigFloat(-1))
@test_throws DomainError ellipe(Float16(2.0))
@test_throws DomainError ellipe(Float32(2.5))
end
@testset "sine and cosine integrals" begin
# Computed via wolframalpha.com: SinIntegral[SetPrecision[Table[x,{x, 1,20,1}],20]] and CosIntegral[SetPrecision[Table[x,{x, 1,20,1}],20]]
sinintvals = [0.9460830703671830149, 1.605412976802694849, 1.848652527999468256, 1.75820313894905306, 1.54993124494467414, 1.4246875512805065, 1.4545966142480936, 1.5741868217069421, 1.665040075829602, 1.658347594218874, 1.578306806945727416, 1.504971241526373371, 1.499361722862824564, 1.556211050077665054, 1.618194443708368739, 1.631302268270032886, 1.590136415870701122, 1.536608096861185462, 1.518630031769363932, 1.548241701043439840]
Expand Down