Difference between revisions of "File:Test8trubes300.jpg"

From TORI
Jump to: navigation, search
(Importing image file)
 
 
Line 1: Line 1:
  +
Test of the [[Discrete Bessel]] transform of the truncated [[Bessel function]]
Importing image file
 
  +
  +
$f(x)=J_0(x)\,$[[UnitStep]]$(\lambda\!-\!x)$
  +
  +
where $~ \lambda=\,$[[BesselJZero]]$[0,1]\approx 2.404825557695773$
  +
==Description==
  +
  +
Small red circle represent the input function, represented at the grid of $M\!=\!8$ points; the coordinates of these circles are $(x_n,f(x_n))$ for $n=1..8$, where $x_n$ are abscissas of the grid for the [[discrete Bessel]].
  +
  +
These circles follow the smooth profile $y\!=\!J_0(x)$, shown with thin black curve
  +
  +
The input array is determined with $f_n\!=\! f(x_n)\, w_n$ for $n=1.8$
  +
  +
The transform matrix $T$ of the [[Discrete Bessel]] determines the output $\displaystyle g_n=\sum_{m=1}^M T_{n,m}\, f_m$
  +
  +
This output is represented by the big blue circles wit coordinates $(x_n,g_n/w_n)$ for $n=1.8$
  +
  +
It follows the expected smooth curve
  +
  +
$\displaystyle
  +
y=g(x)=$ $\displaystyle
  +
\int_0^\infty f(t) \, J_0(t \, x) \, t \, \mathrm d t=$ $\displaystyle
  +
\int_0^\lambda J_0(t) \, J_0(t \, x) \, t \, \mathrm d t=$ $\displaystyle
  +
\lambda \, J_1(\lambda) \, \frac{J_0(\lambda\, x)}{1\!-\!x^2}
  +
$
  +
  +
Curve $y=g(x)$ approaches the ordinate axis $x\!=\!0$ with zero derivative at point
  +
  +
$y\!=\!\lambda \, J_1(\lambda) \approx 1.2484591696955065$
  +
  +
The pole at $x\!=\!1$ due to denominator in the expression for $g(x)$ is compensated by the zero at $x\!=\!1$ in the numerator; so, the first zero of $g(x)$ is
  +
  +
$x=\,$[[BesselJZero]]$[0,2]/\lambda = \,$[[BesselJZero]]$[0,2]/$[[BesselJZero]]$[0,1]\,$ $\approx 2.2954172674276943$
  +
  +
==[[C++]] generator of curves==
  +
<poem><nomathjax><nowiki>
  +
#include<math.h>
  +
#include<stdio.h>
  +
#include<stdlib.h>
  +
#define DB double
  +
#define DO(x,y) for(x=0;x<y;x++)
  +
DB jnp(int n,DB x){ return .5*( jn(n-1,x)-jn(n+1,x) ) ; } // Derivative of n th Bessel
  +
DB jnz(int v, int k){ DB x,t; t=M_PI*(k+.5*v-.25); x= t - (v*v-.25)*.5/t;
  +
x-= jn(v,x)/jnp(v,x); // Newton adjustment of the root
  +
x-= jn(v,x)/jnp(v,x);
  +
x-= jn(v,x)/jnp(v,x);
  +
return x; } // the k th zero of v th Bessel
  +
#include"ado.cin"
  +
  +
int main(){ int m,n,v,k; DB s, x,y; FILE *o;
  +
int M=8;
  +
DB X[M+1],W[M+1],T[M+1][M+1],TT[M+1][M+1], F[M],G[M];
  +
DB S=jnz(0,M+1);
  +
DB jd=jnz(0,1);
  +
DB qs=sqrt(1./S);
  +
DB q=sqrt(2./S);
  +
for(n=1;n<M+1;n++){ x=jnz(0,n); X[n]=x*qs; y=W[n]=q/fabs(j1(x));
  +
printf("%3d %20.16lf %20.16lf\n",n,X[n],W[n]); }
  +
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]; }}
  +
//for(m=1;m<=M;m++){printf("\n");
  +
//for(n=1;n<=M;n++){printf("%14.10lf",T[m][n]); }}
  +
//printf("\n");
  +
//for(m=1;m<=M;m++){printf("\n");
  +
//for(n=1;n<=M;n++){ s=0.;
  +
// for(k=1;k<=M;k++) s+=T[m][k]*T[k][n] ;
  +
// TT[m][n]=s;printf("%14.10lf",TT[m][n]);
  +
// }}
  +
//printf("\n\n");
  +
