Artifact 00085d3d56af1d7aebbf13308c2caa710cb951fcffd9f151bda69521dc651879:
- Executable file
r37/packages/specfn/sfother.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 5835) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/specfn/sfother.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 5835) [annotate] [blame] [check-ins using]
module sfother; % Rulesets for the Struve H and L functions, Lommel % 1 and 2 functions and Whittaker M and W functions. % Author: Chris Cannam, Nov 1992. % The aim is to re-express in terms % of other (more `standard') special % functions. No numerical approximation code. % Neither imports nor exports functions. % This module contains only rulesets. algebraic (operator struveh, struvel); algebraic (struve!*rules := { df(struveh(~n,~z),z) => (2/pi) - struveh(1,z) when numberp n and n = 0, df(struveh(~n,~x),x) => (x*StruveH(-1 + n,x)- n*StruveH(n,x))/x, df((z**n)*struveh(~n,~z),z) => (z**n)*struveh(n-1,z), df((z**(-n))*struveh(~n,~z),z) => (1/(sqrt(pi)*(2**n)*gamma(n+(3/2)))) - (z**(-n))*struveh(n+1,z), struveh(~n,~z) => ((-1)**n)*besselj(-n,z) when numberp n and impart n = 0 and n < 0 and (n*2)=floor(n*2) and not evenp floor(n*2), struveh(~n,~z) => ((2/(pi*z))**(1/2))*(1-cos z) when numberp n and n=1/2, struveh(~n,~z) => ((z/(pi*2))**(1/2)) * (1+(2/(z**2))) - ((2/(pi*z))**(1/2)) * (sin z + ((cos z)/z)) when numberp n and n=3/2, struveh(~n,~x) => (x*0.5)^(n+1)*struve_compute_term(n,x,h) when numberp x and numberp n and symbolic !*rounded, struvel(~n,~x) => struve_compute_term(n,x,l) when numberp x and numberp n and symbolic !*rounded, struvel(~n,~z) => besseli(-n,z) when numberp n and impart n = 0 and n < 0 and (n*2)=floor(n*2) and not evenp floor(n*2), struvel(~n,~z) => -i*(e**((-i*n*pi)/2))*struveh(n,i*z) when symbolic !*complex, df(struvel(~n,~x),x) => (x*StruveL(-1 + n,x)- n*StruveL(n,x))/x })$ algebraic (let struve!*rules); algebraic (operator lommel1, lommel2); algebraic (lommel!*rules := { lommel1(~a,~b,~z) => -(2**a)*besselj(a,z)*gamma(a+1)+z**a when numberp a and numberp b and a = b+1, lommel1(~a,~b,~z) => lommel1(a,-b,z) when numberp b and b < 0 and a neq b and a neq (b+1), lommel1(~a,~b,~z) => (sqrt(pi)*(2**a)*gamma((2*a + 1)/2)*struveh(a,z))/2 when a = b, lommel2(~a,~b,~z) => z**b when numberp a and numberp b and a = b+1, lommel2(~a,~b,~z) => lommel2(a,-b,z) when numberp b and b < 0 and a neq b and a neq (b+1), lommel2(~a,~b,~z) => (sqrt(pi)*(2**a)*gamma((2*a + 1)/2)*(-bessely(a,z)+struveh(a,z)))/2 when a = b })$ algebraic (let lommel!*rules); algebraic (operator whittakerm, whittakerw); algebraic (whittaker!*rules := { whittakerm(~k,~m,~z) => exp(-z/2)*(z**(1/2+m))*kummerm(1/2+m-k,1+2*m,z), whittakerw(~k,~m,~z) => exp(-z/2)*(z**(1/2+m))*kummeru(1/2+m-k,1+2*m,z), df(WhittakerM(~n,~m,~z),z) => 1/(2*z)* ((1+2*m-2*n)*WhittakerM(n-1,m,z) + (2*n-z)*WhittakerM(n,m,z)), df(WhittakerW(~n,~m,~z),z) => 1/(4*z)* ((1-4*m^2-4*n+4*n^2)*WhittakerW(n-1,m,z) + (4*n-2*z)*WhittakerW(n,m,z)) % AS (8.5.4) })$ algebraic (let whittaker!*rules); %Handbook of Mathematical Functions - page 496 algebraic procedure struve_compute_term(n,x,h_or_l); begin scalar dmode!*!*; lisp(dmode!*!* := dmode!*); return begin scalar pre,term,k,precis,result,!*complex,!*rounded, dmode!*,expo,!*msg; lisp (dmode!* := dmode!*!*); if h_or_l = l then << on complex; off rounded; expo := e^(-i*n*pi/2); on rounded; return (-i*expo*struveh(n,i*x))>> else << pre := precision 0; precis := 10.0^(-pre-2); result := 0; << if n > -2 then <<k:=1, term := 2^(n+2)/(pi * (for i:= 1 :n+1 product(2i-1))) ; result := term >> else for kk:=0:-(n+2) do << k:=kk+1; term := (-1)^kk*(1/2*x)^(2*kk)/ (Gamma(kk+3/2) * Gamma(kk+n+3/2)); result := result + term>>; while abs(term) > precis do << term:= term*(-0.25)*(x^2)/((k+0.5)*(k+n+0.5)); result := result + term; k := k+1>>; >>; >>; return result; end; end; symbolic operator struve_compute_term; % Lambert's W (Omega) function. % see: "On Lambert's W function" by R. Corless, G. Gonnet et. al. % only the principal branch is implemented algebraic << % Remove autoload properties. lisp null remprop('lambert_w,'simpfn); lisp null remflag('(lambert_w),'full); operator lambert_w; let { lambert_w(0) => 0, lambert_w(-1/e) => -1, sum((- ~n)^(n-1)/factorial n *~z^n,n,1,infinity) => lambert_w(z), df(lambert_w(~z),z) => 1/((1 + lambert_w(z))*e^lambert_w z), log(lambert_w(~z)) => log(z) - lambert_w z, e^(lambert_w ~z) => ~z/lambert_w z, int(lambert_w(~z),z) => z*(lambert_w z -1 +1/lambert_w z), lambert_w(~z) => num_lambert_w(z) when numberp z and lisp !*rounded}; procedure num_lambert_w(z); if z=0 then 0 else begin scalar wjnew,wj,accu,expwj,oldprec,!*complex,olddmode!*; lisp setq(olddmode!* ,dmode!*); on complex; oldprec := precision 5; accu := 10^(- lisp !:prec!:); if (abs z) <= 1 then % starting point for iteration if z >= -1/e then wj := 0 else wj := log(z) else wj := log(z) - log(log(z)); wjnew := 100; while abs(wjnew) > accu do << expwj := exp(wj); wjnew := - (wj*expwj -z)/ (expwj*(wj+1)-(1/2(wj+2)*(wj*expwj -z))/(wj+1)); wj := wj + wjnew >>; precision oldprec; accu := 10^(- lisp !:prec!:); while abs(wjnew) > accu do << expwj := exp(wj); wjnew := - (wj*expwj -z)/ (expwj*(wj+1)-(1/2(wj+2)*(wj*expwj -z))/(wj+1)); wj := wj + wjnew >>; lisp setq(dmode!*,olddmode!*); return wj; end; >>; endmodule; end;