Beta and Airy Functions

[douganc]

Don't seem to be fully implemented on HP Prime.

XCAS:

Input: Beta(x,y)

Output: Gamma x * Gamma y / Gamma (x+y)

Prime: ??????

XCAS:

Input: Beta(n,2)

Output: Gamma n / Gamma (n+2)

Prime: ??????

XCAS:

Input: Airy_Ai(1)

Output: 0.135292416313

Prime: ????????

XCAS:

Input: Airy_Bi(1)

Output: 1.20742359495

Prime: ?????????

Latest Software Rev.

[Tim Wessman]

All of these seem to be working fine for me and match xcas output.

Airy cant' be implemented on our end due to it using a numerical library we can't use at the moment - hence why it wasn't included.

[Gerald H]

Anyone know an efficient way of caculating the Airy function?

[DrD]

A burrito rich diet?

[parisse]

http://people.math.sfu.ca/~cbm/aands/page_446.htm

[roadrunner]

(06-24-2015 06:48 PM)parisse Wrote:  http://people.math.sfu.ca/~cbm/aands/page_446.htm

That method becomes unstable outside a range for z of about -10 to 10. I believe that is because the calculator is trying to divide one huge number by another huge number to get a small number.

   

How can one calculate the airy function outside that range?

Code:


local f(z)
begin
 local t,t2,k,p,p2;
 t:=1;
 p2:=1;
 for k from 1 to 80 do
  p:=3*k;
  p2:=p2*(3*k-2);
  t:=t+p2*z^(p)/(p)!;
  if t == t2 then return t; end;
  t2:=t;
 end;
 return t;
end;

local g(z)
begin
 local t,t2,k,p,p2;
 t:=z;
 t2:=0;
 p2:=1;
 for k from 1 to 80 do
  p:=3*k+1;
  p2:=p2*(3*k-1);
  t:=t+p2*z^(p)/(p)!;
  if t == t2 then return t; end;
  t2:=t;
 end;
 return t;
end;

export ai(z)
begin
 0.355028053887817*f(z)-
 0.258819403792807*g(z);
end;

export bi(z)
begin
 (0.355028053887817*f(z)+
 0.258819403792807*g(z))*√3;
end;

[Gerald H]

Thank you, parisse & more so roadrunner, here my translation of your programme for the 48G.

Code:



Ai

« DUP f .355028053888
* SWAP g
.258819403793 * -
»

Bi

« DUP f .355028053888
* SWAP g
.258819403793 * + 3.
√ *
»

g

« DUP 0. DUP 1.  Z T
T2 P P2
  « 1. 80.
    FOR K K 3. * 1. +
'P' STO P2 3. K * 1.
- * 'P2' STO T P2 Z P
^ * P ! / + 'T' STO T
T2 ==
      IF
      THEN 81. 'K'
STO
      ELSE T 'T2' STO
      END
    NEXT T
  »
»

f

« 1. 4. NDUPN DROP 
Z T T2 P P2
  « 1. 80.
    FOR K K 3. * 'P'
STO P2 3. K * 2. - *
'P2' STO T P2 Z P ^ *
P ! / + 'T' STO T T2
==
      IF
      THEN 81. 'K'
STO
      ELSE T 'T2' STO
      END
    NEXT T
  »
»

[parisse]

For |z| large, you can use asymptotic series expansion near infinity.
http://people.math.sfu.ca/~cbm/aands/page_448.htm

[Gerald H]

& here's a sys version for the 49G.

Code:



Ai

::
  CK1&Dispatch
  BINT1
  ::
    DUP
    ID f
    % .355028053888
    %*
    SWAP
    ID g
    % .258819403793
    %*
    %-
  ;
;

Bi

::
  CK1&Dispatch
  BINT1
  ::
    DUP
    ID f
    % .355028053888
    %*
    SWAP
    ID g
    % .258819403793
    %*
    %+
    %3
    %SQRT
    %*
  ;
;

f

::
  %1
  BINT4
  NDUPN
  DROP
  '
  NULLLAM
  BINT5
  NDUPN
  DOBIND
  BINT81
  ONE_DO
  INDEX@
  UNCOERCE
  DUP
  %3
  %*
  2PUTLAM
  1GETLAMSWAP_
  %3
  %*
  %2
  %-
  %*
  1PUTLAM
  4GETLAM
  1GETLAM
  5GETLAM
  2GETLAM
  %^
  %*
  2GETLAM
  %NFACT
  %/
  %+
  DUP4PUTLAM
  3GETLAM
  %=
  ITE
  ExitAtLOOP
  ::
    4GETLAM
    3PUTLAM
  ;
  LOOP
  4GETLAM
  ABND
;

g

::
  DUP
  %0
  DUP
  %1
  '
  NULLLAM
  BINT5
  NDUPN
  DOBIND
  BINT81
  ONE_DO
  INDEX@
  UNCOERCE
  DUP
  %3
  %*
  %1+
  2PUTLAM
  1GETLAMSWAP_
  %3
  %*
  %1-
  %*
  1PUTLAM
  4GETLAM
  1GETLAM
  5GETLAM
  2GETLAM
  %^
  %*
  2GETLAM
  %NFACT
  %/
  %+
  DUP4PUTLAM
  3GETLAM
  %=
  ITE
  ExitAtLOOP
  ::
    4GETLAM
    3PUTLAM
  ;
  LOOP
  4GETLAM
  ABND
;

[roadrunner]

(06-27-2015 12:41 PM)parisse Wrote:  For |z| large, you can use asymptotic series expansion near infinity.
http://people.math.sfu.ca/~cbm/aands/page_448.htm

That's much better.

   

Now it matches this web site: http://keisan.casio.com/exec/system/1180573401 pretty close. I only checked out to +/-100.

Code:


export c(k)
begin
 if k == 0 then return 1; end;
 return Gamma(3*k+1/2)/54^k/k!/Gamma(k+1/2);
end;

export d(k)
begin
 return c(k)*(6*k+1)/(1-6*k); 
end;

local f(z)
begin
 local t,t2,k,p,p2;
 t:=1;
 t2:=0;
 p2:=1;
 for k from 1 to 80 do
  p:=3*k;
  p2:=p2*(3*k-2);
  t:=t+p2*z^(p)/(p)!;
  if t == t2 then return t; end;
  t2:=t;
 end;
 return t;
end;

local g(z)
begin
 local t,t2,k,p,p2;
 t:=z;
 t2:=0;
 p2:=1;
 for k from 1 to 80 do
  p:=3*k+1;
  p2:=p2*(3*k-1);
  t:=t+p2*z^(p)/(p)!;
  if t == t2 then return t; end;
  t2:=t;
 end;
 return t;
end;

local aismallz(z)
begin
 return
 0.355028053887817*f(z)-
 0.258819403792807*g(z);
end;

local aibigz(z)
begin
 local s;
 if z ≥ 0 then
  s:=2/3*z^(3/2);
  return π^(−1/2)/2*z^(−1/4)*e^(−s)*
  sum((−1)^K*c(K)*s^(−K),K,0,20);
 end;
 z:=ABS(z);
 s:=2/3*z^(3/2);
 return π^(−1/2)*z^(−1/4)*(SIN(s+π/4)*
  sum((−1)^K*c(2*K)*s^(−2*K),K,0,20)-
  COS(s+π/4)*
  sum((−1)^K*c(2*K+1)*s^(−2*K-1),K,0,20));
end;

local bismallz(z)
begin
  return (0.355028053887817*f(z)+
  0.258819403792807*g(z))*√3;
end;

local bibigz(z)
begin
 local s;
 if z ≥ 0 then
  s:=2/3*z^(3/2);
  return π^(−1/2)*z^(−1/4)*e^(s)*
   sum(c(K)*s^(−K),K,0,20);
 end;
 z:=ABS(z);
 s:=2/3*z^(3/2);
 return π^(−1/2)*z^(−1/4)*(COS(s+π/4)*
  sum((−1)^K*c(2*K)*s^(−2*K),K,0,20)+
  SIN(s+π/4)*
  sum((−1)^K*c(2*K+1)*s^(−2*K-1),K,0,20));
end;

export ai(z)
begin
 HAngle:=0;
 if ABS(z) < 10 then
  return aismallz(z);
 end;
 return aibigz(z);
end;

export bi(z)
begin
 HAngle:=0;
 if ABS(z) < 10 then
  return bismallz(z);
 end;
 return bibigz(z);
end;

[douganc]

Thank You !!
A very useful program.
Request for you: Could you please put the program in the HP Prime Software library section.
So that other people can use the functions as well.
Thank You
Colm

[roadrunner]

Colm,

I had some issues with the version that I posted in this thread so I posted a newer version that clears up those issues and runs faster to boot.

Road

PS. Thanks to Mr. Parisse for teaching me asymptotic series expansion near infinity.