//for(m=1;m<=M;m++){x=X[m]; y=exp(-x*x/2.); F[m]=y*W[m]; /* printf("%14.10lf",F[m]); */ }
  +
for(m=1;m<=M;m++){x=X[m]; if(x<jd) y=j0(x); else y=0.; F[m]=y*W[m]; /*printf("%14.10lf",F[m]);*/ }
  +
/*printf("\n");*/
  +
for(m=1;m<=M;m++){ s=0.; for(n=1;n<=M;n++) s+=T[m][n]*F[n];
  +
G[m]=s;}
  +
//for(m=1;m<=M;m++) printf("%14.10lf",G[m]);
  +
//printf("\n");
  +
o=fopen("test8trube.eps", "w"); ado(o,520,220);
  +
#define M(x,y) fprintf(o,"%8.4lf %8.4lf M\n",0.+x, 0.+y);
  +
#define L(x,y) fprintf(o,"%8.4lf %8.4lf L\n",0.+x, 0.+y);
  +
#define O(x,y) fprintf(o,"%8.4lf %8.4lf O\n",0.+x, 0.+y);
  +
#define o(x,y) fprintf(o,"%8.4lf %8.4lf o\n",0.+x, 0.+y);
  +
fprintf(o,"10 60 translate 100 100 scale\n");
  +
DO(n,6){M(n,1)L(n,0)}
  +
DO(n,3){M(0,n/2.)L(5,n/2.)}
  +
fprintf(o,"2 setlinecap 2 setlinejoin .006 W S\n");
  +
  +
//M(0,1) DO(n,500){x=.02*(n+1); y=exp(-x*x/2.); L(x,y)}
  +
M(0,1) DO(n,127){x=.04*(n+1); y=j0(x); L(x,y)}
  +
fprintf(o,"1 setlinecap 1 setlinejoin S\n");
  +
  +
DO(n,127){x=.04*(n+.2); y=j0(jd*x)/(1.-x*x)*jd*j1(jd); if(n==0)M(x,y) else L(x,y)
  +
printf("%6.4lf %6.3lf\n",x,y);}
  +
fprintf(o,"1 setlinecap 1 setlinejoin S\n");
  +
  +
fprintf(o,"1 0 0 RGB .016 W\n"); for(m=1;m<=M;m++) o(X[m],F[m]/W[m])
  +
fprintf(o,"0 0 1 RGB .01 W\n"); for(m=1;m<=M;m++) O(X[m],G[m]/W[m])
  +
  +
fprintf(o,"showpage\n%c%cTrailer\n",'%','%'); fclose(o);
  +
system("epstopdf test8trube.eps");
  +
system("open test8trube.pdf");
  +
return 0;
  +
}
  +
</nowiki></nomathjax></poem>
  +
  +
==[[Latex]] generator of labels==
  +
<poem><nomathjax><nowiki>
  +
\documentclass[12pt]{article}
  +
\usepackage{geometry}
  +
\usepackage{graphicx}
  +
\usepackage{rotating}
  +
\paperwidth 514pt
  +
\paperheight 188pt
  +
\topmargin -110pt
  +
\oddsidemargin -92pt
  +
\newcommand \ing {\includegraphics}
  +
\newcommand \sx {\scalebox}
  +
\newcommand \rot {\begin{rotate}}
  +
\newcommand \ero {\end{rotate}}
  +
\begin{document}
  +
\begin{picture}(410,212)
  +
%\put(2,4){\ing{bessel8testSte}}
  +
\put(2,4){\ing{test8trube}}
  +
\put(2,190){\sx{1.4}{\rot{-2}$\displaystyle y\!=\!\lambda J_1(\lambda) \frac{J_0(\lambda\, x)}{1\!-\!x^2}$ \ero}}
  +
%\put(1,189){\sx{1.4}{$\frac 3 2$}}
  +
\put(2,159){\sx{1.4}{$1$}}
  +
\put(14,150){\sx{1.4}{\rot{-6}$\displaystyle y\!=\! J_0(x)$ \ero}}
  +
\put(2,110){\sx{1.42}{$\frac 1 2$}}
  +
\put(2,59){\sx{1.4}{0}}
  +
\put(10,50){\sx{1.4}{0}}
  +
\put(108,50){\sx{1.4}{1}}
  +
\put(208,50){\sx{1.4}{2}}
  +
\put(308,50){\sx{1.4}{3}}
  +
\put(409,50){\sx{1.4}{4}}
  +
\put(504,50){\sx{1.5}{$x$}}
  +
\put(366,30){\sx{1.4}{\rot{0}$\displaystyle y\!=\! J_0(x)$ \ero}}
  +
%\put(123,114){\sx{1.4}{\rot{-52}$y\!=\! 2 J_1(2x)/x$\ero }}
  +
%\put(184,141){\sx{1.2}{step }}
  +
%\put(210,134){\sx{1.2}{\rot{-90}step\ero }}
  +
\end{picture}
  +
\end{document}
  +
</nowiki></nomathjax></poem>
  +
  +
==References==
  +
<references/>
  +
  +
[[Category:Bessel function]]
  +
[[Category:Bessel transform]]
  +
[[Category:C++]]
  +
[[Category:Discrete Bessel]]
  +
[[Category:Explicit plot]]
  +
[[Category:Latex]]
  +
[[Category:Test]]

Latest revision as of 08:53, 1 December 2018

Test of the Discrete Bessel transform of the truncated Bessel function

$f(x)=J_0(x)\,$UnitStep$(\lambda\!-\!x)$

where $~ \lambda=\,$BesselJZero$[0,1]\approx 2.404825557695773$

Description

Small red circle represent the input function, represented at the grid of $M\!=\!8$ points; the coordinates of these circles are $(x_n,f(x_n))$ for $n=1..8$, where $x_n$ are abscissas of the grid for the discrete Bessel.

These circles follow the smooth profile $y\!=\!J_0(x)$, shown with thin black curve

The input array is determined with $f_n\!=\! f(x_n)\, w_n$ for $n=1.8$

The transform matrix $T$ of the Discrete Bessel determines the output $\displaystyle g_n=\sum_{m=1}^M T_{n,m}\, f_m$

This output is represented by the big blue circles wit coordinates $(x_n,g_n/w_n)$ for $n=1.8$

It follows the expected smooth curve

$\displaystyle y=g(x)=$ $\displaystyle \int_0^\infty f(t) \, J_0(t \, x) \, t \, \mathrm d t=$ $\displaystyle \int_0^\lambda J_0(t) \, J_0(t \, x) \, t \, \mathrm d t=$ $\displaystyle \lambda \, J_1(\lambda) \, \frac{J_0(\lambda\, x)}{1\!-\!x^2} $

Curve $y=g(x)$ approaches the ordinate axis $x\!=\!0$ with zero derivative at point

$y\!=\!\lambda \, J_1(\lambda) \approx 1.2484591696955065$

The pole at $x\!=\!1$ due to denominator in the expression for $g(x)$ is compensated by the zero at $x\!=\!1$ in the numerator; so, the first zero of $g(x)$ is

$x=\,$BesselJZero$[0,2]/\lambda = \,$BesselJZero$[0,2]/$BesselJZero$[0,1]\,$ $\approx 2.2954172674276943$

C++ generator of curves


#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#define DB double
#define DO(x,y) for(x=0;x<y;x++)
DB jnp(int n,DB x){ return .5*( jn(n-1,x)-jn(n+1,x) ) ; } // Derivative of n th Bessel
DB jnz(int v, int k){ DB x,t; t=M_PI*(k+.5*v-.25); x= t - (v*v-.25)*.5/t;
                x-= jn(v,x)/jnp(v,x); // Newton adjustment of the root
                x-= jn(v,x)/jnp(v,x);
                x-= jn(v,x)/jnp(v,x);
                return x; } // the k th zero of v th Bessel
#include"ado.cin"

int main(){ int m,n,v,k; DB s, x,y; FILE *o;
int M=8;
DB X[M+1],W[M+1],T[M+1][M+1],TT[M+1][M+1], F[M],G[M];
DB S=jnz(0,M+1);
DB jd=jnz(0,1);
DB qs=sqrt(1./S);
DB q=sqrt(2./S);
for(n=1;n<M+1;n++){ x=jnz(0,n); X[n]=x*qs; y=W[n]=q/fabs(j1(x));
                        printf("%3d %20.16lf %20.16lf\n",n,X[n],W[n]); }
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]; }}
//for(m=1;m<=M;m++){printf("\n");
//for(n=1;n<=M;n++){printf("%14.10lf",T[m][n]); }}
//printf("\n");
//for(m=1;m<=M;m++){printf("\n");
//for(n=1;n<=M;n++){ s=0.;
// for(k=1;k<=M;k++) s+=T[m][k]*T[k][n] ;
// TT[m][n]=s;printf("%14.10lf",TT[m][n]);
// }}
//printf("\n\n");
//for(m=1;m<=M;m++){x=X[m]; y=exp(-x*x/2.); F[m]=y*W[m]; /* printf("%14.10lf",F[m]); */ }
for(m=1;m<=M;m++){x=X[m]; if(x<jd) y=j0(x); else y=0.; F[m]=y*W[m]; /*printf("%14.10lf",F[m]);*/ }
/*printf("\n");*/
for(m=1;m<=M;m++){ s=0.; for(n=1;n<=M;n++) s+=T[m][n]*F[n];
                 G[m]=s;}
