airy ai and bi functions

[roadrunner]

Functions for calculation of airy ai and bi follow.

Code:


local c(k)
BEGIN
 LOCAL j, n;
 n:=1/((216^k)*k!);
 FOR j FROM 0 TO 2*k-1 DO
  n:=n*(2*k+2*j+1);
 END;
 RETURN n;
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, k, old1, old2, new1, new2;
 if z ≥ 0 then
  s:=2/3*z^(3/2);
  new1:=1; 
  for k from 1 to 48 do
   old1:=new1;
   new1:=new1+(−1)^k*c(k)*s^(−k);
   if new1 == old1 then break; end;
  end; 
  return 0.282094791774/z^(1/4)/exp(s)*new1;
 end; 
 z:=−z;
 s:=2/3*z^(3/2);
 new1:=1;
 new2:=c(1)/s;
 for k from 1 to 24 do
  old1:=new1;
  old2:=new2;
  new1:=new1+(−1)^k*c(2*k)*s^(−2*k);
  new2:=new2+(−1)^k*c(2*k+1)*s^(−2*k-1);
  if (new1 == old1) and (new2 == old2) then
   break;
  end;
 end;
 s:=s+π/4;
 return (0.564189583547/z^(1/4))*
 (sin(s)*new1-cos(s)*new2);
end;

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

local bibigz(z)
begin
 local s, k, old1, old2, new1, new2;
 if z ≥ 0 then
  s:=2/3*z^(3/2);
  new1:=1; 
  for k from 1 to 48 do
   old1:=new1;
   new1:=new1+c(k)/s^k;
   if new1 == old1 then break; end;
  end; 
  return 0.564189583548/z^(1/4)*exp(s)*new1;
 end;
 z:=−z;
 s:=2/3*z^(3/2);
 new1:=1;
 new2:=c(1)/s;
 for k from 1 to 24 do
  old1:=new1;
  old2:=new2;
  new1:=new1+(−1)^k*c(2*k)*s^(−2*k);
  new2:=new2+(−1)^k*c(2*k+1)*s^(−2*k-1);
  if (new1 == old1) and (new2 == old2) then
   break;
  end;
 end;
 s:=s+π/4;
 return (0.564189583548/z^(1/4))*
 (cos(s)*new1+sin(s)*new2);
end;

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

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

Road

[Gerald H]

I have adapted your programmes for the HP 39gs & everything functions very nicely

EXCEPT

aibigz for positive values of input (negatives produce good results).

I wonder whether you have an update that corrects the deviant results?

[roadrunner]

Gerald,

Sorry, I didn't mean to ignore you, I just noticed this post a couple days ago.

To your question, I don't have an update that is more accurate, but I do have a faster version.

For what it's worth, you have inspired me to start working on a CAS version that may address some of the accuracy issues.

Here is the latest version that I have:

Code:

