forked from chebfun/chebfun
-
Notifications
You must be signed in to change notification settings - Fork 1
/
besselroots.m
118 lines (108 loc) · 6.14 KB
/
besselroots.m
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
function j = besselroots(v, n)
%BESSELROOTS Roots of the function J_v(x) = besselj(v, x).
% BESSELROOTS(V, N) returns the first N roots the Bessel function J_V(X).
% Both V and N must be scalars, and N must be non-negative.
%
% The accuracy of the returned roots can be described roughly as follows:
%
% V = 0 --> Full double precision for N <= 20 (Wolfram Alpha), and very
% accurate approximations for N > 20 (McMahon's expansion);
% -1 <= V <= 5 : V ~= 0 --> 12 decimal figures for the 6 first zeros
% (Piessens's Chebyshev series approximations), and very accurate
% approximations for the others (McMahon's expansion);
% V > 5 --> moderately accurate for the 6 first zeros and good
% approximations for the others (McMahon's expansion).
% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.
% L. L. Peixoto, 2015
% TODO: Improve accuracy to full machine precision for all roots?
% Store to 16 digits for integer (and half-integer) v in [-1,5]?
% Trivial case:
if ( n == 0 )
j = [];
return
end
% Check inputs:
if ( ~isscalar(n) || (round(n) - n ~= 0) || n < 0 )
error('CHEBFUN:besselroots:inputN', 'Input N must be a positive integer');
end
% McMahon's expansion. This expansion gives very accurate approximation
% for the sth zero (s >= 7) in the whole region V >=- 1, and moderate
% approximation in other cases.
s = (1:n)';
mu = 4*v^2;
a1 = 1 / 8;
a3 = (7*mu-31) / 384;
a5 = 4*(3779+mu*(-982+83*mu)) / 61440; % Evaluate via Horner's method.
a7 = 6*(-6277237+mu*(1585743+mu*(-153855+6949*mu))) / 20643840;
a9 = 144*(2092163573+mu*(-512062548+mu*(48010494+mu*(-2479316+70197*mu)))) ...
/ 11890851840;
a11 = 720*(-8249725736393+mu*(1982611456181+mu*(-179289628602+mu*(8903961290 + ...
mu*(-287149133+5592657*mu))))) / 10463949619200;
a13 = 576*(423748443625564327 + mu*(-100847472093088506+mu*(8929489333108377 + ...
mu*(-426353946885548+mu*(13172003634537+mu*(-291245357370 + mu*4148944183)))))) ...
/ 13059009124761600;
b = .25*(2*v+4*s-1)*pi; % beta
j = b - (mu-1)*polyval([a13 0 a11 0 a9 0 a7 0 a5 0 a3 0 a1 0], 1./b);
if ( v == 0 )
% First 20 roots of J0(x) are precomputed (using Wolfram Alpha):
j(1:20) = [
2.4048255576957728
5.5200781102863106
8.6537279129110122
11.791534439014281
14.930917708487785
18.071063967910922
21.211636629879258
24.352471530749302
27.493479132040254
30.634606468431975
33.775820213573568
36.917098353664044
40.058425764628239
43.199791713176730
46.341188371661814
49.482609897397817
52.624051841114996
55.765510755019979
58.906983926080942
62.048469190227170];
elseif ( v >= -1 && v <= 5 )
% Piessens's Chebyshev series approximations (1984). Calculates the 6 first
% zeros to at least 12 decimal figures in region -1 <= V <= 5:
C = [
2.883975316228 8.263194332307 11.493871452173 14.689036505931 17.866882871378 21.034784308088
0.767665211539 4.209200330779 4.317988625384 4.387437455306 4.435717974422 4.471319438161
-0.086538804759 -0.164644722483 -0.130667664397 -0.109469595763 -0.094492317231 -0.083234240394
0.020433979038 0.039764618826 0.023009510531 0.015359574754 0.011070071951 0.008388073020
-0.006103761347 -0.011799527177 -0.004987164201 -0.002655024938 -0.001598668225 -0.001042443435
0.002046841322 0.003893555229 0.001204453026 0.000511852711 0.000257620149 0.000144611721
-0.000734476579 -0.001369989689 -0.000310786051 -0.000105522473 -0.000044416219 -0.000021469973
0.000275336751 0.000503054700 0.000083834770 0.000022761626 0.000008016197 0.000003337753
-0.000106375704 -0.000190381770 -0.000023343325 -0.000005071979 -0.000001495224 -0.000000536428
0.000042003336 0.000073681222 0.000006655551 0.000001158094 0.000000285903 0.000000088402
-0.000016858623 -0.000029010830 -0.000001932603 -0.000000269480 -0.000000055734 -0.000000014856
0.000006852440 0.000011579131 0.000000569367 0.000000063657 0.000000011033 0.000000002536
-0.000002813300 -0.000004672877 -0.000000169722 -0.000000015222 -0.000000002212 -0.000000000438
0.000001164419 0.000001903082 0.000000051084 0.000000003677 0.000000000448 0.000000000077
-0.000000485189 -0.000000781030 -0.000000015501 -0.000000000896 -0.000000000092 -0.000000000014
0.000000203309 0.000000322648 0.000000004736 0.000000000220 0.000000000019 0.000000000002
-0.000000085602 -0.000000134047 -0.000000001456 -0.000000000054 -0.000000000004 0
0.000000036192 0.000000055969 0.000000000450 0.000000000013 0 0
-0.000000015357 -0.000000023472 -0.000000000140 -0.000000000003 0 0
0.000000006537 0.000000009882 0.000000000043 0.000000000001 0 0
-0.000000002791 -0.000000004175 -0.000000000014 0 0 0
0.000000001194 0.000000001770 0.000000000004 0 0 0
-0.000000000512 -0.000000000752 0 0 0 0
0.000000000220 0.000000000321 0 0 0 0
-0.000000000095 -0.000000000137 0 0 0 0
0.000000000041 0.000000000059 0 0 0 0
-0.000000000018 -0.000000000025 0 0 0 0
0.000000000008 0.000000000011 0 0 0 0
-0.000000000003 -0.000000000005 0 0 0 0
0.000000000001 0.000000000002 0 0 0 0];
j(1:6) = chebtech.clenshaw((v-2)/3, C).'; % Evaluate via Clenshaw.
j(1) = j(1) * sqrt(v+1); % Scale the first root.
end
j = j(1:n); % Trim unnecessary points.
end