Difference between revisions of "File:Test8trubes300.jpg"

From TORI
Jump to navigation Jump to 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 uses this file:

Metadata

This file contains additional information, probably added from the digital camera or scanner used to create or digitize it.

If the file has been modified from its original state, some details may not fully reflect the modified file.

Retrieved from "https://mizugadro.mydns.jp/t/index.php?title=File:Test8trubes300.jpg&oldid=21844"