CLIST:={5/72,385/10368,85085/2239488,37182145/644972544,5391411025/46438023168,5849680962125/20061226008576,1267709431363375/1444408272617472,2562040760785380875/831979165027663872,6653619855759634132375/539122498937926189056,4318199286388002551911375/77633639847061371224064,1556514560957130010757145625/5589622068988418728132608,7404339766473067461171741738125/4829433467605993781106573312,3201522602103470169172796946923125/347719209667631552239673278464,2998911957427493414180861368710704375/50071566192138943522512952098816,4537353791587797535655643250859295719375/10815458297502011800862797653344256,39225423528276509695743035903678611493996875​/12459407958722317594593942896652582912,22605380841560292087599677808725844518040​434375/897077373028006866810763888558985969664,2491339022548359790974360491299675324333​23627246875/1162612275444296899386749999572445816684544,161504330214358671291743148270200531​420065637726303125/83708083831989376755845999969216098801287168,44203735179669968332550099681553885​4496719650456891653125/24107928143612940505683647991134236454770704384,95511642084644038718545751097643​2167751840673308640893359375/5207312479020395149227667966084995074230472146944,144457017214747533832767787364​6809270444306661985078049347265625/749852996978936901488784187116239290689187989159936,1144413613334745571289939936​335587465946335290783917701615501171875/53989415782483456907192461472369228929621535219515392,22747509392254737720530136​114542472060615306574911932155011316793359375/93293710472131413535628573424254027590386012859322597376,19658397616786544338082​143630187604354783747942038891768360779972821171875/6717147153993461774565257286546289986507792925871227011072,353926766323760823256​00967051118529224901047737186173987575696558759963671875/967269190175058495537397049262665758057122181325456689594368,8946206872365702329​44215644151123063217823783652854919883950881915775601733984375/1880371305700313715324699863766622233663045520496687804571451392,347892424675238​2615853499149919674414836054459259194824880146708089903950685826953125/541546936041690350013513560764787203294957109903046087716578000896,3507355455251​703811714094987800052443674542214667623417069684458773259372627638711328125/38991379395001705200972976375064678637236911913019318315593616064512,21973581927​15192438038880509856732855962100697489266070794157313421446996951215652647070312​5/16844275898640736646820325794027941171286345946424345512336442139869184,23735012​60744136089926835931372667731673901350175099173565882504513456538481014713835069​3359375/1212787864702133038571063457170011764332616908142552876888223834070581248,847743​44529998308723916798960837573372196734524204017182252625413707127184926402534047​1714716796875/2794263240273714520867730205319707105022349356360441828350467713698619195392,293​01866848864233599927641939000048883493564030406917611694791553223260751800788643​15461369567939453125/603560859899122336507429724349056734684827460973855434923701026158901746204672,6​96298538290334379256633500829285279262640885962557794564490879044447096817939093​4220687820378591181640625/86912763825473616457069880306264169794615154380235182629012947766881851453472768​,8523688820871250422642989240865922225716528216875939630890174775045524532332742​987905296276117735405068359375/62577189954341003849090313820510202252122911153769331492889322392154933046500392​96,38667714335882427542319920691188256176963030255857700135533277866994022040927​4885646323765566081066650926123046875/16220007636165188197684209342276244423750258571057010722956912364046558645652901​855232,5011858314610577236883882261046797841834097410757210068648025423425492743​28853937592993211493312154415634164404296875/11678405498038935502332630726438895985100186171161047720528976902113522224870089​33576704,13352909460205370016226263365076313344814939212102143730791141839959045​69252679746047100492579102872035464054542626953125/16816903917176067123358988246072010218544268086471908717561726739043472003812928​6435045376,548054761590906098427526718746443731645608239646256986956786888580534​4623201171785371933921734091687938173504479303564453125/36324512461100304986455414611515542072055619066779322829933329756333899528235925​869969801216,6156299136950648203636407631678802437575117355946404734485587119425​1436152418762665082933742839051930609502975816016939501953125/20922919177593775672198318816232952233504036582464889950041597939648326128263893​301102605500416,8865821525396366419114919234450607071381044006874523598724426621​1331123112672337599020044459656873688849709834318456590070556640625/15064501807867518483982789547687725608122906339374720764029950516546794812350003​17679387596029952,78522048158571218412266806825462040943487226470599977078588325​3186270805659478130589492242337906742583945065799748764744943471435546875/65078647809987679850805650846010974627090955386098793700609386231482153589352013​7237495441484939264,118734467286475795171720931316241101492708905868485607200625​3837335453815227556454164171513698907774640713770310824754243743475725830078125/46856626423191129492580068609127901731505487877991131464438758086667150584333449​881099671786915627008,7352577916502610672465558798607631918526229582766940096077​27205797991886653231481711826099987531826229104544475113179862009766691057885742​1875/13494708409879045293863059759428835698673580508861445861758362328960139368288033​565756705474631700578304,1048257033555777203573414717917490082624284551615082649​49773667730619703280151212347655047075222412465483434905816886052926732437144122​770263671875/87445710496016213504232627241098855327404801697422169184194187891661703106506457​50610345147561341974740992,33965806706845346128829883414261281916336872264832188​89296470690271492776936203956612605384035195777821805646241740623084506406468767​71715386962890625/12592182311426334744609498322718235167146291444428792362523963056399285247336929​88087889701248833244362702848,56250989290228110779355992240762880862168364393143​46993322752959971535814433687390869263091003139016134380555032305635724015428942​75474042778192138671875/90663712642269610161188387923571293203453298399887305010172534006074853780825895​142328058489915993594114605056,4568761601141617385610073045787001946506176724375​50577144667318161848110384118523573792417514365954029450523060278896039140257154​15948277228487543695068359375/31333379089168377271706706866386238931113459927001052611515627752499469466653429​3611885770141149673861260075073536,789528624693609664002661541424706048621027605​08151980042610192857124107026522457593179633750742787529289350900357624652869958​275597474536140339586478924560546875/22560032944201231635628828943798092030401691147440757880291251981799618015990469​140055775450162776518010725405294592,2785614893643993616534190450454647880744709​59624861815986337282438505274410976534880256383799370702960838687846641771300255​786787963009658410346129014941634521484375/32486447439649773555305513679069252523778435252314691347619402853791449943026275​56168031664823439818593544458362421248,15043139724175696115642432606116996981868​70967550614059801511424368645689252923575569525724394189846195282102221232224365​881323886408476570021283904365689221014404296875/70170726469643510879459909546789585451361420144999733310857910164189531876936755​2132294839601863000816205603006282989568,110482603957317149079476349071864040064​45221489909998352897608174188091113292414171127046737537235688091303790290506852​752403935438663239805840163333624839211907098388671875/20209169223257331133284453949475400609992089001759923193527078127286585180557785​5014100913805336544235067213665809500995584,206833857495263487713485336438857017​34405685310384149935336860012285218648150691065535786874437397400152967044838007​178055590016569049156920288053700482833423483664324951171875/14550601840745278415964806843622288439194304081267144699339496251646341330001605​561015265793984231184924839383938284071682048,2131154017473946398353438861065051​04960795859732605165688730404508583207384950275531961087218140611591956126539897​3745605313828537225117981595720189136649707455486321049993896484375/56572739956817642481271169008003457451587454267966658591031961426400975091046242​421227353407010690846987775524752048470699802624,4143157151243667611708026367625​10151780601769125614678932136336401459262720651058386453440923630634444907433277​7913959069966929388050831639674952385880659453994143179604919952392578125/40732372768908702586515241685762489365142967072935994185543012227008702065553294​54328369445304769740983119837782147489890385788928,65630568673664869017749214824​87276882884575309927683381912848237410830049539913229883127113830968985774794391​2868427050210468995027716638052193870758339560508020080924355935560150146484375/23461846714891412689832779210999193874322349034011132650872775042757012389758697​65693140800495547370806277026562516954176862214422528,39695095054060767120593086​94752581039739193373110313637890940952097580825278832677044627992316386647043589​72138689527964686129763242111205002313844522953014468428612609695318250830230712​890625/50677588904165451410038803095758258768536273913464046525885194092355146761878786​9389718412907038232094155837737503662102202238315266048,162915514083771260808969​30350280929238993052912446987568807506312400212669155584277742367871499608797164​65835273067803720169236225891010956226910209964358301208511789704087809432140168​609619140625/72975728021998250030455876457891892626692234435388226997274679492991411337105453​192119451458613505421558440634200527342717122317398310912,3401813998064130172522​47309892128752450207989771988394493962433079452372915586528180512965197480899149​88040933582395629646388033105270318097336721548640516887116106263758390823493452​94105072021484375/52542524175838740021928231049682162691218408793479523438037769234953816162715926​29832600505020172390352207725662437968675632806852678385664,86708836996656613967​42531681840469771203351451298212187256608456762161533245385016793094969918590638​43130175356081682204056784575820235137983015695553298134935702432556937623700024​62201444180755615234375/45396740887924671378945991626925388565212705197566308250464632619000097164586560​32175366836337428945264307474972346404935746745120714125213696};

