https://mizugadro.mydns.jp/t/index.php?title=Discrete_Bessel&feed=atom&action=historyDiscrete Bessel - Revision history2024-03-29T07:31:45ZRevision history for this page on the wikiMediaWiki 1.31.16https://mizugadro.mydns.jp/t/index.php?title=Discrete_Bessel&diff=28280&oldid=prevT: Text replacement - "\$([^\$]+)\$" to "\\(\1\\)"2019-07-30T09:44:22Z<p>Text replacement - "\$([^\$]+)\$" to "\\(\1\\)"</p>
<a href="https://mizugadro.mydns.jp/t/index.php?title=Discrete_Bessel&diff=28280&oldid=14468">Show changes</a>Thttps://mizugadro.mydns.jp/t/index.php?title=Discrete_Bessel&diff=14468&oldid=prevMaintenance script at 22:00, 30 November 20182018-11-30T22:00:01Z<p></p>
<p><b>New page</b></p><div>[[Discrete Bessel]] is one of the discrete approximations of the [[Bessel transform]].<br />
<br />
This article describes the simple implementation in Mathematica language and the testing.<br />
<br />
==[[Bessel transform]]==<br />
<br />
Let the Bessel transform $f$ of function $g$ be defined with<br />
<br />
$\displaystyle<br />
f(x)=\int_0^\infty J_v(xy)\, g(y)\, y\, \mathrm d y<br />
$<br />
<br />
where $J_v$ is the [[Bessel function]] of order $v$. Often, $v$ is assumed to be integer number or just zero.<br />
<br />
Here $x$ may have sense of the radial coordinate (in polar system of coordinates), measured in some appropriate units;<br />
then $y$ has sense of the corresponding radial wavenumber.<br />
<br />
Such a representation takes 21 symbols and seems to be simplest possible.<br />
<br />
==Naive implementation==<br />
This section represents the Mathematica implementation and its test.<br />
<br />
===Example with $~M\!=\!8$===<br />
<poem><br />
M = 8;<br />
S = 1. BesselJZero[0, M+1];<br />
X[n_] = BesselJZero[0,n]/Sqrt[S];<br />
W[n_] = Sqrt[2./S]/Abs[BesselJ[1, BesselJZero[0,n]]];<br />
T = Table[Table[W[m] BesselJ[0, X[m] X[n]] W[n], {n,1,M}], {m,1,M}];<br />
Id = IdentityMatrix[M];<br />
TT = T.T;<br />
TableForm[Table[{X[n], W[n]}, {n,1,M}]]<br />
TableForm[T]<br />
TableForm[TT]<br />
</poem><br />
<br />
This code produces the output; it begins with table of 8 coordinates $x$ of the grid points and the corresponding weights $w$ :<br />
<br />
$\begin{array}{ll}<br />
0.4586366203331863 &0.5195285071552199\\<br />
1.0527624177874753 & 0.7926530133638695\\<br />
1.650396849184917 & 0.9935886501287517\\<br />
2.2488240306434886 &1.1602517418890865\\<br />
2.8475519209198557 & 1.3058173905056825\\<br />
3.446425324924121 & 1.4367104029426196\\<br />
4.0453800503454875 & 1.556635975360147\\<br />
4.644384788693245 & 1.6679612725451483<br />
\end{array}$<br />
<br />
In this approximation<br />
<br />
$ f_n=f(x_n) \, w_n $ <br><br />
$ g_n=g(x_n) \, w_n $<br />
<br />
$ \displaystyle f_m \approx \sum_{n=1}^M T_{m,n} \, g_n $<br />
<br />
$ \displaystyle g_m \approx \sum_{n=1}^M T_{m,n} \, f_n $<br />
<br />
Matrix $T$ is symmetric and almost unitary; $T\!\approx\! T^{-1}$.<br />
<br />
For $M\!=\!8$, the deviation of $T^2$ from the identity matrix is of order of $10^{-7}$.<br />
This deviation quickly drops at the increasing of number $M$ of the nodes.<br />
<br />
At the grid, the scalar product of two functions $\phi$ and $\psi$ is approximated with<br />
<br />
$\phi^{\dagger} \psi=\,$ $\displaystyle<br />
\displaystyle \int_0^\infty \phi(r)^*\, \psi(r)\, r\, \mathrm dr \approx\,$ $\displaystyle \sum_{n=1}^M \, \phi_n^*\, \psi_n=\,$ $\displaystyle<br />
\sum_{n=1}^M \, \phi(x_n)^*\, \psi(x_n)\, {w_n}{^2}$<br />
<br />
===Test with the eigenfunction===<br />
The Bessel function transform has eigenfunctions; they are self–bessels: <br />
such a function is Bessel transform is itself.<br />
<br />
Simplest among these functions is Gaussian, $f(x)=\exp(-x^2/2)$.<br />
<br />
For the coordinates X, weights W and matrix $T$, calculated in the example 1, the test can be written as follows:<br />
<poem><br />
F = Table[{Exp[-X[n]^2/2] W[n]}, {n, 1, M}]<br />
T.F<br />
</poem><br />
The test prints two almost identical lines,<br />
<br />
{{0.467663}, {0.455425}, {0.25453}, {0.0925536}, {0.0226534}, {0.00378554}, {0.000435062}, {0.0000345345}}<br />
<br />
{{0.467663}, {0.455425}, {0.25453}, {0.0925535}, {0.0226534}, {0.00378551}, {0.000435082}, {0.0000344505}}<br />
<br />
The same test can be performed also in C++ with code below: <poem><nomathjax><nowiki><br />
#include<math.h><br />
#include<stdio.h><br />
#define DB double<br />
#define DO(x,y) for(x=0;x<y;x++)<br />
DB jnp(int n,DB x){ return .5*( jn(n-1,x)-jn(n+1,x) ) ; } // Derivative of n th Bessel<br />
DB jnz(int v, int k){ DB x,t; t=M_PI*(k+.5*v-.25); x= t - (v*v-.25)*.5/t;<br />
x-= jn(v,x)/jnp(v,x); // Newton adjustment of the root<br />
x-= jn(v,x)/jnp(v,x);<br />
x-= jn(v,x)/jnp(v,x);<br />
return x; } // the k th zero of v th Bessel<br />
<br />
int main(){ int m,n,v,k; DB s, x,y; int M=8;<br />
DB X[M+1],W[M+1],T[M+1][M+1],TT[M+1][M+1], F[M],G[M];<br />
DB S=jnz(0,M+1);<br />
DB qs=sqrt(1./S);<br />
DB q=sqrt(2./S);<br />
for(n=1;n<M+1;n++){ x=jnz(0,n); X[n]=x*qs; y=W[n]=q/fabs(j1(x));<br />
printf("%3d %20.16lf %20.16lf\n",n,X[n],W[n]); }<br />
for(m=1;m<=M;m++){ printf("%2d",m); for(n=1;n<=M;n++){ T[m][n]=W[m]*j0(X[m]*X[n])*W[n]; }}<br />
for(m=1;m<=M;m++){printf("\n"); for(n=1;n<=M;n++){printf("%14.10lf",T[m][n]); }}<br />
printf("\n");<br />
for(m=1;m<=M;m++){printf("\n");<br />
for(n=1;n<=M;n++){s=0.; for(k=1;k<=M;k++) s+=T[m][k]*T[k][n] ;<br />
TT[m][n]=s;printf("%14.10lf",TT[m][n]); }}<br />
printf("\n\n");<br />
for(m=1;m<=M;m++){x=X[m]; y=exp(-x*x/2.); F[m]=y*W[m]; printf("%14.10lf",F[m]); }<br />
printf("\n");<br />
for(m=1;m<=M;m++){ s=0.; for(n=1;n<=M;n++) s+=T[m][n]*F[n]; G[m]=s; }<br />
for(m=1;m<=M;m++) printf("%14.10lf",G[m]);<br />
printf("\n");<br />
return 0;}<br />
</nowiki></nomathjax></poem><br />
[[File:BesselTestExp200.jpg|300px|thumb| $y=\exp(-x^2/2)$, black solid line; the representation at nodes of the mesh, small red circles; the discrete transform, big blue curcles]]<br />
<br />
The [[C++]] algorithm above is used in the code that generates figure at right.<br />
<br />
The [[discrete Bessel]] on the grid of 8 nodes reproduces 7 decimal digits of the simplest self-Bessel function.<br />
<br />
The same Gaussian is also [[self-Fourier]], and in the similar way can be used for testing of the [[discrete Fourier]] transform.<br />
<br />
===Test with step function===<br />
[[File:Bessel8TestStep.jpg|300px|thumb|$y=\mathrm{UnitStep}(2\!-\!x)$, green line.<br />
The discrete representation, small red circles.<br />
The [[Discrete Bessel]] transform, big blue circles.<br />
The exact [[Bessel transform]] $y\!=\!g(x)$, black line.]]<br />
<br />
Another example with [[discrete Bessel]] transform of the step function is considered in the this section. Let <br />
<br />
$~f(x)=\mathrm{UnitStep}(2\!-\!x)$<br />
<br />
With the same arrays of coordinates X, corresponding weights W and kernel matrix T, as in the previous tests, the [[discrete Bessel]] can be calculated with the [[Mathematica]] code below:<br />
<poem><br />
F = Table[{ W[n] UnitStep[2-X[n]]}, {n, 1, M}]<br />
T.F<br />
G = Table[{ W[n] BesselJ[1, 2 X[n]] 2/X[n]}, {n, 1, M}]<br />
</poem><br />
It produces 3 lines of output: <br><br />
{{0.519529}, {0.792653}, {0.993589}, {0.}, {0.}, {0.}, {0.}, {0.}}<br><br />
{{0.888362}, {0.85255}, {0.316075}, {-0.208464}, {-0.342645}, {-0.0956198}, {0.204077}, {0.240501}}<br><br />
{{0.93354} , {0.85489} , {0.265299}, {-0.237771}, {-0.297814}, {-0.0309065}, {0.189716}, {0.145389}}<br />
<br />
The last two lines show some similarity, but, for the abrupt function, the agreement is worse than for the smooth function in the previous example.<br />
The last table corresponds to the exact integration,<br />
<br />
$\displaystyle<br />
g(x)=\int_0^2 J_0(xy) \,y \, \mathrm d y = 2 J_1(2x)/x$<br />
<br />
The [[C++]] implementation becomes straightforward, as the example with self-Fourier from the previous section is considered; the [[C++]] implementation is used to generate figure at right.<br />
<br />
===Test with truncated Bessel===<br />
[[File:Test8trubes300.jpg|300px|thumb|$y\!=\!J_0(x)$, $y\!=\!g(x)$, black solid lines, and their discrete representations, colored circles]]<br />
Let $\lambda$ be the first zero of the Bessel function,<br />
<br />
$\lambda=\,$ [[BesselJZero]]$[0,1] \approx\, 2.404825557695773$ <br />
<br />
sp, that $~ J_0(\lambda)\!=\!0\,$. Let <br />
<br />
$f(x)=J_0(x)\,$[[UnitStep]]$(\lambda\!-\!x)$<br />
<br />
The [[Bessel transform]] is<br />
<br />
$\displaystyle<br />
g(x)=\int_0^\infty f(y)\, J_0(xy)\, y\, \mathrm d y=$ $\displaystyle<br />
\int_0^\lambda J_0(y)\, J_0(xy)\, y\, \mathrm d y=$ $\displaystyle<br />
\lambda\, J_1(\lambda)\, \frac{ J_0(\lambda x) }{1\!-\!x^2}$<br />
<br />
Note, that the singularity (pole) in denominator of the last expression is compensated by the zero at its numerator;<br />
so, the function $g$ has no singularities at the positive part of the real axis, and its first is<br />
<br />
[[BesselJZero]]$[0,2]/\lambda=\,$[[BesselJZero]]$[0,2]/$[[BesselJZero]]$[0,1] \approx\, 2.2954172674276943$<br />
<br />
==Problems==<br />
==References==<br />
<references/><br />
<br />
==Keywords==<br />
[[Bessel function]]<br />
[[Bessel transform]]<br />
<br />
[[Category:Bessel function]]<br />
[[Category:Bessel transform]]<br />
[[Category:English]]<br />
[[Category:Mathematica]]</div>Maintenance script