{"id":2814,"date":"2017-11-06T12:00:00","date_gmt":"2017-11-06T17:00:00","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=2814"},"modified":"2018-03-19T12:22:56","modified_gmt":"2018-03-19T17:22:56","slug":"three-term-recurrence-relations-and-bessel-functions","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2017\/11\/06\/three-term-recurrence-relations-and-bessel-functions\/","title":{"rendered":"Three-Term Recurrence Relations and Bessel Functions"},"content":{"rendered":"<div class=\"content\">\r\n<p>Three-term recurrence relations are the basis for computing Bessel functions.<\/p>\r\n<h3>Contents<\/h3>\r\n<div>\r\n<ul>\r\n \t<li><a href=\"#c3a6c2ca-9bbb-4551-a4ce-7e1ece2f4b32\">A Familiar Three-Term Recurrence<\/a><\/li>\r\n \t<li><a href=\"#165c9689-8afb-455b-b9de-0727a37e1d54\">Friedrich Bessel<\/a><\/li>\r\n \t<li><a href=\"#2aca0abf-9cf5-4ba7-9b59-d37efbd348b8\">Bessel Functions<\/a><\/li>\r\n \t<li><a href=\"#ea1545b8-fbc6-4896-be9f-6cc936f7d117\">Miller's Algorithm<\/a><\/li>\r\n \t<li><a href=\"#107cc989-46d5-476c-9284-0f2e15fc0bc7\">Matrix Formulation<\/a><\/li>\r\n \t<li><a href=\"#2bdbb3c9-0d83-4b7b-bfc7-8c44f0c108cf\">Lower<\/a><\/li>\r\n \t<li><a href=\"#44ed1600-aec1-4b74-b0ae-6f0686dee44d\">Upper<\/a><\/li>\r\n \t<li><a href=\"#1605c92d-d833-477c-9ab9-f2f13f574541\">Center<\/a><\/li>\r\n \t<li><a href=\"#a2dc4154-57a7-4378-a774-56f9f61917fc\">Relative error<\/a><\/li>\r\n \t<li><a href=\"#e618929e-98d7-472d-b67f-4788f45d098d\">Triangular Factors<\/a><\/li>\r\n \t<li><a href=\"#956e6b1e-7083-4e19-8fbf-78d708405bfa\">References<\/a><\/li>\r\n<\/ul>\r\n<\/div>\r\n<h4>A Familiar Three-Term Recurrence<a name=\"c3a6c2ca-9bbb-4551-a4ce-7e1ece2f4b32\"><\/a><\/h4>\r\nStart with two large integers, say 514229 and 317811. For reasons that will soon become clear, I'll label them $f_{30}$ and $f_{29}$.\r\n<pre class=\"codeinput\">   clear\r\n   n = 30\r\n   f(30) = 514229;\r\n   f(29) = 317811\r\n<\/pre>\r\n<pre class=\"codeoutput\">n =\r\n    30\r\nf =\r\n  Columns 1 through 6\r\n           0           0           0           0           0           0\r\n  Columns 7 through 12\r\n           0           0           0           0           0           0\r\n  Columns 13 through 18\r\n           0           0           0           0           0           0\r\n  Columns 19 through 24\r\n           0           0           0           0           0           0\r\n  Columns 25 through 30\r\n           0           0           0           0      317811      514229\r\n<\/pre>\r\nNow take repeated differences, and count backwards.\r\n\r\n$$ f_{n-1} = f_{n+1} - f_n, \\ n = 29, 28, ..., 2, 1 $$\r\n<pre class=\"codeinput\">   <span class=\"keyword\">for<\/span> n = 29:-1:2\r\n       f(n-1) = f(n+1) - f(n);\r\n   <span class=\"keyword\">end<\/span>\r\n   n\r\n   f\r\n<\/pre>\r\n<pre class=\"codeoutput\">n =\r\n     2\r\nf =\r\n  Columns 1 through 6\r\n           0           1           1           2           3           5\r\n  Columns 7 through 12\r\n           8          13          21          34          55          89\r\n  Columns 13 through 18\r\n         144         233         377         610         987        1597\r\n  Columns 19 through 24\r\n        2584        4181        6765       10946       17711       28657\r\n  Columns 25 through 30\r\n       46368       75025      121393      196418      317811      514229\r\n<\/pre>\r\nRecognize those integers? They're the venerable <a href=\"https:\/\/blogs.mathworks.com\/cleve\/2012\/06\/03\/fibonacci-matrices\/\">Fibonacci numbers<\/a>. I've just reversed the usual procedure for computing them. And, I didn't pick the first two values at random. In fact, their choice is very delicate. Try picking any other six-digit integers, and see if you can take 30 differences before you hit a negative value. Of course, we usually start with $f_0 = 0$, $f_1 = 1$ and go forward with additions.\r\n\r\nThe relation\r\n\r\n$$ f_n = f_{n+1} - f_{n-1} $$\r\n\r\nis an example of a <i>three-term recurrence<\/i>. I've written it in such a way that I can go either forward or backward.\r\n<h4>Friedrich Bessel<a name=\"165c9689-8afb-455b-b9de-0727a37e1d54\"><\/a><\/h4>\r\nFriedrich Bessel was a German astronomer, mathematician and physicist who lived from 1784 until 1846. He was a contemporary of Gauss, but the two had some sort of falling out. Although he did not have a university education, he made major contributions to precise measurements of planetary orbits and distances to stars. He has a crater on the Moon and an asteroid named after him.\r\n\r\nIn mathematics, the special functions that bear his name were actually discovered by someone else, Daniel Bernoulli, and christened after Bessel's death. They play the same role in polar coordinates that the trig functions play in Cartesian coordinates. They are applicable in a wide range of fields, from physics to finance.\r\n<h4>Bessel Functions<a name=\"2aca0abf-9cf5-4ba7-9b59-d37efbd348b8\"><\/a><\/h4>\r\nBessel functions are defined by an ordinary differential equation that generalizes the harmonic oscillator to polar coordinates. The equation contains a parameter $n$, the <i>order<\/i>.\r\n\r\n$$ x^2 y'' + x y' + (x^2 - n^2) y = 0 $$\r\n\r\nSolutions to this equation are also sometimes called <i>cylindrical harmonics<\/i>.\r\n\r\nFor my purposes here, their most important property is the three-term recurrence relation:\r\n\r\n$$ \\frac{2n}{x} J_n(x) = J_{n-1}(x) + J_{n+1}(x) $$\r\n\r\nThis has a variable coefficient that depends upon the two parameters, the order $n$ and the argument $x$.\r\n\r\nTo complete the specification, it would be necessary to supply two initial conditions that depend upon $x$. But, as we shall see, using the recurrence in the forward direction is disastrous numerically. Alternatively, if we could somehow supply two end conditions, we could use the recurrence in the reverse direction. This the idea behind a classical method known as Miller's algorithm.\r\n<h4>Miller's Algorithm<a name=\"ea1545b8-fbc6-4896-be9f-6cc936f7d117\"><\/a><\/h4>\r\nMiller's algorithm also employs this interesting identity, which I'll call the <i>normalizer<\/i>.\r\n\r\n$$ J_0(x) + 2 J_2(x) + 2 J_4(x) + 2 J_6(x) + \\ ... = 1 $$\r\n\r\nIt's not obvious where this identity comes from and I won't try to derive it here.\r\n\r\nWe don't know the end conditions, but we can pick arbitrary values, run the recurrence backwards, and then use the normalizer to rescale.\r\n\r\nStep 1. Pick a large $N$ and temporarily let\r\n\r\n$$ \\tilde{J}_N = 0 $$\r\n\r\n$$ \\tilde{J}_{N-1} = 1 $$\r\n\r\nThen for $n = N-1, ..., 1$ compute\r\n\r\n$$ \\tilde{J}_{n-1} = 2n\/x \\tilde{J}_n \\ - \\tilde{J}_{n+1} $$\r\n\r\nStep 2. Normalize:\r\n\r\n$$ \\sigma = \\tilde{J}_0 + 2 \\tilde{J}_2 + 2 \\tilde{J}_4 + ... $$\r\n\r\n$$ J_n = \\tilde{J}_n\/\\sigma, \\ n = 0, ..., N $$\r\n\r\nI'll show some results in a moment, after I introduce an alternative formulation.\r\n<h4>Matrix Formulation<a name=\"107cc989-46d5-476c-9284-0f2e15fc0bc7\"><\/a><\/h4>\r\nYou know that I like to describe algorithms in matrix terms. The three-term recurrence generates a tridiagonal matrix. But which three diagonals? I can put the three below the diagonal, above the diagonal, or centered on the diagonal. This is accomplished with the second argument of the sparse matrix generator, <tt>spdiags<\/tt>. Here is the code.\r\n<pre class=\"codeinput\">   type <span class=\"string\">blt<\/span>\r\n<\/pre>\r\n<pre class=\"codeoutput\">function B = blt(n,x,loc)\r\n% BesseL Three term recurence.\r\n% B = blt(n,x,loc) is an n-by-n sparse tridiagonal coefficent matrix\r\n% that generates the Bessel functions besselj((0:n-1)',x). \r\n% loc specifies the location of the diagonals.\r\n%   loc = 'center', the default, centers the three diagonals.\r\n%   loc = 'lower' puts the three diagonals in the lower triangle.\r\n%   loc = 'upper' puts the three diagonals in the upper triangle.\r\n\r\n   if nargin == 2\r\n       loc = 'center';\r\n   end\r\n   \r\n   switch loc\r\n       case 'center', locd = -1:1;\r\n       case 'lower', locd = -2:0;\r\n       case 'upper', locd = 0:2;\r\n   end\r\n   \r\n   % Three term recurrence: 2*n\/x * J_n(x) = J_(n-1)(x) + J_n+1(x).\r\n   d = -2*(0:n-1)'\/x;\r\n   e = ones(n,1);\r\n   B = spdiags([e d e],locd,n,n);\r\nend\r\n<\/pre>\r\n<h4>Lower<a name=\"2bdbb3c9-0d83-4b7b-bfc7-8c44f0c108cf\"><\/a><\/h4>\r\nIf there were no roundoff error, and if we knew the values of the first two Bessel functions, $J_0(x)$ and $J_1(x)$, we could use backslash, the linear equation solver, with the diagonals in the lower triangle to compute $J_n(x)$ for $n &gt; 1$. Here, for example, is $x = 1$.\r\n<pre class=\"codeinput\">   format <span class=\"string\">long<\/span>\r\n   J = @(n,x) besselj(n,x);\r\n   n = 8\r\n   x = 1.0\r\n   B = blt(n,x,<span class=\"string\">'lower'<\/span>);\r\n   rhs = zeros(n,1);\r\n   rhs(1:2) = J((0:1)',x);\r\n   v = B\\rhs\r\n<\/pre>\r\n<pre class=\"codeoutput\">n =\r\n     8\r\nx =\r\n     1\r\nv =\r\n   0.765197686557967\r\n   0.440050585744934\r\n   0.114903484931901\r\n   0.019563353982668\r\n   0.002476638964110\r\n   0.000249757730213\r\n   0.000020938338022\r\n   0.000001502326053\r\n<\/pre>\r\nThese values are accurate to almost all the digits shown. Except for the last one. It's beginning to show the effects of severe numerical instability. Let's try a larger value of <tt>n<\/tt>.\r\n<pre class=\"codeinput\">   n = 18\r\n   B = blt(n,x,<span class=\"string\">'lower'<\/span>);\r\n   rhs = zeros(n,1);\r\n   rhs(1:2) = J((0:1)',x);\r\n   v = B\\rhs;\r\n   semilogy(0:n-1,v,<span class=\"string\">'o-'<\/span>)\r\n   title(<span class=\"string\">'J_n(1)'<\/span>)\r\n   xlabel(<span class=\"string\">'n'<\/span>)\r\n<\/pre>\r\n<pre class=\"codeoutput\">n =\r\n    18\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/bessel_blog2_01.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n\r\nThe values beyond about <tt>n = 9<\/tt> are suspect and, in fact, totally wrong. The forward elimination that backslash uses to solve this lower triangular system is just the three-term recurrence for increasing <tt>n<\/tt>. And that recurrence has two solutions -- the one we're trying to compute and a second one corresponding to $Y_n(x)$, the <i>Bessel function of second kind<\/i>. $Y_n(x)$ is stimulated by roundoff error and completely takes over for $n &gt; 9$.\r\n\r\nSo the lower triangular option -- the forward recurrence -- is unsatisfactory for two reasons: it requires two Bessel values to get started and it is unstable.\r\n<h4>Upper<a name=\"44ed1600-aec1-4b74-b0ae-6f0686dee44d\"><\/a><\/h4>\r\nWe need to make one modification to the upper triangular coefficient matrix.\r\n<pre class=\"codeinput\">   n = 30;\r\n   B = blt(n,x,<span class=\"string\">'upper'<\/span>);\r\n   B(n-1,n) = 0;\r\n<\/pre>\r\nNow the last two rows and columns are a 2-by-2 identity.\r\n<pre class=\"codeinput\">   full(B(n-1:n,n-1:n))\r\n<\/pre>\r\n<pre class=\"codeoutput\">ans =\r\n     1     0\r\n     0     1\r\n<\/pre>\r\nThe back substitution process just runs the three-term recursion backwards. We need to start with $J_{n-1}(x)$ and $J_{n-2}(x)$.\r\n<pre class=\"codeinput\">   rhs = zeros(n,1);\r\n   rhs(n-1:n) = J((n-2:n-1)',x);\r\n   format <span class=\"string\">long<\/span> <span class=\"string\">e<\/span>\r\n   v = B\\rhs\r\n<\/pre>\r\n<pre class=\"codeoutput\">v =\r\n     7.651976865579656e-01\r\n     4.400505857449330e-01\r\n     1.149034849319004e-01\r\n     1.956335398266838e-02\r\n     2.476638964109952e-03\r\n     2.497577302112342e-04\r\n     2.093833800238925e-05\r\n     1.502325817436807e-06\r\n     9.422344172604491e-08\r\n     5.249250179911870e-09\r\n     2.630615123687451e-10\r\n     1.198006746303136e-11\r\n     4.999718179448401e-13\r\n     1.925616764480172e-14\r\n     6.885408200044221e-16\r\n     2.297531532210343e-17\r\n     7.186396586807488e-19\r\n     2.115375568053260e-20\r\n     5.880344573595754e-22\r\n     1.548478441211652e-23\r\n     3.873503008524655e-25\r\n     9.227621982096665e-27\r\n     2.098223955943776e-28\r\n     4.563424055950103e-30\r\n     9.511097932712488e-32\r\n     1.902951751891381e-33\r\n     3.660826744416801e-35\r\n     6.781552053554108e-37\r\n     1.211364502417112e-38\r\n     2.089159981718163e-40\r\n<\/pre>\r\nThese values are accurate to almost all the significant digits shown. The backwards recursion suppresses the unwanted $Y_n(x)$ and the process is stable.\r\n<h4>Center<a name=\"1605c92d-d833-477c-9ab9-f2f13f574541\"><\/a><\/h4>\r\nCan we avoid having to know those two starting values? Yes, by using a matrix formulation of the Miller algorithm. Insert the normalization condition in the first row. Now $B$ is tridiagonal, except for first row, which is\r\n<pre class=\"language-matlab\">[1 0 2 0 2 0 2 <span class=\"keyword\">...<\/span><span class=\"comment\"> ]<\/span>\r\n<\/pre>\r\nOther rows have $1$ s on the super- and subdiagonal and $-2n\/x$ on the diagonal of the $(n+1)$ -st row.\r\n<pre class=\"codeinput\">   n = 30;\r\n   B = blt(n,x);\r\n   B(1,1:2) = [1 0];\r\n   B(1,3:2:n) = 2;\r\n   spy(B)\r\n   title(<span class=\"string\">'B'<\/span>)\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/bessel_blog2_02.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n\r\nThe beauty of this approach is that it does not require any values of the Bessel function <i>a priori<\/i>. The right hand side is simply a unit vector.\r\n<pre class=\"codeinput\">   rhs = eye(n,1);\r\n   format <span class=\"string\">long<\/span> <span class=\"string\">e<\/span>\r\n   v = B\\rhs\r\n   semilogy(0:n-1,v,<span class=\"string\">'o-'<\/span>)\r\n   xlabel(<span class=\"string\">'n'<\/span>)\r\n   title(<span class=\"string\">'J_n(1)'<\/span>)\r\n<\/pre>\r\n<pre class=\"codeoutput\">v =\r\n     7.651976865579665e-01\r\n     4.400505857449335e-01\r\n     1.149034849319005e-01\r\n     1.956335398266840e-02\r\n     2.476638964109954e-03\r\n     2.497577302112343e-04\r\n     2.093833800238926e-05\r\n     1.502325817436807e-06\r\n     9.422344172604494e-08\r\n     5.249250179911872e-09\r\n     2.630615123687451e-10\r\n     1.198006746303136e-11\r\n     4.999718179448401e-13\r\n     1.925616764480171e-14\r\n     6.885408200044220e-16\r\n     2.297531532210343e-17\r\n     7.186396586807487e-19\r\n     2.115375568053260e-20\r\n     5.880344573595753e-22\r\n     1.548478441211652e-23\r\n     3.873503008524655e-25\r\n     9.227621982096663e-27\r\n     2.098223955943775e-28\r\n     4.563424055950102e-30\r\n     9.511097932712487e-32\r\n     1.902951751891381e-33\r\n     3.660826744416762e-35\r\n     6.781552053355331e-37\r\n     1.211364395117367e-38\r\n     2.088559301926495e-40\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/bessel_blog2_03.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n<h4>Relative error<a name=\"a2dc4154-57a7-4378-a774-56f9f61917fc\"><\/a><\/h4>\r\n<pre class=\"codeinput\">   w = J((0:n-1)',x);\r\n   relerr = (v - w).\/w;\r\n   semilogy(0:n-1,abs(relerr) + eps*(v==w),<span class=\"string\">'o-'<\/span>)\r\n   xlabel(<span class=\"string\">'n'<\/span>)\r\n   title(<span class=\"string\">'relative error'<\/span>)\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/bessel_blog2_04-1.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n\r\nThe relative error deteriorates for the last few values of <tt>n<\/tt>. This is not unexpected because we have effectively truncated a semi-infinite matrix.\r\n<h4>Triangular Factors<a name=\"e618929e-98d7-472d-b67f-4788f45d098d\"><\/a><\/h4>\r\nThe LU decomposition of $B$ shows no pivoting for this value of $x$, but there is fill-in. The normalization condition is mixing in with the recurrence.\r\n<pre class=\"codeinput\">    [L,U] = lu(B);\r\n    spy(L)\r\n    title(<span class=\"string\">'L'<\/span>)\r\n    snapnow\r\n    spy(U)\r\n    title(<span class=\"string\">'U'<\/span>)\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/bessel_blog2_05.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/> <img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/bessel_blog2_06.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n<h4>References<a name=\"956e6b1e-7083-4e19-8fbf-78d708405bfa\"><\/a><\/h4>\r\nJ. C. P. Miller, \"On the choice of standard solutions for a homogeneous linear equation of the second order\", <i>Quart. J. Mech. Appl. Math.<\/i> 3, 1950, 225-235.\r\n\r\nWalter Gautschi, \"Computational Aspects of Three-Term Recurrence Relations\", <i>SIAM Review<\/i> 9, 1967, 24-82.\r\n\r\n<script language=\"JavaScript\"> <!-- \r\n    function grabCode_5da534de2d7a43c0ad4c7a4c28cb0e39() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='5da534de2d7a43c0ad4c7a4c28cb0e39 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 5da534de2d7a43c0ad4c7a4c28cb0e39';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2017 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<\/p>\r\n<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\r\n<p>\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }   \r\n     --> <\/script>\r\n<p style=\"text-align: right; font-size: xx-small; font-weight: lighter; font-style: italic; color: gray;\">\r\n<a><span style=\"font-size: x-small; font-style: italic;\">Get\r\nthe MATLAB code<noscript>(requires JavaScript)<\/noscript><\/span><\/a>\r\n\r\nPublished with MATLAB\u00ae R2017a<\/p>\r\n\r\n<\/div>\r\n<!--\r\n5da534de2d7a43c0ad4c7a4c28cb0e39 ##### SOURCE BEGIN #####\r\n%% Three-Term Recurrence Relations and Bessel Functions\r\n% Three-term recurrence relations are the basis for computing\r\n% Bessel functions.\r\n\r\n%% A Familiar Three-Term Recurrence\r\n% Start with two large integers, say 514229 and 317811.\r\n% For reasons that will soon become clear, I'll label them $f_{30}$\r\n% and $f_{29}$.\r\n\r\nclear\r\nn = 30\r\nf(30) = 514229;\r\nf(29) = 317811\r\n\r\n%%\r\n% Now take repeated differences, and count backwards.\r\n%\r\n% $$ f_{n-1} = f_{n+1} - f_n, \\ n = 29, 28, ..., 2, 1 $$\r\n\r\nfor n = 29:-1:2\r\nf(n-1) = f(n+1) - f(n);\r\nend\r\nn\r\nf\r\n\r\n%%\r\n% Recognize those integers?  They're the venerable\r\n% <https:\/\/blogs.mathworks.com\/cleve\/2012\/06\/03\/fibonacci-matrices\/ % Fibonacci numbers>.  I've just reversed the usual procedure for\r\n% computing them.  And, I didn't pick the first two values\r\n% at random.  In fact, their choice is very delicate.  Try picking\r\n% any other six-digit integers, and see if you can take 30 differences\r\n% before you hit a negative value.  Of course, we usually start with\r\n% $f_0 = 0$, $f_1 = 1$ and go forward with additions.\r\n\r\n%%\r\n% The relation\r\n%\r\n% $$ f_n = f_{n+1} - f_{n-1} $$\r\n%\r\n% is an example of a _three-term recurrence_.  I've written it in such\r\n% a way that I can go either forward or backward.\r\n\r\n%% Friedrich Bessel\r\n% Friedrich Bessel was a German astronomer, mathematician and physicist\r\n% who lived from 1784 until 1846.\r\n% He was a contemporary of Gauss, but the two had some sort of falling\r\n% out.  Although he did not have a university education, he made major\r\n% contributions to precise measurements of planetary orbits and\r\n% distances to stars.  He has a crater on the Moon and an asteroid\r\n% named after him.\r\n\r\n%%\r\n% In mathematics, the special functions that bear his name were actually\r\n% discovered by someone else, Daniel Bernoulli, and christened after\r\n% Bessel's death.  They play the same role in polar coordinates that the\r\n% trig functions play in Cartesian coordinates.  They are applicable\r\n% in a wide range of fields, from physics to finance.\r\n\r\n%% Bessel Functions\r\n% Bessel functions are defined by an ordinary differential\r\n% equation that generalizes the harmonic oscillator to polar\r\n% coordinates.  The equation contains a parameter $n$, the _order_.\r\n%\r\n% $$ x^2 y'' + x y' + (x^2 - n^2) y = 0 $$\r\n%\r\n% Solutions to this equation are also sometimes called\r\n% _cylindrical harmonics_.\r\n\r\n%%\r\n% For my purposes here, their most important property is the\r\n% three-term recurrence relation:\r\n%\r\n% $$ \\frac{2n}{x} J_n(x) = J_{n-1}(x) + J_{n+1}(x) $$\r\n%\r\n% This has a variable coefficient that depends upon the two parameters,\r\n% the order $n$ and the argument $x$.\r\n\r\n%%\r\n% To complete the specification, it would be necessary to supply two\r\n% initial conditions that depend upon $x$.\r\n% But, as we shall see, using the recurrence in the forward direction\r\n% is disastrous numerically.\r\n% Alternatively, if we could somehow supply two end conditions,\r\n% we could use the recurrence in the reverse direction.\r\n% This the idea behind a classical method known as Miller's algorithm.\r\n\r\n%% Miller's Algorithm\r\n% Miller's algorithm also employs this interesting identity,\r\n% which I'll call the _normalizer_.\r\n%\r\n% $$ J_0(x) + 2 J_2(x) + 2 J_4(x) + 2 J_6(x) + \\ ... = 1 $$\r\n%\r\n% It's not obvious where this identity comes from and I won't try to\r\n% derive it here.\r\n\r\n%%\r\n% We don't know the end conditions, but we can pick arbitrary values,\r\n% run the recurrence backwards, and then use the normalizer to rescale.\r\n\r\n%%\r\n% Step 1. Pick a large $N$ and temporarily let\r\n%\r\n% $$ \\tilde{J}_N = 0 $$\r\n%\r\n% $$ \\tilde{J}_{N-1} = 1 $$\r\n%\r\n% Then for $n = N-1, ..., 1$ compute\r\n%\r\n% $$ \\tilde{J}_{n-1} = 2n\/x \\tilde{J}_n \\ - \\tilde{J}_{n+1} $$\r\n%\r\n\r\n%%\r\n% Step 2.  Normalize:\r\n%\r\n% $$ \\sigma = \\tilde{J}_0 + 2 \\tilde{J}_2 + 2 \\tilde{J}_4 + ... $$\r\n%\r\n% $$ J_n = \\tilde{J}_n\/\\sigma, \\ n = 0, ..., N $$\r\n\r\n%%\r\n% I'll show some results in a moment, after I introduce an\r\n% alternative formulation.\r\n\r\n%% Matrix Formulation\r\n% You know that I like to describe algorithms in matrix terms.\r\n% The three-term recurrence generates a tridiagonal matrix.\r\n% But which three diagonals?  I can put the three below the diagonal,\r\n% above the diagonal, or centered on the diagonal.  This is accomplished\r\n% with the second argument of the sparse matrix generator, |spdiags|.\r\n% Here is the code.\r\n\r\ntype blt\r\n\r\n%% Lower\r\n% If there were no roundoff error, and if we knew the values of the\r\n% first two Bessel functions, $J_0(x)$ and $J_1(x)$, we could use\r\n% backslash, the linear equation solver, with the diagonals in the\r\n% lower triangle to compute $J_n(x)$ for $n > 1$.  Here, for example,\r\n% is $x = 1$.\r\n\r\nformat long\r\nJ = @(n,x) besselj(n,x);\r\nn = 8\r\nx = 1.0\r\nB = blt(n,x,'lower');\r\nrhs = zeros(n,1);\r\nrhs(1:2) = J((0:1)',x);\r\nv = B\\rhs\r\n\r\n%%\r\n% These values are accurate to almost all the digits shown.\r\n% Except for the last one.  It's beginning to show the effects of\r\n% severe numerical instability.  Let's try a larger value of |n|.\r\n\r\nn = 18\r\nB = blt(n,x,'lower');\r\nrhs = zeros(n,1);\r\nrhs(1:2) = J((0:1)',x);\r\nv = B\\rhs;\r\nsemilogy(0:n-1,v,'o-')\r\ntitle('J_n(1)')\r\nxlabel('n')\r\n\r\n%%\r\n% The values beyond about |n = 9| are suspect and, in fact, totally\r\n% wrong.  The forward elimination that backslash uses to solve this\r\n% lower triangular system is just the three-term recurrence for\r\n% increasing |n|.  And that recurrence has two solutions REPLACE_WITH_DASH_DASH the one\r\n% we're trying to compute and a second one corresponding to $Y_n(x)$,\r\n% the _Bessel function of second kind_.  $Y_n(x)$ is stimulated\r\n% by roundoff error and completely takes over for $n > 9$.\r\n\r\n%%\r\n% So the lower triangular option REPLACE_WITH_DASH_DASH the forward recurrence REPLACE_WITH_DASH_DASH is\r\n% unsatisfactory for two reasons: it requires two Bessel values to get\r\n% started and it is unstable.\r\n\r\n%% Upper\r\n% We need to make one modification to the upper triangular coefficient\r\n% matrix.\r\n\r\nn = 30;\r\nB = blt(n,x,'upper');\r\nB(n-1,n) = 0;\r\n\r\n%%\r\n% Now the last two rows and columns are a 2-by-2 identity.\r\n\r\nfull(B(n-1:n,n-1:n))\r\n\r\n%%\r\n% The back substitution process just runs the three-term recursion\r\n% backwards.  We need to start with $J_{n-1}(x)$ and $J_{n-2}(x)$.\r\n\r\nrhs = zeros(n,1);\r\nrhs(n-1:n) = J((n-2:n-1)',x);\r\nformat long e\r\nv = B\\rhs\r\n\r\n%%\r\n% These values are accurate to almost all the significant digits\r\n% shown.  The backwards recursion suppresses the unwanted $Y_n(x)$ and\r\n% the process is stable.\r\n\r\n%% Center\r\n% Can we avoid having to know those two starting values?\r\n% Yes, by using a matrix formulation of the Miller algorithm.\r\n% Insert the normalization condition in the first row.\r\n% Now $B$ is tridiagonal, except for first row, which is\r\n%\r\n%   [1 0 2 0 2 0 2 ... ]\r\n%\r\n% Other rows have $1$ s on the super- and subdiagonal and $-2n\/x$\r\n% on the diagonal of the $(n+1)$ -st row.\r\n\r\nn = 30;\r\nB = blt(n,x);\r\nB(1,1:2) = [1 0];\r\nB(1,3:2:n) = 2;\r\nspy(B)\r\ntitle('B')\r\n\r\n%%\r\n% The beauty of this approach is that it does not require any\r\n% values of the Bessel function _a priori_.  The right hand side\r\n% is simply a unit vector.\r\n\r\nrhs = eye(n,1);\r\nformat long e\r\nv = B\\rhs\r\nsemilogy(0:n-1,v,'o-')\r\nxlabel('n')\r\ntitle('J_n(1)')\r\n\r\n%% Relative error\r\n\r\nw = J((0:n-1)',x);\r\nrelerr = (v - w).\/w;\r\nsemilogy(0:n-1,abs(relerr) + eps*(v==w),'o-')\r\nxlabel('n')\r\ntitle('relative error')\r\n\r\n%%\r\n% The relative error deteriorates for the last few values of |n|.\r\n% This is not unexpected because we have effectively truncated\r\n% a semi-infinite matrix.\r\n\r\n%% Triangular Factors\r\n% The LU decomposition of $B$ shows no pivoting for this value of $x$,\r\n% but there is fill-in.  The normalization condition is mixing in\r\n% with the recurrence.\r\n\r\n[L,U] = lu(B);\r\nspy(L)\r\ntitle('L')\r\nsnapnow\r\nspy(U)\r\ntitle('U')\r\n\r\n%% References\r\n% J. C. P. Miller, \"On the choice of standard solutions for a\r\n% homogeneous linear equation of the second order\",\r\n% _Quart. J. Mech. Appl. Math._ 3, 1950, 225-235.\r\n%\r\n% Walter Gautschi, \"Computational Aspects of Three-Term Recurrence\r\n% Relations\", _SIAM Review_ 9, 1967, 24-82.\r\n%\r\n\r\n##### SOURCE END ##### 5da534de2d7a43c0ad4c7a4c28cb0e39\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/bessel_blog2_06.png\" class=\"img-responsive attachment-post-thumbnail size-post-thumbnail wp-post-image\" alt=\"\" decoding=\"async\" loading=\"lazy\" \/><\/div><p>\r\nThree-term recurrence relations are the basis for computing Bessel functions.\r\nContents\r\n\r\n\r\n \tA Familiar Three-Term Recurrence\r\n \tFriedrich Bessel\r\n \tBessel Functions\r\n \tMiller's Algorithm\r\n... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2017\/11\/06\/three-term-recurrence-relations-and-bessel-functions\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":2821,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[11,24,8,17],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/2814"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/users\/78"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/comments?post=2814"}],"version-history":[{"count":11,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/2814\/revisions"}],"predecessor-version":[{"id":3112,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/2814\/revisions\/3112"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media\/2821"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=2814"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=2814"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=2814"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}