f(z)
BEGIN
 A:=1;
 B:=1;
 C:=1;
 D:=1;
 FOR K FROM 3 TO 105 STEP 3 DO
  E:=D;
  A:=A*(K-2);
  B:=B*z^3;
  C:=C*(K-2)*(K-1)*K;
  D:=D+A*B/C;
  IF D==E THEN RETURN D;END;
 END;
 RETURN D;
END;

g(z)
BEGIN
 A:=1;
 B:=z;
 C:=1;
 D:=z;
 FOR K FROM 3 TO 105 STEP 3 DO
  E:=D;
  A:=A*(K-1);
  B:=B*z^3;
  C:=C*(K-1)*(K)*(K+1);
  D:=D+A*B/C;
  IF D==E THEN RETURN D;END;
 END;
 RETURN D;
END;

aismallz(z)
BEGIN
 RETURN 0.355028053887817*f(z)-
  0.258819403792807*g(z);
END;

aibigz(z)
BEGIN
 IF z ≥ 0 THEN
  S:=2/3*z^(3/2);
  A:=1;
  FOR K FROM 1 TO 48 DO
   B:=A;
   A:=A+(−1)^K*CLIST(K)*S^(−K);
   IF A == B THEN BREAK;END;
  END;
  RETURN 0.282094791774/z^(0.25)/exp(S)*A;
 END;
 z:=−z;
 S:=2/3*z^(3/2);
 A:=1;
 C:=CLIST(1)/S;
 FOR K FROM 1 TO 24 DO
  B:=A;
  D:=C;
  A:=A+(−1)^K*CLIST(2*K)*S^(−2*K);
  C:=C+(−1)^K*CLIST(2*K+1)*S^(−2*K-1);
  IF (A == B) AND (C == D) THEN BREAK;END;
 END;
 S:=S+π/4;
 RETURN (0.564189583547/z^(0.25))*
 (sin(S)*A-cos(S)*C);