//for(m=1;m<=M;m++) printf("%14.10lf",G[m]);
//printf("\n");
o=fopen("test8trube.eps", "w"); ado(o,520,220);
#define M(x,y) fprintf(o,"%8.4lf %8.4lf M\n",0.+x, 0.+y);
#define L(x,y) fprintf(o,"%8.4lf %8.4lf L\n",0.+x, 0.+y);
#define O(x,y) fprintf(o,"%8.4lf %8.4lf O\n",0.+x, 0.+y);
#define o(x,y) fprintf(o,"%8.4lf %8.4lf o\n",0.+x, 0.+y);
fprintf(o,"10 60 translate 100 100 scale\n");
DO(n,6){M(n,1)L(n,0)}
DO(n,3){M(0,n/2.)L(5,n/2.)}
fprintf(o,"2 setlinecap 2 setlinejoin .006 W S\n");

//M(0,1) DO(n,500){x=.02*(n+1); y=exp(-x*x/2.); L(x,y)}
M(0,1) DO(n,127){x=.04*(n+1); y=j0(x); L(x,y)}
fprintf(o,"1 setlinecap 1 setlinejoin S\n");

DO(n,127){x=.04*(n+.2); y=j0(jd*x)/(1.-x*x)*jd*j1(jd); if(n==0)M(x,y) else L(x,y)
                        printf("%6.4lf %6.3lf\n",x,y);}
fprintf(o,"1 setlinecap 1 setlinejoin S\n");

fprintf(o,"1 0 0 RGB .016 W\n"); for(m=1;m<=M;m++) o(X[m],F[m]/W[m])
fprintf(o,"0 0 1 RGB .01 W\n"); for(m=1;m<=M;m++) O(X[m],G[m]/W[m])

fprintf(o,"showpage\n%c%cTrailer\n",'%','%'); fclose(o);
system("epstopdf test8trube.eps");
system("open test8trube.pdf");
return 0;
}

Latex generator of labels


\documentclass[12pt]{article}
\usepackage{geometry}
\usepackage{graphicx}
\usepackage{rotating}
\paperwidth 514pt
\paperheight 188pt
\topmargin -110pt
\oddsidemargin -92pt
\newcommand \ing {\includegraphics}
\newcommand \sx {\scalebox}
\newcommand \rot {\begin{rotate}}
\newcommand \ero {\end{rotate}}
\begin{document}
\begin{picture}(410,212)
%\put(2,4){\ing{bessel8testSte}}
\put(2,4){\ing{test8trube}}
\put(2,190){\sx{1.4}{\rot{-2}$\displaystyle y\!=\!\lambda J_1(\lambda) \frac{J_0(\lambda\, x)}{1\!-\!x^2}$ \ero}}
%\put(1,189){\sx{1.4}{$\frac 3 2$}}
\put(2,159){\sx{1.4}{$1$}}
\put(14,150){\sx{1.4}{\rot{-6}$\displaystyle y\!=\! J_0(x)$ \ero}}
\put(2,110){\sx{1.42}{$\frac 1 2$}}
\put(2,59){\sx{1.4}{0}}
\put(10,50){\sx{1.4}{0}}
\put(108,50){\sx{1.4}{1}}
\put(208,50){\sx{1.4}{2}}
\put(308,50){\sx{1.4}{3}}
\put(409,50){\sx{1.4}{4}}
\put(504,50){\sx{1.5}{$x$}}
\put(366,30){\sx{1.4}{\rot{0}$\displaystyle y\!=\! J_0(x)$ \ero}}
%\put(123,114){\sx{1.4}{\rot{-52}$y\!=\! 2 J_1(2x)/x$\ero }}
%\put(184,141){\sx{1.2}{step }}
%\put(210,134){\sx{1.2}{\rot{-90}step\ero }}
\end{picture}
\end{document}

References

File history

Click on a date/time to view the file as it appeared at that time.

Date/TimeThumbnailDimensionsUserComment
current06:14, 1 December 2018Thumbnail for version as of 06:14, 1 December 20182,133 × 780 (135 KB)Maintenance script (talk | contribs)Importing image file
  • You cannot overwrite this file.

The following page links to this file:

Metadata