From 9e3a88c2f1458ce7cbf44ce416a27a0fdc64debe Mon Sep 17 00:00:00 2001 From: elaineVRC Date: Sun, 26 Feb 2017 23:49:43 -0800 Subject: [PATCH 1/7] recNM3.ipynb Matlab doesn't calculate associated Legendre type 3 (z>1) function but does do type 2 (-1<=x<=1) approx. .0007 sec compared to average .00014sec recNM3(n,m,z) which is a recursive implementation. --- recNM3.ipynb | 58 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) create mode 100644 recNM3.ipynb diff --git a/recNM3.ipynb b/recNM3.ipynb new file mode 100644 index 00000000..7114b8ea --- /dev/null +++ b/recNM3.ipynb @@ -0,0 +1,58 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "function recNM3(n,m,z) #n,m positive integers, z number \n", + " #reliable for values of n <= 17 for z=2. \n", + " # 0 <= m <= n\n", + " # z > 1.\n", + " M = (2*m -1) # M must be odd\n", + " # (2m-1)!!= (2m)!/( m! 2^m) double factorial \n", + " dblfac=1\n", + " for j=1:M\n", + " if iseven(j)\n", + " continue\n", + " end\n", + " dblfac=j*dblfac \n", + " end\n", + " if n == m\n", + " return( dblfac*(z*z -1)^(m/2)) \n", + " elseif n == m+1\n", + " return((2.*m +1.)*z*dblfac*(z*z -1)^(m/2))\n", + " end\n", + " pj2=dblfac*(z*z -1.)^(m/2)\n", + " pj1=z*(2.*m+1.)*pj2 \n", + " for j = m+2 :n \n", + " pjj=(z*(2.*j-1.)*pj1 - pj2*(j +m-1.)) /(j-m)\n", + " pj2=pj1\n", + " pj1=pjj\n", + " end\n", + "return pj1 \n", + "end\n", + "\n", + "\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 0.4.2", + "language": "julia", + "name": "julia-0.4" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "0.4.2" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} From 390ab128833a6a36d83df3f0fa7be63d1ae76c51 Mon Sep 17 00:00:00 2001 From: elaineVRC Date: Tue, 7 Mar 2017 22:41:50 -0800 Subject: [PATCH 2/7] associatedLegendre.jl Julia file that calculates some associated Legendre functions --- src/associatedLegendre.jl | 105 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 105 insertions(+) create mode 100644 src/associatedLegendre.jl diff --git a/src/associatedLegendre.jl b/src/associatedLegendre.jl new file mode 100644 index 00000000..93680502 --- /dev/null +++ b/src/associatedLegendre.jl @@ -0,0 +1,105 @@ + +# solutions, w, of the differential equation are +# associated Legendre functions + +# (1-x*x)d/dx(dw/dx)⁢-2⁢x⁢(dw/dx)+(n(n+1)-m*m/(1-x*x))w=0 + + #for m =0 they are a.k.a. Legendre polynomials +#n is the degree and m is the order +#The recursion equation 8.5.3 Abramowitz & Stegun (AS) Handbook of +# Mathematical Functions used is +# (n-m+1)P(n+1,m,z)=(2n+1)zP(n,m,z)- (n+m)P(n-1,m,z) +# (for a recent comparison of transforms and recursions +# see https:doi.org/10.1515/jag-2016-0032) +#this recursion is valid for legendre functions of the first kind +#(P(n,m,z))and second kind Q(n,m,z) , +# (for types 1,2,3 ,definition from Mathematica documentation +# type 1 (?) is m=0, type 2 is -1<=x<=1 and type 3 is z> 1 +# http://mpmath.org/doc/0.18/functions/orthogonal.html +#Thus we define functions LegendreP2,....P3,....Q2,...Q3 +# for n and m positive integers 0<=m<=n +# refer to http://dlmf.nist.gov/14.6 +#P(n,m,z) and Q(n,m,⁡z) (z> 1) are often referred to as the prolate +#spheroidal harmonics of the first and second kinds, respectively +# We can define spherical harmonics (there are other definitions) +# Y(n,m,theta,phi)= LegendreP2(n,|m|,theta)*exp(im*m*phi)*((-)^(|m|-m)/2)) +# *sqrt((2n+1)*factorial(n-|m|)/(4*pi*factorial(n+|m|))) +# refer to https://en.wikipedia.org/wiki/Spherical_harmonics + + + +function LegendreP2(n::Integer,m::Integer,x::Number) + + + # 0 <= m <= n + # -1<=x <= 1 + M = (2*m -1) # M must be odd + # (2m-1)!!= (2m)!/( m! 2^m) double factorial + dblfac=1 + for j=1:M + if iseven(j) + continue + end + dblfac=j*dblfac + end + + pj2= ((-1)^m)*dblfac*(1-x*x)^(m/2) + pj1=x*(2.*m+1.)*pj2 + if n == m + return pj2 + elseif n == m+1 + return pj1 + end + + + for j = m+2 :n + pjj=(x*(2.*j-1.)*pj1 - pj2*(j +m-1.)) /(j-m) + pj2=pj1 + pj1=pjj + end +return pj1 +end + + + + + + +function LegendreP3(n::Integer,m::Integer,z::Number) + + #reliable for values of n <= 17 for z=2. + # 0 <= m <= n + # z > 1. + M = (2*m -1) # M must be odd + # (2m-1)!!= (2m)!/( m! 2^m) double factorial + dblfac=1 + for j=1:M + if iseven(j) + continue + end + dblfac=j*dblfac + end + + + pj2=dblfac*(z*z -1.)^(m/2) + pj1=z*(2.*m+1.)*pj2 + + if n==m + return (pj2) + elseif n== m+1 + return (pj1) + end + + + for j = m+2 :n + pjj=(z*(2.*j-1.)*pj1 - pj2*(j +m-1.)) /(j-m) + pj2=pj1 + pj1=pjj + end +return pj1 +end + + + + + From 5c805031ea55159bba3cc9826b90ae75c30a1794 Mon Sep 17 00:00:00 2001 From: elaineVRC Date: Fri, 31 Mar 2017 17:39:26 -0700 Subject: [PATCH 3/7] Update associatedLegendre.jl On line 46 changed 1 to one(x) On line 84 changed 1. to one(z) for arbitrary precision, etc. --- src/associatedLegendre.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/associatedLegendre.jl b/src/associatedLegendre.jl index 93680502..7b0e32e3 100644 --- a/src/associatedLegendre.jl +++ b/src/associatedLegendre.jl @@ -43,7 +43,7 @@ function LegendreP2(n::Integer,m::Integer,x::Number) dblfac=j*dblfac end - pj2= ((-1)^m)*dblfac*(1-x*x)^(m/2) + pj2= ((-1)^m)*dblfac*(one(x)-x*x)^(m/2) pj1=x*(2.*m+1.)*pj2 if n == m return pj2 @@ -81,7 +81,7 @@ function LegendreP3(n::Integer,m::Integer,z::Number) end - pj2=dblfac*(z*z -1.)^(m/2) + pj2=dblfac*(z*z - one(z))^(m/2) pj1=z*(2.*m+1.)*pj2 if n==m From 82a1469d246f611db0f3d611adb4d456c3f7b621 Mon Sep 17 00:00:00 2001 From: elaineVRC Date: Mon, 22 May 2017 17:37:30 -0700 Subject: [PATCH 4/7] associatedLegendre.jl Thanks for editting and making the PR more readable; congratulations on your new position AA. Thanks to SJ for mentoring and encouragement. I can't benchmark against what does not exist: matlab, maxima, scipy.special, gsl do not have associated Legendre type 3 (z > 1). I have used SLegendreP3 (which I derived and used a long time ago) ; similarly SLegendreP2 gives reliable results compared to scipy.special for |x| <1 . function BenchmarkZZLegendreP2 , shows that ZZLegendreP2 is faster than scipy.special. ZZLegendreP3 comes from a modification of work published by Selezneva, etal. If the general direction of the PR is OK; then I can extract all the test... functions to test.jl and do additional work. All suggestions, comments,etc. are welcome. Thank you for helping me to learn Julia scientific computing. --- associatedLegendre.jl | 596 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 596 insertions(+) create mode 100644 associatedLegendre.jl diff --git a/associatedLegendre.jl b/associatedLegendre.jl new file mode 100644 index 00000000..d015035f --- /dev/null +++ b/associatedLegendre.jl @@ -0,0 +1,596 @@ + +# solutions, w, of the differential equation are +# associated Legendre functions +# (1-x*x)d/dx(dw/dx)⁢-2⁢x⁢(dw/dx)+(n(n+1)-m*m/(1-x*x))w=0 +#for m =0 they are a.k.a. Legendre polynomials +#n is the degree and m is the order +# +#There are more than a handful of programs. Some fast and inaccurate +#and some slow and more accurate; pick your poison (trade off) . +# +# SSLegendreP2 and 3 : The finite sum is arranged so that +# each term is a multiple of the preceeding term. SLegendreP2 and 3 +# are very slow and used for testing. +# +# The recursion equation 8.5.3 Abramowitz & Stegun (AS) Handbook of +# Mathematical Functions used in LegendreP2 is +# (n-m+1)P(n+1,m,z)=(2n+1)zP(n,m,z)- (n+m)P(n-1,m,z) +# +# +# ZZLegendreP2 and 3 uses a different recursion +# P(n,m+2,x)+2(m+1)x(s(1-x*x))^(-1/2)P(n,m+1,x)+s(n-m)(n+m+1)P(n,m,x)=0 +#from http://dlmf.nist.gov/14.10 equation 14.10.1 where s=sign(1-x*x) +#refer to Lebedev, Special Functions, Dover Pub. eqn. 7.12.8 page 194 +# See function BenchmarkZZLegemdreP2() for speed. +# +#These recursions are valid for Legendre functions of the first kind +#(P(n,m,z))and second kind Q(n,m,z) , +#(for types 1,2,3 ,definition from Mathematica documentation +# type 1 (?) is m=0, type 2 is -1<=x<=1 and type 3 is z> 1 +# http://mpmath.org/doc/0.18/functions/orthogonal.html +#Thus we define functions LegendreP2,....P3,....Q2,...Q3 +# for n and m positive integers 0<=m<=n +# refer to http://dlmf.nist.gov/14.6 +#P(n,m,z) and Q(n,m,⁡z) (z> 1) are often referred to as the prolate +#spheroidal harmonics of the first and second kinds, respectively +# We can define spherical harmonics (there are other definitions) +# Y(n,m,theta,phi)= LegendreP2(n,|m|,theta)*exp(im*m*phi)*((-)^(|m|-m)/2)) +# *sqrt((2n+1)*factorial(n-|m|)/(4*pi*factorial(n+|m|))) +# refer to https://en.wikipedia.org/wiki/Spherical_harmonics + + + + + + +""" + LegendreP2(n::Integer,m::Integer,x::Number) + associated Legendre function of first kind type 2 (-1.<=x<=1.). + + +""" + +function LegendreP2(n::Integer,m::Integer,x::Number) + + + # 0 <= m <= n + # -1<=x <= 1 + M = (2*m -1) # M must be odd + # (2m-1)!!= (2m)!/( m! 2^m) double factorial + dblfac=1 + for j=1:M + if iseven(j) + continue + end + dblfac=j*dblfac + end + + pj2= ((-1)^m)*dblfac*(one(x)-x*x)^(m/2) + pj1=x*(2.*m+1.)*pj2 + if n == m + return pj2 + elseif n == m+1 + return pj1 + end + + + for j = m+2 :n + pjj=(x*(2.*j-1.)*pj1 - pj2*(j +m-1.)) /(j-m) + pj2=pj1 + pj1=pjj + end +return pj1 +end + + + + +#checking LegendreP2 against scipy.special +# for n< 18 the relative error is less than 2e-14 +# the error depends on the n amd m values and x values to some extent + +using PyCall + +@pyimport scipy.special as s + +""" + function testLegendreP2() + relative error LegendreP2 versus scipy special +""" + + +function testLegendreP2() + x=.6 + rerr = 2e-14 + for n=0:17 + for m=0:n #m=rand(0:n) + a=LegendreP2(n,m,x) + b=s.lpmv(m,n,x) + df=norm(a-b)/norm(b) + if ( df > rerr) + #if (norm(a-b)/norm(b)) > rerr) + println(n," ", m," ", a," ",b," ",df) + #@test_approx_eq_eps a b 1e-14 + end + end + + end + +end + +testLegendreP2() + +""" + SLegendreP2(n,m,x) + Legendre function 1st kind type 2 , -1.<=x<=1. + series (to be optimized for in SSLegendreP2) +""" + +function SLegendreP2(n,m,x)#n,m integers + + v=(big(1.)-x*x)^(-m*big(1)/2) + v=v*factorial(big(n+m))/(((big(2))^n)* factorial(big(n))); + term = big(0.); + #if iseven(n+m); MXP=Int(trunc((n+m)/2));end; + #if isodd(n+m);MXP=Int(trunc((n+m-1)/2));end; + #MXP=Int(floor((n+m)/2)); + + MXP=Int(trunc(((n+m)/2))) + + + for p=0:MXP + + term= term + + ((-big(1.))^p)*binomial(big(2*(n-p)),big(n-m))* + (x^(big(1)*(n + m -2*p)))* + binomial(big(n),big(p)) + + end + return term * v +end + +# I derived this finite sum using 8.6.6,8.6.18,8.2.5 (AS), +#binomial expansion and differentiating + + + + + +""" + function testLegendreP2F() + compare LegendreP2 against Float64(SLegendreP2) +""" +function testLegendreP2F() +#= + checking LegendreP2 against Float64SLegendreP2 + for n < 18 the relative error is less than 7e-15 + the error depends on the n amd m values and + x values to some extent + =# + + + x=.6 + rerr = 7e-15 + for n=1:17 + for m=0:n #rand(1:n) + a=LegendreP2(n,m,x) + b=Float64(SLegendreP2(n,m,big(x))) + df=norm(a-b)/norm(b) + if ( df > rerr) + println(n," ", m," ", a," ",b," ",df) + end + end +end +end +testLegendreP2F() + + + + +""" + SSLegendreP2(n::Integer,m::Integer,x::Number) +associated Legendre function 1st kind type 2 , -1.<=x<=1. +each term in sum is multiple of previous, from SLegendreP2 + + +""" + +function SSLegendreP2(n::Integer,m::Integer,x::Number)#n,m integers + + M = (2*n -1) # M must be odd + # (2m-1)!!= (2m)!/( m! 2^m) double factorial + dblfac=1 + for j=1:M + if iseven(j) + continue + end + dblfac=j*dblfac + end + + MXP=Int(trunc(((n+m)/2))) + sum= (x^(n+m))*dblfac/(((one(x)-x*x)^(m/2))*prod(1:n-m)) + prev=sum + for p=1:MXP + term=(-prev)*(n-p+1)*(n+m-2*p+1)*(n+m-2*p+2 + )/(p*(2*n-2*p+1)*(2*n-2*p+2)*x*x) + sum=sum+term + prev=term + end + return sum +end + + + +#checking SSLegendreP2 against scipy.special +# for n< 17 the relative error is less than 3e-11 +# the error depends on the n amd m values and x values to some extent + +using PyCall + +@pyimport scipy.special as s + +""" + function testSSLegendreP2() + compare SSLegendreP2 with scipy special +""" +function testSSLegendreP2() + x=.6 + rerr = 3e-11 + for n=1:17 + for m=0:n #m=rand(1:n) + a=SSLegendreP2(n,m,x) + b=s.lpmv(m,n,x) + df=norm(a-b)/norm(b) + if ( df > rerr) + #if (norm(a-b)/norm(b)) > rerr) + println(n," ", m," ", a," ",b," ",df) + #@test_approx_eq_eps a b 1e-14 + end + end + + end + +end +testSSLegendreP2() + +#checking SSLegendreP2 against Float64(SLegendreP2) +# for n< 17 the relative error is less than 3e-11 +# the error depends on the n amd m values and x values to some extent +""" + function testSSLegendreP2F() + compare SSLegendreP2 to float64SLegendreP2 +""" +function testSSLegendreP2F() + x=.6 + rerr = 3e-11 + for n=0:17 + for m=0:n #m=rand(0:n) + a=SSLegendreP2(n,m,x) + b=Float64(SLegendreP2(n,m,big(x))) + df=norm(a-b)/norm(b) + if ( df > rerr) + println(n," ", m," ", a," ",b," ",df) + end + end + end +end +testSSLegendreP2F() + + +""" +SLegendreP3(n::Integer,m::Integer,z::Number) + associated Legendre function z>1 1st kind type 3 + series expansion +""" +function SLegendreP3(n::Integer,m::Integer,z::Number) +#function mylegenp3(n,m,z)#n,m integers + v=(z*z - big(1.))^(-m*big(1)/2) + + v=v*factorial(big(n+m))/(((big(2))^n)* factorial(big(n))); + term = big(0.); + #if iseven(n+m); MXP=Int(trunc((n+m)/2));end; + #if isodd(n+m);MXP=Int(trunc((n+m-1)/2));end; + #MXP=Int(floor((n+m)/2)); + MXP=Int(trunc(((n+m)/2))) + + for p=0:MXP + + term= term + + ((-big(1.))^p)*binomial(big(2*(n-p)),big(n-m))* + (z^(big(1)*(n + m -2*p)))* + binomial(big(n),big(p)) + + end + return term * v +end + +#print mylegenp3(2,1,2.)#agrees to 16 decimals +#print mylegenp3(2,1,mpf(2))#agrees +#print mylegenp3(2,1,mpf('2.'))#agrees 10.392304845413263761 +#println(mylegenp3(2,1,big(2.))) AGREES with mpmath + +# I derived this finite sum using 8.6.6,8.6.18,8.2.5 (AS), +#binomial expansion and differentiating + + +function timeS2() + for i=1:10 + SLegendreP2(100,50,big(.6)) + end +end +TS2= @elapsed(timeS2()) +println("time SLegendreP2(100,50,big(.6)) $TS2") + +function timeP2() + for i=1:10 + SSLegendreP2(100,50,big(.6)) + end +end +TP2= @elapsed(timeP2()) +println("time SSLegendreP2(100,50,big(.6)) $TP2") +#time SLegendreP2(100,50,big(.6)) 0.0085234 +#time SSLegendreP2(100,50,big(.6)) 0.004890456 +# the series is double speed when each term is a multiple of previous term + + + +""" + SSLegendreP3(n::Integer,m::Integer,x::Number) + associated Legendre function 1st kind type 3 , x >1. + each term in sum is multiple of previous. + +""" + +function SSLegendreP3(n::Integer,m::Integer,x::Number)#n,m integers + + + M = (2*n -1) # M must be odd + # (2m-1)!!= (2m)!/( m! 2^m) double factorial + dblfac=1 + for j=1:M + if iseven(j) + continue + end + dblfac=j*dblfac + end + + MXP=Int(trunc(((n+m)/2))) + sum= (x^(n+m))*dblfac/(((x*x - one(x))^(m/2))*prod(1:n-m)) + prev=sum + for p=1:MXP + term=(-prev)*(n-p+1)*(n+m-2*p+1)*(n+m-2*p+2 + )/(p*(2*n-2*p+1)*(2*n-2*p+2)*x*x) + sum=sum+term + prev=term + end + return sum +end + + + + +#checking SSLegendreP3 against Float64(SLegendreP3) +# for n< 17 the relative error is less than 4e-14 +# the error depends on the n amd m values and x values to some extent +""" + function testSSLegendreP3F() + compare SSLegendreP3 with Float64 SLegendreP3 +""" +function testSSLegendreP3F() + x=2. + rerr = 4e-14 + for n=0:17 + for m=rand(0:n) + a=SSLegendreP3(n,m,x) + b=Float64(SLegendreP3(n,m,big(x))) + df=norm(a-b)/norm(b) + if ( df > rerr) + println(n," ", m," ", a," ",b," ",df) + end + end +end +end +testSSLegendreP3F() + +# based on code published at +# i. a. selezneva, yu. l. ratis, e. hernández, j. pérez-quiles +#and p. fernández de córdoba: +#a code to calculate high order legendre polynomials +#rev. acad. colomb. cienc.: volumen xxxvii, número 145 - diciembre 2013 +#www.scielo.org.co/pdf/racefn/v37n145/v37n145a09.pdf + #Selezneva.....pdf + +""" +ZZLegendreP2(n::Integer,m::Integer,z::Number) + associated Legendre function First kind type 2 |x|<=1. + uses different three term recursion +""" + +function ZZLegendreP2(n::Integer,m::Integer,z::Number) + nz=1 + pnm = zeros(typeof(z),2*n+1) + fac = prod(2.:n) + sqz2 = sqrt((one(z)-z.*z)) + hsqz2 = 0.5*sqz2 + ihsqz2 = z./hsqz2 + if(n==0) + pnm[1]=one(z) + return ((-one(z))^m)*pnm[n+1+m] #pnm + end + if(n==1) + pnm[1]=-.5*sqz2 + pnm[2]=z + pnm[3]=sqz2 + return ((-one(z))^m)*pnm[n+1+m] #pnm + end + pnm[1] = (1-2*abs(n-2*floor(n/2)))*hsqz2.^n/fac + pnm[2] = -pnm[1]*n.*ihsqz2 + for mr=1:2*n-1 + pnm[mr+2]=(mr-n).*ihsqz2.*pnm[mr+1]-(2*n-mr+1)*mr*pnm[mr] + end + return ((-one(z))^m)*pnm[n+1+m] +end + + + + +#Benchmarking ZZLegendreP2 against scipy.special + +using PyCall + +@pyimport scipy.special as s + + +function timePY() + for i=1:100 + s.lpmv(15,30,.6) + end +end + +function timeZ2() + for i=1:100 + ZZLegendreP2(30,15,.6) + end +end +""" + function BenchmarkZZLegendreP2() + timing for ZZLegendreP2 LegendreP2 + more than five times as fast + compared to s.lpmv from scipy.special + speed may depend on n,m,x values +""" + +function BenchmarkZZLegendreP2() + +TPY= @elapsed(timePY()) +println("time s.lpmv(15,30,.6) $TPY") + + +TZ2= @elapsed(timeZ2()) +println("time ZZLegendreP2(30,15,.6) $TZ2") + +#time s.lpmv(15,30,.6) 0.000674941 +#time ZZLegendreP2(30,15,.6) 0.000112901 + +end +BenchmarkZZLegendreP2() + +#checking ZZLegendreP2 against scipy special +# for n< 51 the relative error is less than 6e-7 +# the error depends on the n amd m values and x values to some extent +# +""" + function testZZLegendreP2() +compare ZZLegendreP2 with scipy special +""" +function testZZLegendreP2() + x=.6 + rerr = 6e-7 + for n=1:50 + for m= 0:n #rand(1:n) + a=ZZLegendreP2(n,m,x) + b=s.lpmv(m,n,x) + df=norm(a-b)/norm(b) + if ( df > rerr) + println(n," ", m," ", a," ",b," ",df) + end + end + end +end +testZZLegendreP2() + + +#checking ZZLegendreP2 against Float64(SLegendreP2) +# for n< 51 the relative error is less than 6e-7 +# the error depends on the n amd m values and x values to some extent +# largest errors for n=m or n= m+1 (could be fixed?) +""" + function testZZLegendreP2F() + compare ZZLegendreP2 with Float64 SLegendreP2 +""" +function testZZLegendreP2F() + x=.6 + rerr = 6e-7 + for n=1:50 + for m= 0:n #rand(1:n) + a=ZZLegendreP2(n,m,x) + b=Float64(SLegendreP2(n,m,big(x))) + df=norm(a-b)/norm(b) + if ( df > rerr) + println(n," ", m," ", a," ",b," ",df) + end + end + end +end +testZZLegendreP2F() + + +# derived from work published at +# i. a. selezneva, yu. l. ratis, e. hernández, j. pérez-quiles +#and p. fernández de córdoba: +#a code to calculate high order legendre polynomials +#rev. acad. colomb. cienc.: volumen xxxvii, número 145 - diciembre 2013 +#www.scielo.org.co/pdf/racefn/v37n145/v37n145a09.pdf + #Selezneva.....pdf + + +""" +ZZLegendreP3(n::Integer,m::Integer,z::Number) + associated Legendre function First kind type 3 ( z > 1.) + uses different three term recursion +""" + +function ZZLegendreP3(n::Integer,m::Integer,z::Number) + + nz=1 + pnm = zeros(typeof(z),2*n+1)#was nz,n+1 + fac = prod(2.:n) + sqz2 = sqrt((z.*z - one(z))) + hsqz2 = 0.5*sqz2 + ihsqz2 = z./hsqz2 + if(n==0) + pnm[1]=one(z) + return ((-one(z))^m)*pnm[n+1+m] + end + if(n==1) + pnm[1]=-.5*sqz2 + pnm[2]=z + pnm[3]=-sqz2 # sign + return ((-one(z))^m)*pnm[n+1+m] + end + pnm[1] = (1-2*abs(n-2*floor(n/2)))*hsqz2.^n/fac + pnm[2] = -pnm[1]*n.*ihsqz2 + for mr=1:2*n-1 + pnm[mr+2]=(mr-n).*ihsqz2.*pnm[mr+1]+(2*n-mr+1)*mr*pnm[mr] + + end + return ((-one(z))^m)*pnm[n+1+m] +end + + + +#checking ZZLegendreP3 against Float64(SLegendreP3) +# for n< 18 the relative error is less than 4e-6 +# the error depends on the n amd m values and x values to some extent +#n=1 may have error +""" + function testZZLegendreP3F() + compare ZZLegendreP@ with float64 SLegendreP3 +""" +function testZZLegendreP3F() + x=2. + rerr = 4e-6 + for n=0:17 + for m=0:n #rand(1:n) + a=ZZLegendreP3(n,m,x) + b=Float64(SLegendreP3(n,m,big(x))) + df=norm(a-b)/norm(b) + if ( df > rerr) + println(n," ", m," ", a," ",b," ",df) + end + end + end +end +testZZLegendreP3F() + + From f41bf17217ffc4290d065acea9c57da78c7abb35 Mon Sep 17 00:00:00 2001 From: elaineVRC Date: Thu, 17 Aug 2017 21:37:25 -0700 Subject: [PATCH 5/7] < As requested, @btime from BenchmarkTools and @time both gave similar results in that ZZLegendreP2 is light years ahead of scipy.special.lpmv. Accuracy of ZZLegendreP2 ( |z|<=1) against lpmv and ZZLegendreP3 (z > 1) against Fortran table is good. --- associatedLegendre.jl | 883 ++++++++++++++++++++++++++++++------------ 1 file changed, 642 insertions(+), 241 deletions(-) diff --git a/associatedLegendre.jl b/associatedLegendre.jl index d015035f..d4e39031 100644 --- a/associatedLegendre.jl +++ b/associatedLegendre.jl @@ -1,23 +1,28 @@ # solutions, w, of the differential equation are # associated Legendre functions -# (1-x*x)d/dx(dw/dx)⁢-2⁢x⁢(dw/dx)+(n(n+1)-m*m/(1-x*x))w=0 +# (1-x*x)d/dx(dw/dx)⁢-2⁢x⁢(dw/dx)+(n(n+1)-m*m/(1-x*x))w=0 +#n is the degree and m is the order. #for m =0 they are a.k.a. Legendre polynomials -#n is the degree and m is the order +# see https://github.com/pjabardo/Jacobi.jl/blob/master/src/legendre.jl +# for m=0 also. # -#There are more than a handful of programs. Some fast and inaccurate -#and some slow and more accurate; pick your poison (trade off) . +# SSLegendreP2 and 3 : The finite sum I derived is arranged so that +# each term is a multiple of the preceeding term. # -# SSLegendreP2 and 3 : The finite sum is arranged so that -# each term is a multiple of the preceeding term. SLegendreP2 and 3 -# are very slow and used for testing. -# -# The recursion equation 8.5.3 Abramowitz & Stegun (AS) Handbook of -# Mathematical Functions used in LegendreP2 is +# LegendreP2 uses the recursion equation 8.5.3 Abramowitz & Stegun (AS) +# https://www.math.hkbu.edu.hk/support/aands/frameindex.htm +#Handbook of Mathematical Functions is # (n-m+1)P(n+1,m,z)=(2n+1)zP(n,m,z)- (n+m)P(n-1,m,z) +# By normalizing so that that P*P integrated is unity, the new P denoted +# by PP satisfies an altered 3 term recursion. This is discussed in +# Numerical Recipes 3rd ed. see https://github.com/milthorpe/ +# SphericalHarmonics.jl/blob/master/src/SphericalHarmonics.jl +# T.Fukushima has Fortran codes (https://static-content.springer. +# com/esm/art%3A10.1007%2Fs00190-011-0519-2/MediaObjects +# /190_2011_519_MOESM_ESM.txt);JamesBremerJr(github) and arXiv:1707.03287 # -# -# ZZLegendreP2 and 3 uses a different recursion +# ZZLegendreP2 and 3 uses a different recursion # P(n,m+2,x)+2(m+1)x(s(1-x*x))^(-1/2)P(n,m+1,x)+s(n-m)(n+m+1)P(n,m,x)=0 #from http://dlmf.nist.gov/14.10 equation 14.10.1 where s=sign(1-x*x) #refer to Lebedev, Special Functions, Dover Pub. eqn. 7.12.8 page 194 @@ -31,7 +36,7 @@ #Thus we define functions LegendreP2,....P3,....Q2,...Q3 # for n and m positive integers 0<=m<=n # refer to http://dlmf.nist.gov/14.6 -#P(n,m,z) and Q(n,m,⁡z) (z> 1) are often referred to as the prolate +#P(n,m,z) and Q(n,m,⁡z) (z> 1) are often referred to as the prolate #spheroidal harmonics of the first and second kinds, respectively # We can define spherical harmonics (there are other definitions) # Y(n,m,theta,phi)= LegendreP2(n,|m|,theta)*exp(im*m*phi)*((-)^(|m|-m)/2)) @@ -40,9 +45,6 @@ - - - """ LegendreP2(n::Integer,m::Integer,x::Number) associated Legendre function of first kind type 2 (-1.<=x<=1.). @@ -57,11 +59,8 @@ function LegendreP2(n::Integer,m::Integer,x::Number) # -1<=x <= 1 M = (2*m -1) # M must be odd # (2m-1)!!= (2m)!/( m! 2^m) double factorial - dblfac=1 - for j=1:M - if iseven(j) - continue - end + dblfac=one(x) + for j=1:2:M dblfac=j*dblfac end @@ -85,11 +84,118 @@ end +using BenchmarkTools + +using PyCall +pyimport_conda("scipy.special","scipy") #### +@pyimport scipy.special as s + + +function timelpmv() + for k=1:10^6 + k2=k-1 + x=.1 +k2*.9e-6 + s.lpmv(15,30,x) + end +end + +@btime timelpmv() + +function timeLegendreP2() + for k=1:10^6 + k2=k-1 + x2= .1 + k2*.9e-6 + LegendreP2(30,15,x2) + end +end + +@btime timeLegendreP2() + + +""" + SSLegendreP2(n::Integer,m::Integer,x::Number) + associated Legendre function 1st kind type 2 , -1.<=x<=1. + each term in sum is multiple of previous, from SLegendreP2 + + +""" + +function SSLegendreP2(n::Integer,m::Integer,x::Number)#n,m integers + + M = (2*n -1) # M must be odd + # (2m-1)!!= (2m)!/( m! 2^m) double factorial + + DT=one(x) #1.0 + if n > 0 + for j=1:M + if isodd(j) + DT=DT*j + end + if j <= (n-m) + DT=DT/j + end + end + end + MXP=div(m+n,2) + + sum= (x^(n+m))*DT/((one(x)-x*x)^(m/2)) + prev=sum + for p=1:MXP + term=(-prev)*(n-p+1)*(n+m-2*p+1)*(n+m-2*p+2 + )/(p*(2*n-2*p+1)*(2*n-2*p+2)*x*x) + sum=sum+term + prev=term + end + return sum +end + + + +""" + ZZLegendreP2(n::Integer,m::Integer,z::Number) + associated Legendre function First kind type 2 |x|<=1. + uses different three term recursion +""" + +function ZZLegendreP2(n::Integer,m::Integer,z::Number) + nz=1 + pnm = zeros(typeof(z),2*n+1) + #fac = prod(2.:n) + fac= one(z) #1.0 + for k=1:n + fac=fac*k + end + sqz2 = sqrt((one(z)-z.*z)) + hsqz2 = 0.5*sqz2 + ihsqz2 = z./hsqz2 + if(n==0) + pnm[1]=one(z) + return ((-one(z))^m)*pnm[n+1+m] #pnm + end + if(n==1) + pnm[1]=-.5*sqz2 + pnm[2]=z + pnm[3]=sqz2 + return ((-one(z))^m)*pnm[n+1+m] #pnm + end + pnm[1] = (1-2*abs(n-2*floor(n/2)))*hsqz2.^n/fac + pnm[2] = -pnm[1]*n.*ihsqz2 + for mr=1:2*n-1 + pnm[mr+2]=(mr-n).*ihsqz2.*pnm[mr+1]-(2*n-mr+1)*mr*pnm[mr] + end + return ((-one(z))^m)*pnm[n+1+m] +end + + #checking LegendreP2 against scipy.special # for n< 18 the relative error is less than 2e-14 -# the error depends on the n amd m values and x values to some extent +# for n < 45 rel error less than 2e-13 +# for n < 62 the rel.error is less than 5e-12 +# for n < 75 the rel.error is less than 5e-11 +# the error depends on the n amd m values and x values to s(ome extent using PyCall +pyimport_conda("scipy.special","scipy") #### @pyimport scipy.special as s @@ -101,8 +207,8 @@ using PyCall function testLegendreP2() x=.6 - rerr = 2e-14 - for n=0:17 + rerr = 5e-11 + for n=0:75 #17 for m=0:n #m=rand(0:n) a=LegendreP2(n,m,x) b=s.lpmv(m,n,x) @@ -120,104 +226,70 @@ end testLegendreP2() -""" - SLegendreP2(n,m,x) - Legendre function 1st kind type 2 , -1.<=x<=1. - series (to be optimized for in SSLegendreP2) -""" - -function SLegendreP2(n,m,x)#n,m integers - - v=(big(1.)-x*x)^(-m*big(1)/2) - v=v*factorial(big(n+m))/(((big(2))^n)* factorial(big(n))); - term = big(0.); - #if iseven(n+m); MXP=Int(trunc((n+m)/2));end; - #if isodd(n+m);MXP=Int(trunc((n+m-1)/2));end; - #MXP=Int(floor((n+m)/2)); - - MXP=Int(trunc(((n+m)/2))) - - - for p=0:MXP - - term= term + - ((-big(1.))^p)*binomial(big(2*(n-p)),big(n-m))* - (x^(big(1)*(n + m -2*p)))* - binomial(big(n),big(p)) - - end - return term * v -end - -# I derived this finite sum using 8.6.6,8.6.18,8.2.5 (AS), -#binomial expansion and differentiating - - - - """ function testLegendreP2F() - compare LegendreP2 against Float64(SLegendreP2) + compare LegendreP2 against Float64(SSLegendreP2) """ function testLegendreP2F() #= - checking LegendreP2 against Float64SLegendreP2 + checking LegendreP2 against Float64SSLegendreP2 for n < 18 the relative error is less than 7e-15 - the error depends on the n amd m values and + error depends on the n amd m values and x values to some extent =# x=.6 rerr = 7e-15 - for n=1:17 + for n=1:75 + if n > 18 + rerr=2e-13 + end + + + if n > 44 + rerr= 5e-12 + end + if n > 62 + rerr= 5e-11 + end for m=0:n #rand(1:n) a=LegendreP2(n,m,x) - b=Float64(SLegendreP2(n,m,big(x))) + b=Float64(SSLegendreP2(n,m,big(x))) df=norm(a-b)/norm(b) if ( df > rerr) println(n," ", m," ", a," ",b," ",df) end end -end + end + end testLegendreP2F() - - +#checking SSLegendreP2 against Float64(SSLegendreP2) +# for n< 17 the relative error is less than 3e-11 +# the error depends on the n amd m values and x values to some extent """ - SSLegendreP2(n::Integer,m::Integer,x::Number) -associated Legendre function 1st kind type 2 , -1.<=x<=1. -each term in sum is multiple of previous, from SLegendreP2 - - + function testSSLegendreP2F() + compare SSLegendreP2 to float64SSLegendreP2 """ - -function SSLegendreP2(n::Integer,m::Integer,x::Number)#n,m integers - - M = (2*n -1) # M must be odd - # (2m-1)!!= (2m)!/( m! 2^m) double factorial - dblfac=1 - for j=1:M - if iseven(j) - continue - end - dblfac=j*dblfac - end - - MXP=Int(trunc(((n+m)/2))) - sum= (x^(n+m))*dblfac/(((one(x)-x*x)^(m/2))*prod(1:n-m)) - prev=sum - for p=1:MXP - term=(-prev)*(n-p+1)*(n+m-2*p+1)*(n+m-2*p+2 - )/(p*(2*n-2*p+1)*(2*n-2*p+2)*x*x) - sum=sum+term - prev=term +function testSSLegendreP2F() + x=.6 + rerr = 3e-11 + for n=0:17 + for m=0:n #m=rand(0:n) + a=SSLegendreP2(n,m,x) + b=Float64(SSLegendreP2(n,m,big(x))) + df=norm(a-b)/norm(b) + if ( df > rerr) + println(n," ", m," ", a," ",b," ",df) + end + end end - return sum end +testSSLegendreP2F() @@ -253,12 +325,12 @@ function testSSLegendreP2() end testSSLegendreP2() -#checking SSLegendreP2 against Float64(SLegendreP2) +#checking SSLegendreP2 against Float64(SSLegendreP2) # for n< 17 the relative error is less than 3e-11 # the error depends on the n amd m values and x values to some extent """ function testSSLegendreP2F() - compare SSLegendreP2 to float64SLegendreP2 + compare SSLegendreP2 to float64SSLegendreP2 """ function testSSLegendreP2F() x=.6 @@ -266,7 +338,7 @@ function testSSLegendreP2F() for n=0:17 for m=0:n #m=rand(0:n) a=SSLegendreP2(n,m,x) - b=Float64(SLegendreP2(n,m,big(x))) + b=Float64(SSLegendreP2(n,m,big(x))) df=norm(a-b)/norm(b) if ( df > rerr) println(n," ", m," ", a," ",b," ",df) @@ -276,65 +348,56 @@ function testSSLegendreP2F() end testSSLegendreP2F() +# derived from work published at +# i. a. selezneva, yu. l. ratis, e. hernández, j. pérez-quiles +#and p. fernández de córdoba: +#a code to calculate high order legendre polynomials +#rev. acad. colomb. cienc.: volumen xxxvii, número 145 - diciembre 2013 +#www.scielo.org.co/pdf/racefn/v37n145/v37n145a09.pdf + #Selezneva.....pdf + """ -SLegendreP3(n::Integer,m::Integer,z::Number) - associated Legendre function z>1 1st kind type 3 - series expansion +ZZLegendreP3(n::Integer,m::Integer,z::Number) + associated Legendre function First kind type 3 ( z > 1.) + uses different three term recursion """ -function SLegendreP3(n::Integer,m::Integer,z::Number) -#function mylegenp3(n,m,z)#n,m integers - v=(z*z - big(1.))^(-m*big(1)/2) - - v=v*factorial(big(n+m))/(((big(2))^n)* factorial(big(n))); - term = big(0.); - #if iseven(n+m); MXP=Int(trunc((n+m)/2));end; - #if isodd(n+m);MXP=Int(trunc((n+m-1)/2));end; - #MXP=Int(floor((n+m)/2)); - MXP=Int(trunc(((n+m)/2))) - - for p=0:MXP - - term= term + - ((-big(1.))^p)*binomial(big(2*(n-p)),big(n-m))* - (z^(big(1)*(n + m -2*p)))* - binomial(big(n),big(p)) - - end - return term * v -end - -#print mylegenp3(2,1,2.)#agrees to 16 decimals -#print mylegenp3(2,1,mpf(2))#agrees -#print mylegenp3(2,1,mpf('2.'))#agrees 10.392304845413263761 -#println(mylegenp3(2,1,big(2.))) AGREES with mpmath - -# I derived this finite sum using 8.6.6,8.6.18,8.2.5 (AS), -#binomial expansion and differentiating - -function timeS2() - for i=1:10 - SLegendreP2(100,50,big(.6)) +function ZZLegendreP3(n::Integer,m::Integer,z::Number) + + nz=1 + pnm = zeros(typeof(z),2*n+1)#was nz,n+1 + #fac = prod(2.:n) + fac= one(z) #1.0 + for k=1:n + fac=fac*k end -end -TS2= @elapsed(timeS2()) -println("time SLegendreP2(100,50,big(.6)) $TS2") - -function timeP2() - for i=1:10 - SSLegendreP2(100,50,big(.6)) + sqz2 = sqrt((z.*z - one(z))) + hsqz2 = 0.5*sqz2 + ihsqz2 = z./hsqz2 + if(n==0) + pnm[1]=one(z) + return ((-one(z))^m)*pnm[n+1+m] end -end -TP2= @elapsed(timeP2()) -println("time SSLegendreP2(100,50,big(.6)) $TP2") -#time SLegendreP2(100,50,big(.6)) 0.0085234 -#time SSLegendreP2(100,50,big(.6)) 0.004890456 -# the series is double speed when each term is a multiple of previous term - + if(n==1) + pnm[1]=-.5*sqz2 + pnm[2]=z + pnm[3]=-sqz2 # sign + return ((-one(z))^m)*pnm[n+1+m] + end + pnm[1] = (1-2*abs(n-2*floor(n/2)))*hsqz2.^n/fac + pnm[2] = -pnm[1]*n.*ihsqz2 + for mr=1:2*n-1 + pnm[mr+2]=(mr-n).*ihsqz2.*pnm[mr+1]+(2*n-mr+1)*mr*pnm[mr] + + end + return ((-one(z))^m)*pnm[n+1+m] +end + -""" + + """ SSLegendreP3(n::Integer,m::Integer,x::Number) associated Legendre function 1st kind type 3 , x >1. each term in sum is multiple of previous. @@ -346,16 +409,21 @@ function SSLegendreP3(n::Integer,m::Integer,x::Number)#n,m integers M = (2*n -1) # M must be odd # (2m-1)!!= (2m)!/( m! 2^m) double factorial - dblfac=1 - for j=1:M - if iseven(j) - continue - end - dblfac=j*dblfac + # dblfac=1 + DT=one(x) #1.0 + if n > 0 + for j=1:M + if isodd(j) + DT=DT*j + end + if j <= (n-m) + DT=DT/j + end + end end - MXP=Int(trunc(((n+m)/2))) - sum= (x^(n+m))*dblfac/(((x*x - one(x))^(m/2))*prod(1:n-m)) + MXP=div(m+n,2) + sum= (x^(n+m))*DT/((x*x - one(x))^(m/2)) prev=sum for p=1:MXP term=(-prev)*(n-p+1)*(n+m-2*p+1)*(n+m-2*p+2 @@ -367,77 +435,37 @@ function SSLegendreP3(n::Integer,m::Integer,x::Number)#n,m integers end - + -#checking SSLegendreP3 against Float64(SLegendreP3) -# for n< 17 the relative error is less than 4e-14 -# the error depends on the n amd m values and x values to some extent -""" - function testSSLegendreP3F() - compare SSLegendreP3 with Float64 SLegendreP3 -""" -function testSSLegendreP3F() - x=2. - rerr = 4e-14 - for n=0:17 - for m=rand(0:n) - a=SSLegendreP3(n,m,x) - b=Float64(SLegendreP3(n,m,big(x))) - df=norm(a-b)/norm(b) - if ( df > rerr) - println(n," ", m," ", a," ",b," ",df) - end - end -end -end -testSSLegendreP3F() -# based on code published at -# i. a. selezneva, yu. l. ratis, e. hernández, j. pérez-quiles -#and p. fernández de córdoba: -#a code to calculate high order legendre polynomials -#rev. acad. colomb. cienc.: volumen xxxvii, número 145 - diciembre 2013 -#www.scielo.org.co/pdf/racefn/v37n145/v37n145a09.pdf - #Selezneva.....pdf - -""" -ZZLegendreP2(n::Integer,m::Integer,z::Number) - associated Legendre function First kind type 2 |x|<=1. - uses different three term recursion -""" -function ZZLegendreP2(n::Integer,m::Integer,z::Number) - nz=1 - pnm = zeros(typeof(z),2*n+1) - fac = prod(2.:n) - sqz2 = sqrt((one(z)-z.*z)) - hsqz2 = 0.5*sqz2 - ihsqz2 = z./hsqz2 - if(n==0) - pnm[1]=one(z) - return ((-one(z))^m)*pnm[n+1+m] #pnm - end - if(n==1) - pnm[1]=-.5*sqz2 - pnm[2]=z - pnm[3]=sqz2 - return ((-one(z))^m)*pnm[n+1+m] #pnm - end - pnm[1] = (1-2*abs(n-2*floor(n/2)))*hsqz2.^n/fac - pnm[2] = -pnm[1]*n.*ihsqz2 - for mr=1:2*n-1 - pnm[mr+2]=(mr-n).*ihsqz2.*pnm[mr+1]-(2*n-mr+1)*mr*pnm[mr] - end - return ((-one(z))^m)*pnm[n+1+m] -end + +#A representation of associated Legendre function was derived from 8.6.6, +# 8.6.18,8.2.5,3.1.1 in AS (https://www.math.hkbu.edu.hk +# /support/aands/frameindex.htm). +# 14.7.8,14.7.10,14.7.11,14.7.14.1.2.2 in DLMF(dlmf.nist.gov), +#similar formula for type 2. +#(8.6.18 AS) Rodrigues' formula for integer n +# P(n,z)=(1/((2^n)n!)) (d/dz)^n (z^2 - 1)^n +#(8.6.6 AS) P(n,m,z)= ((z^2 -1)^(m/2))(d/dz)^m P(n,z) for z > 1 +# P(n,m,x) = ((-)^m)((1-x^2)^(m/2))(d/dz)^m P(n,x) for |x| <= 1 +# expand using binomial theorem where a=z^2 and b = -1 +#(3.1.1 AS) define binomial coefficient = C(n,k) = n!/((n-k)! k!) +# defining n!=n*(n-1)*(n-2)*...3*2*1 = factorial(n) +# (a + b)^n = a^n + C(n,1)*(a^(n-1))*b + C(n,2)*(a^(n-1))*b^2 +# +C(n,3)*(a^(n-3))*b^3 +.......+b^n where n is a positive integer +#differentiate to obtain the finite sum expression + + + #Benchmarking ZZLegendreP2 against scipy.special using PyCall - +using BenchmarkTools @pyimport scipy.special as s @@ -475,18 +503,67 @@ println("time ZZLegendreP2(30,15,.6) $TZ2") end BenchmarkZZLegendreP2() +function timeZZ2() + for k=1:10^6 + k2=k-1 + x2=.1 +k2*.9e-6 + ZZLegendreP2(30,15,x2) + end +end + +function LPMV2() + for k=1:10^6 + k2=k-1 + x2=.1 +k2*.9e-6 + s.lpmv(15,30,x2) + end +end + +println("lpmv time") +@btime LPMV2() +println("ZZLegendreP2 time") +@btime timeZZ2() + +#= + for k=1:10^6 + k2=k-1 + x2=.1 +k2*.9e-6 + s.lpmv(15,30,x2) + end + + +println("lpmv time") +@btime LPMV2() +println("ZZLegendreP2 time") +@btime timeZZ2() + +=# + + + + + + + + + #checking ZZLegendreP2 against scipy special -# for n< 51 the relative error is less than 6e-7 +# for n< 41 the relative error is less than 2.4e-9 +# for n < 61 the relative error is less than 6e-5 # the error depends on the n amd m values and x values to some extent -# +# +using PyCall +@pyimport scipy.special as s + """ - function testZZLegendreP2() -compare ZZLegendreP2 with scipy special + function testZZLegendreP2( + compare ZZLegendreP2 with scipy special """ function testZZLegendreP2() x=.6 - rerr = 6e-7 - for n=1:50 + rerr = 2.4e-9# 5e-7 + #for n =51: 80 + for n=1:40 for m= 0:n #rand(1:n) a=ZZLegendreP2(n,m,x) b=s.lpmv(m,n,x) @@ -500,21 +577,22 @@ end testZZLegendreP2() -#checking ZZLegendreP2 against Float64(SLegendreP2) -# for n< 51 the relative error is less than 6e-7 +#checking ZZLegendreP2 against Float64(SSLegendreP2) +# for n< 41 the relative error is less than 2.4e-9 +# for n <61 the relative error is less than 6e-5 # the error depends on the n amd m values and x values to some extent # largest errors for n=m or n= m+1 (could be fixed?) """ function testZZLegendreP2F() - compare ZZLegendreP2 with Float64 SLegendreP2 + compare ZZLegendreP2 with Float64 SSLegendreP2 """ function testZZLegendreP2F() x=.6 - rerr = 6e-7 - for n=1:50 + rerr = 2.4e-9 + for n=1:40 for m= 0:n #rand(1:n) a=ZZLegendreP2(n,m,x) - b=Float64(SLegendreP2(n,m,big(x))) + b=Float64(SSLegendreP2(n,m,big(x))) df=norm(a-b)/norm(b) if ( df > rerr) println(n," ", m," ", a," ",b," ",df) @@ -526,16 +604,16 @@ testZZLegendreP2F() # derived from work published at -# i. a. selezneva, yu. l. ratis, e. hernández, j. pérez-quiles -#and p. fernández de córdoba: +# i. a. selezneva, yu. l. ratis, e. hernández, j. pérez-quiles +#and p. fernández de córdoba: #a code to calculate high order legendre polynomials -#rev. acad. colomb. cienc.: volumen xxxvii, número 145 - diciembre 2013 +#rev. acad. colomb. cienc.: volumen xxxvii, número 145 - diciembre 2013 #www.scielo.org.co/pdf/racefn/v37n145/v37n145a09.pdf #Selezneva.....pdf """ -ZZLegendreP3(n::Integer,m::Integer,z::Number) + ZZLegendreP3(n::Integer,m::Integer,z::Number) associated Legendre function First kind type 3 ( z > 1.) uses different three term recursion """ @@ -544,7 +622,11 @@ function ZZLegendreP3(n::Integer,m::Integer,z::Number) nz=1 pnm = zeros(typeof(z),2*n+1)#was nz,n+1 - fac = prod(2.:n) + #fac = prod(2.:n) + fac= one(z) #1.0 + for k=1:n + fac=fac*k + end sqz2 = sqrt((z.*z - one(z))) hsqz2 = 0.5*sqz2 ihsqz2 = z./hsqz2 @@ -569,21 +651,25 @@ end -#checking ZZLegendreP3 against Float64(SLegendreP3) -# for n< 18 the relative error is less than 4e-6 + + +#checking ZZLegendreP3 against Float64(SSLegendreP3) +# for n< 17 the relative error is less than 6.4e-7 for z=2. +# for n<24 the relative error is less than .0054 +# for z=3. n<17 the relative eror is less than 2e-7 # the error depends on the n amd m values and x values to some extent #n=1 may have error """ function testZZLegendreP3F() - compare ZZLegendreP@ with float64 SLegendreP3 + compare ZZLegendreP3 with float64 SSLegendreP3 """ function testZZLegendreP3F() - x=2. - rerr = 4e-6 - for n=0:17 + x=3.#2. + rerr = 2e-7 #6.4e-7 + for n=1:16 for m=0:n #rand(1:n) a=ZZLegendreP3(n,m,x) - b=Float64(SLegendreP3(n,m,big(x))) + b=Float64(SSLegendreP3(n,m,big(x))) df=norm(a-b)/norm(b) if ( df > rerr) println(n," ", m," ", a," ",b," ",df) @@ -594,3 +680,318 @@ end testZZLegendreP3F() + + + + +using PyCall +@pyimport scipy.special as s + +""" + function ZZSS2(n,m,x,ESTOP) + Values near .99999 < x <= 1.000000 and -1.000000 ESTOP + println() + println("n=$n, m=$m , x=$x , zz=$valuezz") + println("SS=$valueSS,ZZ=$valueZZ, py=$valuepy ") + println(" zzrerr= $zzrerr, prerr= $prerr") + end + +end + + +ESTOP=.000000000001 + +n=10 +for ix=1:10 + x=1. - (.1)^ix + for m=0:n + ZZSS2(n,m,x,ESTOP) + end +end + +n=10 +for ix=1:10 + x=-1. + (.1)^ix + for m=0:n + ZZSS2(n,m,x,ESTOP) + end +end + + + + + + + + + + + + +""" + function SSZZ1(n,m,x,ESTOP) + Values near 1.< x < 1.00001 causing problems + Identify n,m, x ,when zz goes negative, its underflow + The cutoff does change with n values + The user has the option to accept reduced accuracy + for certain values of n,m,x or to increase the accuracy + by using higher precision +""" +function SSZZ1(n,m,x,ESTOP) + valueZZ=Float64(ZZLegendreP3(n,m,big(x)))#high precision + #valueS=Float64(SLegendreP3(n,m,big(x))) + valueSS=Float64(SSLegendreP3(n,m,big(x))) + valuezz=ZZLegendreP3(n,m,x) + #valuess=SSLegendreP3(n,m,x) + zzrerr=abs((valuezz-valueZZ)/valueZZ) + #ssrerr=abs((valuess-valueSS)/valueSS) + + if zzrerr > ESTOP + #if sign(valuezz) < 1 + println() + println("n= $n, m= $m, x= $x") + println("SS=$valueSS,ZZ=$valueZZ, zz= $valuezz, zzrerr= $zzrerr") + + + #println("S=$valueS,SS=$valueSS,ZZ=$valueZZ, n= $n, m= $m, x= $x") + #println(" ssrerr= $ssrerr,ss= $valuess,zz= $valuezz, zzrerr= $zzrerr") + + end +end + + + + + +ESTOP = .000000000001 #e-12 + +n=10 + +for ix =1:7 + for m=0:n + x=1.+.1^(ix) + + SSZZ1(n,m,x,ESTOP) + end +end + +println(" ") + + n=80 +ESTOP=.000000000001 # e-12 +for ix=1:7 + for m=0:n + x= 1.+ (.1)^(ix) + + SSZZ1(n,m,x,ESTOP) + end +end + + + + + +""" + function TableSSZZ(n,m,x,T) + compare SSLegendreP3,ZZLegendreP3 and Table1 value is T + with and without arbitrary precision to + Table 1 from Computer Physics Comm. 181 (2010)2091-7 + User has the option to accept reduced accuracy for certain + values of n,m,x or to increase the accuracy by using + higher values of precision in arbitrary precision arithmetic +""" +function TableSSZZ(n,m,x,T) + valueZZ=Float64(ZZLegendreP3(n,m,big(x))) #higher precision + #valueS=Float64(SLegendreP3(n,m,big(x))) + valueSS=Float64(SSLegendreP3(n,m,big(x))) + valuezz=ZZLegendreP3(n,m,x) + #valuess=SSLegendreP3(n,m,x) + zzrerr=abs((valuezz-T)/T) + ZZrerr=abs((valueZZ-T)/T) + println() + println("n=$n, m=$m, x=$x") + println("SS=$valueSS,ZZ=$valueZZ, T=$T ,zz=$valuezz " ) + # println("SS=$valueSS,ZZ=$valueZZ, T=$T , n= $n, m= $m, x= $x") + # println(" ZZrerr=$ZZrerr, ss=$valuess,zz=$valuezz, zzrerr=$zzrerr") + println(" ZZrerr=$ZZrerr, zzrerr=$zzrerr") + +end +x=5.0 +m=0 +n=0 +T=1.0 +x=5.0 +m=0 +n=0 +T=1.0 +TableSSZZ(n,m,x,T) +n=1 +T=5.0 +TableSSZZ(n,m,x,T) +n=2 +T=37.0 +TableSSZZ(n,m,x,T) +n=3 +T=305. +TableSSZZ(n,m,x,T) +n=4 +T=2641. +TableSSZZ(n,m,x,T) +m=1 +n=1 +T=4.898979485566 +TableSSZZ(n,m,x,T) +n=2 +T=73.4846922835 +TableSSZZ(n,m,x,T) +n=3 +T=9.112101843153e2 +TableSSZZ(n,m,x,T) +n=4 +T=1.053280589397e4 +TableSSZZ(n,m,x,T) +m=2 +n=2 +T=72. +TableSSZZ(n,m,x,T) +n=3 +T=1800. +TableSSZZ(n,m,x,T) +n=4 +T=3.132e4 +TableSSZZ(n,m,x,T) +m=3 +n=3 +T=1.763632614804e3 +TableSSZZ(n,m,x,T) +n=4 +T=6.172714151814e4 +TableSSZZ(n,m,x,T) +m=4 +n=4 +T=6.048e4 +TableSSZZ(n,m,x,T) +x=2.0 +m=0 +n=0 +T=1.0 +TableSSZZ(n,m,x,T) +n=1 +T=2.0 +TableSSZZ(n,m,x,T) +n=2 +T=5.5 +TableSSZZ(n,m,x,T) +n=3 +T=17. +TableSSZZ(n,m,x,T) +n=4 +T=55.375 +TableSSZZ(n,m,x,T) +m=1 +n=1 +T=1.732050807569 +TableSSZZ(n,m,x,T) +n=2 +T=1.039230484541e1 +TableSSZZ(n,m,x,T) +n=3 +T=4.936344801571e1 +TableSSZZ(n,m,x,T) +n=4 +T=2.165063509461e2 +TableSSZZ(n,m,x,T) +m=2 +n=2 +T=9.0 +TableSSZZ(n,m,x,T) +n=3 +T=90. +TableSSZZ(n,m,x,T) +n=4 +T=607.5 +TableSSZZ(n,m,x,T) +m=3 +n=3 +T=7.794228634060e1 +TableSSZZ(n,m,x,T) +n=4 +T=1.091192008768e3 +TableSSZZ(n,m,x,T) +m=4 +n=4 +T=945.0 +TableSSZZ(n,m,x,T) +x=1.000001 +m=0 +n=0 +T=1.0 +TableSSZZ(n,m,x,T) +n=1 +T=1.0000010 +TableSSZZ(n,m,x,T) +n=2 +T=1.000003000001 +TableSSZZ(n,m,x,T) +n=3 +T=1.000006000007 +TableSSZZ(n,m,x,T) +n=4 +T=1.000010000022 +TableSSZZ(n,m,x,T) +m=1 +n=1 +T=1.414213915900e-3 +TableSSZZ(n,m,x,T) +n=2 +T=4.242645990341e-3 +TableSSZZ(n,m,x,T) +n=3 +T=8.485304708618e-3 +TableSSZZ(n,m,x,T) +n=4 +T=1.414220279870e-2 +TableSSZZ(n,m,x,T) +m=2 +n=2 +T=6.000002999773e-6 +TableSSZZ(n,m,x,T) +n=3 +T=3.000004499888e-5 +TableSSZZ(n,m,x,T) +n=4 +T=9.000025499681e-5 +TableSSZZ(n,m,x,T) +m=3 +n=3 +T=4.242643868860e-8 +TableSSZZ(n,m,x,T) +n=4 #underflow +T=2.969853678052e-7 +TableSSZZ(n,m,x,T) +m=4 #underflow +n=4 +T=4.200004199683e-10 +TableSSZZ(n,m,x,T) + + + + + + + + From 8362ec84e752336a32e831f0b9e053c53369736a Mon Sep 17 00:00:00 2001 From: elaineVRC Date: Thu, 24 Aug 2017 18:30:06 -0700 Subject: [PATCH 6/7] Update associatedLegendre.jl This PR is the file entitled associatedLegendre.jl. The first part with all the comments is like a readme. The rest is the source code and numerous tests(requested) and exploration of area near |x| ~ 1. As requested @btime from BenchmarkTools as well as @time gave similar results in that ZZLegendreP2 is light years ahead of scipy.special.lpmv. Accuracy of ZZLegendreP2 (|z|=<1.) against lpmv and ZZLegendreP3 (z>1.) against Fortran is good. --- associatedLegendre.jl | 30 ++++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/associatedLegendre.jl b/associatedLegendre.jl index d4e39031..e2859ac2 100644 --- a/associatedLegendre.jl +++ b/associatedLegendre.jl @@ -1,7 +1,7 @@ # solutions, w, of the differential equation are # associated Legendre functions -# (1-x*x)d/dx(dw/dx)⁢-2⁢x⁢(dw/dx)+(n(n+1)-m*m/(1-x*x))w=0 +# (1-x*x)d/dx(dw/dx)⁢-2⁢x⁢(dw/dx)+(n(n+1)-m*m/(1-x*x))w=0 #n is the degree and m is the order. #for m =0 they are a.k.a. Legendre polynomials # see https://github.com/pjabardo/Jacobi.jl/blob/master/src/legendre.jl @@ -36,7 +36,7 @@ #Thus we define functions LegendreP2,....P3,....Q2,...Q3 # for n and m positive integers 0<=m<=n # refer to http://dlmf.nist.gov/14.6 -#P(n,m,z) and Q(n,m,⁡z) (z> 1) are often referred to as the prolate +#P(n,m,z) and Q(n,m,⁡z) (z> 1) are often referred to as the prolate #spheroidal harmonics of the first and second kinds, respectively # We can define spherical harmonics (there are other definitions) # Y(n,m,theta,phi)= LegendreP2(n,|m|,theta)*exp(im*m*phi)*((-)^(|m|-m)/2)) @@ -65,7 +65,7 @@ function LegendreP2(n::Integer,m::Integer,x::Number) end pj2= ((-1)^m)*dblfac*(one(x)-x*x)^(m/2) - pj1=x*(2.*m+1.)*pj2 + pj1=x*(2*m+1)*pj2 # pj1=x*(2.*m+1.)*pj2 if n == m return pj2 elseif n == m+1 @@ -74,7 +74,8 @@ function LegendreP2(n::Integer,m::Integer,x::Number) for j = m+2 :n - pjj=(x*(2.*j-1.)*pj1 - pj2*(j +m-1.)) /(j-m) + pjj=(x*(2*j-1)*pj1 - pj2*(j +m-1)) /(j-m) + #pjj=(x*(2.*j-1.)*pj1 - pj2*(j +m-1.)) /(j-m) pj2=pj1 pj1=pjj end @@ -166,14 +167,14 @@ function ZZLegendreP2(n::Integer,m::Integer,z::Number) fac=fac*k end sqz2 = sqrt((one(z)-z.*z)) - hsqz2 = 0.5*sqz2 + hsqz2 = sqz2 /2 # 0.5*sqz2 ihsqz2 = z./hsqz2 if(n==0) pnm[1]=one(z) return ((-one(z))^m)*pnm[n+1+m] #pnm end if(n==1) - pnm[1]=-.5*sqz2 + pnm[1]= -sqz2/2 # -.5*sqz2 pnm[2]=z pnm[3]=sqz2 return ((-one(z))^m)*pnm[n+1+m] #pnm @@ -349,10 +350,10 @@ end testSSLegendreP2F() # derived from work published at -# i. a. selezneva, yu. l. ratis, e. hernández, j. pérez-quiles -#and p. fernández de córdoba: +# i. a. selezneva, yu. l. ratis, e. hernández, j. pérez-quiles +#and p. fernández de córdoba: #a code to calculate high order legendre polynomials -#rev. acad. colomb. cienc.: volumen xxxvii, número 145 - diciembre 2013 +#rev. acad. colomb. cienc.: volumen xxxvii, número 145 - diciembre 2013 #www.scielo.org.co/pdf/racefn/v37n145/v37n145a09.pdf #Selezneva.....pdf @@ -373,14 +374,14 @@ function ZZLegendreP3(n::Integer,m::Integer,z::Number) fac=fac*k end sqz2 = sqrt((z.*z - one(z))) - hsqz2 = 0.5*sqz2 + hsqz2 = sqz2 / 2 # 0.5*sqz2 ihsqz2 = z./hsqz2 if(n==0) pnm[1]=one(z) return ((-one(z))^m)*pnm[n+1+m] end if(n==1) - pnm[1]=-.5*sqz2 + pnm[1]= - sqz2 / 2 # -.5*sqz2 pnm[2]=z pnm[3]=-sqz2 # sign return ((-one(z))^m)*pnm[n+1+m] @@ -604,10 +605,10 @@ testZZLegendreP2F() # derived from work published at -# i. a. selezneva, yu. l. ratis, e. hernández, j. pérez-quiles -#and p. fernández de córdoba: +# i. a. selezneva, yu. l. ratis, e. hernández, j. pérez-quiles +#and p. fernández de córdoba: #a code to calculate high order legendre polynomials -#rev. acad. colomb. cienc.: volumen xxxvii, número 145 - diciembre 2013 +#rev. acad. colomb. cienc.: volumen xxxvii, número 145 - diciembre 2013 #www.scielo.org.co/pdf/racefn/v37n145/v37n145a09.pdf #Selezneva.....pdf @@ -995,3 +996,4 @@ TableSSZZ(n,m,x,T) + From 10d485f92ec01942ccd895b665f6e0e09b9b313e Mon Sep 17 00:00:00 2001 From: elaineVRC Date: Thu, 24 Aug 2017 18:35:59 -0700 Subject: [PATCH 7/7] Update recNM3.ipynb delete --- recNM3.ipynb | 59 +--------------------------------------------------- 1 file changed, 1 insertion(+), 58 deletions(-) diff --git a/recNM3.ipynb b/recNM3.ipynb index 7114b8ea..8b137891 100644 --- a/recNM3.ipynb +++ b/recNM3.ipynb @@ -1,58 +1 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "function recNM3(n,m,z) #n,m positive integers, z number \n", - " #reliable for values of n <= 17 for z=2. \n", - " # 0 <= m <= n\n", - " # z > 1.\n", - " M = (2*m -1) # M must be odd\n", - " # (2m-1)!!= (2m)!/( m! 2^m) double factorial \n", - " dblfac=1\n", - " for j=1:M\n", - " if iseven(j)\n", - " continue\n", - " end\n", - " dblfac=j*dblfac \n", - " end\n", - " if n == m\n", - " return( dblfac*(z*z -1)^(m/2)) \n", - " elseif n == m+1\n", - " return((2.*m +1.)*z*dblfac*(z*z -1)^(m/2))\n", - " end\n", - " pj2=dblfac*(z*z -1.)^(m/2)\n", - " pj1=z*(2.*m+1.)*pj2 \n", - " for j = m+2 :n \n", - " pjj=(z*(2.*j-1.)*pj1 - pj2*(j +m-1.)) /(j-m)\n", - " pj2=pj1\n", - " pj1=pjj\n", - " end\n", - "return pj1 \n", - "end\n", - "\n", - "\n" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Julia 0.4.2", - "language": "julia", - "name": "julia-0.4" - }, - "language_info": { - "file_extension": ".jl", - "mimetype": "application/julia", - "name": "julia", - "version": "0.4.2" - } - }, - "nbformat": 4, - "nbformat_minor": 0 -} +