END;

bismallz(z)
BEGIN
 RETURN (0.355028053887817*f(z)+
  0.258819403792807*g(z))*√3;
END;

bibigz(z)
BEGIN
 if z ≥ 0 then
  S:=2/3*z^(3/2);
  A:=1;
  FOR K FROM 1 TO 48 DO
   B:=A;
   A:=A+CLIST(K)/S^K;
   IF A == B THEN BREAK;END;
  END;
  RETURN 0.564189583548/z^(0.25)*exp(S)*A;
 END;
 z:=−z;
 S:=2/3*z^(3/2);
 A:=1;
 C:=CLIST(1)/S;
 FOR K FROM 1 TO 24 DO
  B:=A;
  D:=C;
  A:=A+(−1)^K*CLIST(2*K)*S^(−2*K);
  C:=C+(−1)^K*CLIST(2*K+1)*S^(−2*K-1);
  IF (A == B) AND (C == D) THEN BREAK;END;
 END;
 S:=S+π/4;
 RETURN (0.564189583548/z^(0.25))*
 (cos(S)*A+sin(S)*C);
END;

EXPORT Airy_Ai(z)
BEGIN
 HAngle:=0;
 IF (z < 6.5)  AND (z > −6.5) THEN
  RETURN aismallz(z);
 END;
 RETURN aibigz(z);
END;

EXPORT Airy_Bi(z)
BEGIN
 HAngle:=0;
 IF ABS(z) < 6.5 THEN RETURN bismallz(z);END;
 RETURN bibigz(z);
END;

It's good for making cool plots like this:

   

but not much else

-road

[Gerald H]

You may like some of the modifications I incorporated in your programmes for the 38G, including my attempt to improve accuracy.

Code:


AIRYAI

Ans►Z:
IF ABS(Ans)<7
THEN
RUN AIRYF:
Ans*.355028053888►A:
RUN AIRYG:
A-Ans*.258819403793: 
ELSE
IF Z>=0
THEN
(2/3)*Z^1.5►S:
1►Q:
FOR K=1 TO 48 STEP 1;
 Q►O:
K:
RUN AIRYC:
Q+IFTE(K MOD 2,-1,1)*Ans*S^(-K)►Q:
IF O==Q
THEN
BREAK:
END:
END:
.282094791774/Z^.25/e^S*Q:
ELSE
-Z►Z:
(2/3)*Z^1.5►S:
1►Q:
RUN AIRYC:
Ans/S►R:
FOR K=1 TO 24 STEP 1;
 Q►O:
R►P:
IFTE(K MOD 2,-1,1)►L:
2*K:
RUN AIRYC:
Q+L*Ans*S^(-2*K)►Q:
2*K+1:
RUN AIRYC:
R+L*Ans*S^(-2*K-1)►R:
IF O==Q AND P==R
THEN
BREAK:
END:
END:
S+π/4:
(.564189583547/Z^.25)*(SIN(Ans)*Q-COS(Ans)*R):
END: 
END:
MSGBOX "Ai:
"Ans: 

AIRYBI

Ans►Z:
IF ABS(Ans)<7 OR Z>=0
THEN
RUN AIRYF:
Ans*.355028053888►A:
RUN AIRYG:
(A+Ans*.258819403793)*√3:
ELSE
-Z►Z:
(2/3)*Z^1.5►S:
1►Q:
RUN AIRYC:
Ans/S►R:
FOR K=1 TO 24 STEP 1;
 Q►O:
R►P:
IFTE(K MOD 2,-1,1)►L:
2*K:
RUN AIRYC:
Q+L*Ans*S^(-2*K)►Q:
2*K+1:
RUN AIRYC:
R+L*Ans*S^(-2*K-1)►R:
IF O==Q AND P==R
THEN
BREAK:
END:
END:
S+π/4:
(.564189583548/Z^.25)*(SIN(Ans)*R+COS(Ans)*Q):
END:
MSGBOX "Bi:
"Ans: 

AIRYC

Ans►M:
1/(M!*216^M):
FOR J=0 TO 2*M-1 STEP 1;
Ans*(2*M+2*J+1):
END: 

AIRYF

1►T:
0►U:
1►Q:
FOR K=1 TO 80 STEP 1;
 3*K►P:
Q*(P-2)►Q:
T+Q*Z^P/P!►T:
IF T==U
THEN
BREAK:
END:
T►U:
END:

AIRYG

Z►T:
0►U:
1►Q:
FOR K=1 TO 80 STEP 1;
3*K+1►P:
Q*(P-2)►Q:
T+Q*Z^P/P!►T:
IF T==U
THEN
BREAK:
END:
T►U:
END: