## Monday, July 23, 2018

### Solution of the differential equation

```''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Cephes Math Library Release 2.8:  June, 2000
'
'Contributors:
'    * Sergey Bochkanov (ALGLIB project). Translation from C to
'      pseudocode.
'
'
'This program is free software; you can redistribute it and/or modify
'the Free Software Foundation (www.fsf.org); either version 2 of the
'
'This program is distributed in the hope that it will be useful,
'but WITHOUT ANY WARRANTY; without even the implied warranty of
'MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
'GNU General Public License for more details.
'
'A copy of the GNU General Public License is available at
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Routines
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Airy function
'
'Solution of the differential equation
'
'y"(x) = xy.
'
'The function returns the two independent solutions Ai, Bi
'and their first derivatives Ai'(x), Bi'(x).
'
'Evaluation is by power series summation for small x,
'by rational minimax approximations for large x.
'
'
'
'ACCURACY:
'Error criterion is absolute when function <= 1, relative
'when function > 1, except * denotes relative error criterion.
'For large negative x, the absolute error increases as x^1.5.
'For large positive x, the relative error increases as x^1.5.
'
'Arithmetic  domain   function  # trials      peak         rms
'IEEE        -10, 0     Ai        10000       1.6e-15     2.7e-16
'IEEE          0, 10    Ai        10000       2.3e-14*    1.8e-15*
'IEEE        -10, 0     Ai'       10000       4.6e-15     7.6e-16
'IEEE          0, 10    Ai'       10000       1.8e-14*    1.5e-15*
'IEEE        -10, 10    Bi        30000       4.2e-15     5.3e-16
'IEEE        -10, 10    Bi'       30000       4.9e-15     7.3e-16
'
'Cephes Math Library Release 2.8:  June, 2000
'Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Sub Airy(ByVal x As Double, _
ByRef Ai As Double, _
ByRef Aip As Double, _
ByRef Bi As Double, _
ByRef Bip As Double)
Dim z As Double
Dim zz As Double
Dim t As Double
Dim f As Double
Dim g As Double
Dim uf As Double
Dim ug As Double
Dim k As Double
Dim zeta As Double
Dim theta As Double
Dim domflg As Long
Dim c1 As Double
Dim c2 As Double
Dim sqrt3 As Double
Dim sqpii As Double
Dim AFN As Double
Dim AFD As Double
Dim AGN As Double
Dim AGD As Double
Dim APFN As Double
Dim APFD As Double
Dim APGN As Double
Dim APGD As Double
Dim AN As Double
Dim APN As Double
Dim APD As Double
Dim BN16 As Double
Dim BD16 As Double
Dim BPPN As Double
Dim BPPD As Double

sqpii = 0.564189583547756
c1 = 0.355028053887817
c2 = 0.258819403792807
sqrt3 = 1.73205080756888
domflg = 0#
If x > 25.77 Then
Ai = 0#
Aip = 0#
Bi = MaxRealNumber
Bip = MaxRealNumber
Exit Sub
End If
If x < -2.09 Then
domflg = 15#
t = Sqr(-x)
zeta = -(2# * x * t / 3#)
t = Sqr(t)
k = sqpii / t
z = 1# / zeta
zz = z * z
AFN = -0.131696323418332
AFN = AFN * zz - 0.626456544431912
AFN = AFN * zz - 0.693158036036933
AFN = AFN * zz - 0.279779981545119
AFN = AFN * zz - 0.04919001326095
AFN = AFN * zz - 4.06265923594885E-03
AFN = AFN * zz - 1.59276496239262E-04
AFN = AFN * zz - 2.77649108155233E-06
AFN = AFN * zz - 1.67787698489115E-08
AFD = 1#
AFD = AFD * zz + 13.3560420706553
AFD = AFD * zz + 32.6825032795225
AFD = AFD * zz + 26.73670409415
AFD = AFD * zz + 9.1870740290726
AFD = AFD * zz + 1.47529146771666
AFD = AFD * zz + 0.115687173795188
AFD = AFD * zz + 4.40291641615211E-03
AFD = AFD * zz + 7.54720348287414E-05
AFD = AFD * zz + 4.5185009297058E-07
uf = 1# + zz * AFN / AFD
AGN = 1.97339932091686E-02
AGN = AGN * zz + 0.391103029615688
AGN = AGN * zz + 1.06579897599596
AGN = AGN * zz + 0.93916922981665
AGN = AGN * zz + 0.351465656105548
AGN = AGN * zz + 6.33888919628925E-02
AGN = AGN * zz + 5.85804113048388E-03
AGN = AGN * zz + 2.82851600836737E-04
AGN = AGN * zz + 6.98793669997261E-06
AGN = AGN * zz + 8.11789239554389E-08
AGN = AGN * zz + 3.41551784765924E-10
AGD = 1#
AGD = AGD * zz + 9.30892908077442
AGD = AGD * zz + 19.8352928718312
AGD = AGD * zz + 15.5646628932865
AGD = AGD * zz + 5.47686069422975
AGD = AGD * zz + 0.954293611618962
AGD = AGD * zz + 8.64580826352392E-02
AGD = AGD * zz + 4.12656523824223E-03
AGD = AGD * zz + 1.01259085116509E-04
AGD = AGD * zz + 1.17166733214414E-06
AGD = AGD * zz + 4.9183457006293E-09
ug = z * AGN / AGD
theta = zeta + 0.25 * PI()
f = Sin(theta)
g = Cos(theta)
Ai = k * (f * uf - g * ug)
Bi = k * (g * uf + f * ug)
APFN = 0.185365624022536
APFN = APFN * zz + 0.886712188052584
APFN = APFN * zz + 0.987391981747399
APFN = APFN * zz + 0.401241082318004
APFN = APFN * zz + 7.10304926289631E-02
APFN = APFN * zz + 5.90618657995662E-03
APFN = APFN * zz + 2.33051409401777E-04
APFN = APFN * zz + 4.08718778289035E-06
APFN = APFN * zz + 2.48379932900442E-08
APFD = 1#
APFD = APFD * zz + 14.7345854687503
APFD = APFD * zz + 37.542393343549
APFD = APFD * zz + 31.4657751203046
APFD = APFD * zz + 10.9969125207299
APFD = APFD * zz + 1.78885054766999
APFD = APFD * zz + 0.141733275753663
APFD = APFD * zz + 5.44066067017226E-03
APFD = APFD * zz + 9.39421290654511E-05
APFD = APFD * zz + 5.65978713036027E-07
uf = 1# + zz * APFN / APFD
APGN = -3.55615429033082E-02
APGN = APGN * zz - 0.637311518129436
APGN = APGN * zz - 1.70856738884312
APGN = APGN * zz - 1.50221872117317
APGN = APGN * zz - 0.563606665822103
APGN = APGN * zz - 0.102101031120217
APGN = APGN * zz - 9.48396695961445E-03
APGN = APGN * zz - 4.60325307486781E-04
APGN = APGN * zz - 1.14300836484517E-05
APGN = APGN * zz - 1.33415518685547E-07
APGN = APGN * zz - 5.63803833958894E-10
APGD = 1#
APGD = APGD * zz + 9.8586580169613
APGD = APGD * zz + 21.6401867356586
APGD = APGD * zz + 17.3130776389749
APGD = APGD * zz + 6.17872175280829
APGD = APGD * zz + 1.08848694396321
APGD = APGD * zz + 9.95005543440888E-02
APGD = APGD * zz + 4.78468199683887E-03
APGD = APGD * zz + 1.18159633322839E-04
APGD = APGD * zz + 1.37480673554219E-06
APGD = APGD * zz + 5.79912514929148E-09
ug = z * APGN / APGD
k = sqpii * t
Aip = -(k * (g * uf + f * ug))
Bip = k * (f * uf - g * ug)
Exit Sub
End If
If x >= 2.09 Then
domflg = 5#
t = Sqr(x)
zeta = 2# * x * t / 3#
g = Exp(zeta)
t = Sqr(t)
k = 2# * t * g
z = 1# / zeta
AN = 0.346538101525629
AN = AN * z + 12.0075952739646
AN = AN * z + 76.2796053615235
AN = AN * z + 168.089224934631
AN = AN * z + 159.756391350164
AN = AN * z + 70.5360906840444
AN = AN * z + 14.026469116339
AN = AN * z + 1#
Ai = sqpii * f / k
k = -(0.5 * sqpii * t / g)
APN = 0.613759184814036
APN = APN * z + 14.7454670787755
APN = APN * z + 82.0584123476061
APN = APN * z + 171.184781360976
APN = APN * z + 159.317847137142
APN = APN * z + 69.9778599330103
APN = APN * z + 13.9470856980482
APN = APN * z + 1#
APD = 0.334203677749737
APD = APD * z + 11.1810297306158
APD = APD * z + 71.172735214786
APD = APD * z + 158.778084372838
APD = APD * z + 153.206427475809
APD = APD * z + 68.675230459278
APD = APD * z + 13.8498634758259
APD = APD * z + 1#
f = APN / APD
Aip = f * k
If x > 8.3203353 Then
BN16 = -0.253240795869364
BN16 = BN16 * z + 0.575285167332467
BN16 = BN16 * z - 0.329907036873225
BN16 = BN16 * z + 0.06444040689482
BN16 = BN16 * z - 3.82519546641337E-03
BD16 = 1#
BD16 = BD16 * z - 7.15685095054035
BD16 = BD16 * z + 10.6039580715665
BD16 = BD16 * z - 5.23246636471251
BD16 = BD16 * z + 0.957395864378384
BD16 = BD16 * z - 0.055082814716355
f = z * BN16 / BD16
k = sqpii * g
Bi = k * (1# + f) / t
BPPN = 0.465461162774652
BPPN = BPPN * z - 1.08992173800494
BPPN = BPPN * z + 0.638800117371828
BPPN = BPPN * z - 0.126844349553103
BPPN = BPPN * z + 7.6248784434211E-03
BPPD = 1#
BPPD = BPPD * z - 8.70622787633159
BPPD = BPPD * z + 13.8993162704553
BPPD = BPPD * z - 7.14116144616431
BPPD = BPPD * z + 1.34008595960681
BPPD = BPPD * z - 7.84273211323342E-02
f = z * BPPN / BPPD
Bip = k * t * (1# + f)
Exit Sub
End If
End If
f = 1#
g = x
t = 1#
uf = 1#
ug = x
k = 1#
z = x * x * x
Do While t > MachineEpsilon
uf = uf * z
k = k + 1#
uf = uf / k
ug = ug * z
k = k + 1#
ug = ug / k
uf = uf / k
f = f + uf
k = k + 1#
ug = ug / k
g = g + ug
t = Abs(uf / f)
Loop
uf = c1 * f
ug = c2 * g
If domflg Mod 2# = 0# Then
Ai = uf - ug
End If
If domflg \ 2# Mod 2# = 0# Then
Bi = sqrt3 * (uf + ug)
End If
k = 4#
uf = x * x / 2#
ug = z / 3#
f = uf
g = 1# + ug
uf = uf / 3#
t = 1#
Do While t > MachineEpsilon
uf = uf * z
ug = ug / k
k = k + 1#
ug = ug * z
uf = uf / k
f = f + uf
k = k + 1#
ug = ug / k
uf = uf / k
g = g + ug
k = k + 1#
t = Abs(ug / g)
Loop
uf = c1 * f
ug = c2 * g
If domflg \ 4# Mod 2# = 0# Then
Aip = uf - ug
End If
If domflg \ 8# Mod 2# = 0# Then
Bip = sqrt3 * (uf + ug)
End If
End Sub
```