1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
|
/* Copyright 2019 David Grote
*
* This file is part of WarpX.
*
* License: BSD-3-Clause-LBNL
*/
/* -------------------------------------------------------------------------
! program to calculate the first zeroes (root abscissas) of the first
! kind bessel function of integer order n using the subroutine rootj.
! --------------------------------------------------------------------------
! sample run:
!
! (calculate the first 10 zeroes of 1st kind bessel function of order 2).
!
! zeroes of bessel function of order: 2
!
! number of calculated zeroes: 10
!
! table of root abcissas (5 items per line)
! 5.135622 8.417244 11.619841 14.795952 17.959819
21.116997 24.270112 27.420574 30.569204 33.716520
!
! table of error codes (5 items per line)
! 0 0 0 0 0
! 0 0 0 0 0
!
! --------------------------------------------------------------------------
! reference: from numath library by tuan dang trong in fortran 77
! [bibli 18].
!
! c++ release 1.0 by j-p moreau, paris
! (www.jpmoreau.fr)
! ------------------------------------------------------------------------ */
#include <AMReX_REAL.H>
#include <cmath>
void SecantRootFinder(int n, int nitmx, amrex::Real tol, amrex::Real *zeroj, int *ier);
/*----------------------------------------------------------------------
* calculate the first nk zeroes of bessel function j(n, x)
* including the trivial root (when n > 0)
*
* inputs:
* n order of function j (integer >= 0) i*4
* nk number of first zeroes (integer > 0) i*4
* outputs:
* roots(nk) table of first zeroes (abcissas) r*8
* ier(nk) table of error codes (must be zeroes) i*4
*
* reference :
* abramowitz m. & stegun irene a.
* handbook of mathematical functions
*/
void GetBesselRoots(int n, int nk, amrex::Vector<amrex::Real>& roots, amrex::Vector<int> &ier) {
using namespace amrex::literals;
amrex::Real zeroj;
int ierror, ik, k;
const amrex::Real tol = 1e-14_rt;
const amrex::Real nitmx = 10;
const amrex::Real c1 = 1.8557571_rt;
const amrex::Real c2 = 1.033150_rt;
const amrex::Real c3 = 0.00397_rt;
const amrex::Real c4 = 0.0908_rt;
const amrex::Real c5 = 0.043_rt;
const amrex::Real t0 = 4.0_rt*n*n;
const amrex::Real t1 = t0 - 1.0_rt;
const amrex::Real t3 = 4.0_rt*t1*(7.0_rt*t0 - 31.0_rt);
const amrex::Real t5 = 32.0_rt*t1*((83.0_rt*t0 - 982.0_rt)*t0 + 3779.0_rt);
const amrex::Real t7 = 64.0_rt*t1*(((6949.0_rt*t0 - 153855.0_rt)*t0 + 1585743.0_rt)*t0 - 6277237.0_rt);
roots.resize(nk);
ier.resize(nk);
// first zero
if (n == 0) {
zeroj = c1 + c2 - c3 - c4 + c5;
SecantRootFinder(n, nitmx, tol, &zeroj, &ierror);
ier[0] = ierror;
roots[0] = zeroj;
ik = 1;
}
else {
// Include the trivial root
ier[0] = 0;
roots[0] = 0.;
const amrex::Real f1 = std::pow(n, (1.0_rt/3.0_rt));
const amrex::Real f2 = f1*f1*n;
const amrex::Real f3 = f1*n*n;
zeroj = n + c1*f1 + (c2/f1) - (c3/n) - (c4/f2) + (c5/f3);
SecantRootFinder(n, nitmx, tol, &zeroj, &ierror);
ier[1] = ierror;
roots[1] = zeroj;
ik = 2;
}
// other zeroes
// k counts the nontrivial roots
// ik counts all roots
k = 2;
while (ik < nk) {
// mac mahon's series for k >> n
const amrex::Real b0 = (k + 0.5_rt*n - 0.25_rt)*MathConst::pi;
const amrex::Real b1 = 8.0_rt*b0;
const amrex::Real b2 = b1*b1;
const amrex::Real b3 = 3.0_rt*b1*b2;
const amrex::Real b5 = 5.0_rt*b3*b2;
const amrex::Real b7 = 7.0_rt*b5*b2;
zeroj = b0 - (t1/b1) - (t3/b3) - (t5/b5) - (t7/b7);
const amrex::Real errj = std::abs(jn(n, zeroj));
// improve solution using procedure SecantRootFinder
if (errj > tol) SecantRootFinder(n, nitmx, tol, &zeroj, &ierror);
roots[ik] = zeroj;
ier[ik] = ierror;
k++;
ik++;
}
}
void SecantRootFinder(int n, int nitmx, amrex::Real tol, amrex::Real *zeroj, int *ier) {
using namespace amrex::literals;
amrex::Real p0, p1, q0, q1, dp, p;
amrex::Real c[2];
c[0] = 0.95_rt;
c[1] = 0.999_rt;
*ier = 0;
p = *zeroj;
for (int ntry=0 ; ntry <= 1 ; ntry++) {
p0 = c[ntry]*(*zeroj);
p1 = *zeroj;
q0 = jn(n, p0);
q1 = jn(n, p1);
for (int it=1; it <= nitmx; it++) {
if (q1 == q0) break;
p = p1 - q1*(p1 - p0)/(q1 - q0);
dp = p - p1;
if (it > 1 && std::abs(dp) < tol) {
*zeroj = p;
return;
}
p0 = p1;
q0 = q1;
p1 = p;
q1 = jn(n, p1);
}
}
*ier = 3;
*zeroj = p;
}
|
