Difference between revisions of "File:Koritestplot.jpg"

From TORI
Jump to navigation Jump to search
(Importing image file)
 
 
Line 1: Line 1:
  +
Test of implementation of function [[kori]] through the built-in C++ function j0.
−
Importing image file
  
  +
  +
for various values of the argument $x$, function [[kori]] can be expressed with
  +
  +
[[kori]]$(x)=\displaystyle \frac{J_0\big(L_1 \sqrt{x}\big)}{1-x}$
  +
  +
where $L_1$ is first zero of the Bessel function function $j_0$. This function, [[BesselJ0]], is denoted in C++ with identifier j0.
  +
The straightforward implementation is good for positive values of argument, larger than 2; for smaller, but still positive values of the argument, the representation is not good; for other values (negative and complex) this implementation fails.
  +
  +
The following special expansion is used:
  +
  +
[[kori]]$(x)=\,$ $\displaystyle
  +
\big(1\!-\!x/M_2\big)\big(1\!-\!x/M_3\big) \sum_{n=0}^{\infty}$ $ c_n x^n$
  +
  +
where $M_m=(L_m/L_1)^2$, and $L_n$ is $n$th zero of function [[BesselJ0]].
  +
The coefficients $c$ are evaluated using the [[Mathematica]] routine [[Series]]; approximations for these coefficients can be extracted from the C++ code below.
  +
  +
The series is truncated at the term with $n=15$, as it is sufficient for the numerical implementation for $-M_2<x\le M_2$ with at least 15 significant figures.
  +
  +
The thick red curve shows the implementation of function [[kori]] through the truncated series above.
  +
The thin black curve shows the implementation through the built-in function j0; it fully overlaps with the red curve, but can be seen at the strong zooming–in. In order to make the difference visible, deviation
  +
  +
$\mathcal D(x)= (1\!-\!x/M_2)(1\!-\!x/M_3) \sum_{n=0}^{15}$ $ c_n x^n - \, $[[j0]]$(L1*\mathrm{sqrt}(x))/(1\!-\!x)$
  +
  +
is shown with another thin black curve. In order to make it visible, it is scaled with factor $10^{15}$, as it is marked in the figure.
  +
  +
Deviation $\mathcal D(x)$ is of order of
  +
$10^{-15}$ or even $10^{-14}$, while $0\!<\!x\!<\!2$; this can be attributed to the rounding errors at evaluation of function double j0(double). For $x\!=\!1$, the representation through j0 just fails, and this is the main motivation to built-up the independent implementation of function [[kori]]. In order to keep the code transparent and easy for the reading by human, it is preferably to cover the range from $0$ to $M_2$ with single–peace approximation; such a single-peace approximations have advantages also for the numerical evaluation of the contour integrals.
  +
  +
Deviation is of order of $10^{-16}$ for $2\!<\!x\!<6$, and this verifies each of the representations mentioned.
  +
  +
Deviation grows–up at $x\!>\!6$, but still can be used to plot graphics at least until the next zero of function [[kori]], which is $M_2$.
  +
  +
The comparison shows, that in the range $0<x<M_1\approx 5.3$, the implementation [[kori0]] with 15 terms,
  +
  +
[[kori0]]$(x)= (1\!-\!x/M_2)(1\!-\!x/M_3) \sum_{n=0}^{15} c_n x^n$
  +
  +
is better than the straightforward implementation through the built-in function j0.
  +
  +
Approximation [[kori0]] is expected to be used for evaluation of integrals expressed through function [[kori]].
  +
  +
==[[C++]] generator of curves==
  +
  +
<poem><nomathjax><nowiki>
  +
#include <stdio.h>
  +
#include <stdlib.h>
  +
#define DB double
  +
#define DO(x,y) for(x=0;x<y;x++)
  +
//using namespace std;
  +
#include <complex>
  +
typedef std::complex<double> z_type;
  +
#define Re(x) x.real()
  +
#define Im(x) x.imag()
  +
#define I z_type(0.,1.)
  +
#define M(x,y) fprintf(o,"%6.4lf %6.4lf M\n",0.+x,0.+y);
  +
#define L(x,y) fprintf(o,"%6.4lf %6.4lf L\n",0.+x,0.+y);
  +
  +
// #include "korifit76.cin"
  +
#include "ado.cin"
  +
  +
long double L= 2.4048255576957727686216318793265 ;
  +
  +
DB LN[9]={ 0., 2.404825557695773, 5.5200781102863115,
  +
8.653727912911013,11.791534439014281,14.930917708487787,
  +
18.071063967910924,21.21163662987926, 24.352471530749302};
  +
  +
long double kori0(long double x){ int n,m; long double s;
  +
long double c[16]={1,
  +
-0.256005011219941993570670508019,
  +
0.0281978125768042007057113441199,
  +
-0.00181258393780241662780680613867,
  +
0.0000775822332475845672943781286224,
  +
-2.38609544935449493722623606533e-6,
  +
5.5465254159392191445565645766e-8,
  +
-1.01054597146460752951273870914e-9,
  +
1.48348772753372628782330712313e-11,
  +
-1.79332434839872451210764178471e-13,
  +
1.81697031752603393402790223753e-15,
  +
-1.56569115453044785387543614077e-17,
  +
1.16170076283001586118141461154e-19,
  +
-7.50065787627627390138399503073e-22,
  +
4.25298321710799220467323756624e-24,
  +
-2.13477668517290665487301104539e-26};
  +
m=15;
  +
s=c[m]*x; for(n=m-1;n>0;n--) {s+=c[n]; s*=x;}
  +
return (1+s)*(1-(long double)0.189791479516754136723328431595*x); }
  +
  +
/*
  +
double korid(double x){ int n,m; double s;
  +
double c[16]={1,
  +
-0.256005011219941993570670508019,
  +
0.0281978125768042007057113441199,
  +
-0.00181258393780241662780680613867,
  +
0.0000775822332475845672943781286224,
  +
-2.38609544935449493722623606533e-6,
  +
5.5465254159392191445565645766e-8,
  +
-1.01054597146460752951273870914e-9,
  +
1.48348772753372628782330712313e-11,
  +
-1.79332434839872451210764178471e-13,
  +
1.81697031752603393402790223753e-15,
  +
-1.56569115453044785387543614077e-17,
  +
1.16170076283001586118141461154e-19,
  +
-7.50065787627627390138399503073e-22,
  +
4.25298321710799220467323756624e-24,
  +
-2.13477668517290665487301104539e-26};
  +
m=15;
  +
s=c[m]*x; for(n=m-1;n>0;n--) {s+=c[n]; s*=x;}
  +
return (1+s)*(1-0.189791479516754136723328431595*x); }
  +
*/
  +
  +
int main(){int m,n,i; double x,y,z,t,u; FILE *o;
  +
long double X,Y,Z,T,U;
  +
  +
//o=fopen("61.eps","w");
  +
o=fopen("koritestplo.eps","w");
  +
ado(o,1040,500);
  +
fprintf(o,"20 200 translate 100 100 scale 2 setlinecap 1 setlinejoin\n");
  +
  +
for(m=0;m<11;m+=1) {M(m,0)L(m,2)}
  +
for(n=0;n<3;n+=1){M(0,n)L(10,n)}
  +
fprintf(o,"0.007 W S\n");
  +
  +
for(m=2;m<4;m++){ x=LN[m]/LN[1]; x*=x; M(x,0)L(x,1) }
  +
fprintf(o,"0.004 W 1 0 1 RGB S\n");
  +
  +
for(n=0;n<1060;n+=10){ x=-.3+.01*n; y=kori0(x); if(n==0) M(x,y) else L(x,y);} fprintf(o,"1 0 0 RGB .03 W S\n");
  +
  +
for(n=0;n<1040;n+=10){ x=.01*(n+.5); y=j0(L*sqrt(x))/(1.-x); if(n==0) M(x,y) else L(x,y);} fprintf(o,"0 0 0 RGB .01 W S\n");
  +
  +
//for(n=0;n<720;n+=1){ x=.01*(n+.5); y=korid(x)-j0(L*sqrt(x))/(1.-x); y*=1.e15; if(n==0) M(x,y) else L(x,y);}
  +
//fprintf(o,"0 0 1 RGB .0009 W S\n");
  +
  +
/*
  +
for(n=0;n<720;n+=1){ X=(long double).01*(n+(long double).5); x=(double)X;
  +
Y=kori0(X);
  +
Z=korid(x);
  +
T=Y-Z;
  +
U=T*1.e15; u=(double)U; y=u;if(n==0) M(x,y) else L(x,y);}
  +
fprintf(o,"0 0 1 RGB .0008 W S\n");
  +
*/
  +
  +
for(n=0;n<730;n+=1){ X=(long double).01*(n+(long double).5); x=(double)X;
  +
Y=kori0(X);
  +
Z=(long double)j0(L*sqrt(x))/(1.-x);
  +
T=Y-Z;
  +
U=T*1.e15; u=(double)U; y=u;if(n==0) M(x,y) else L(x,y);}
  +
fprintf(o,"0 0 0 RGB .0008 W S\n");
  +
  +
fprintf(o,"showpage\n%c%cTrailer\n",'%','%'); fclose(o);
  +
system("epstopdf koritestplo.eps");
  +
system( "open koritestplo.pdf");
  +
return 0;}
  +
  +
</nowiki></nomathjax></poem>
  +
==[[Latex]] generator of labels==
  +
<poem><nomathjax><nowiki>
  +
\documentclass[12pt]{article}
  +
\usepackage{geometry}
  +
\usepackage{graphics}
  +
\usepackage{rotating}
  +
%\paperwidth 570pt
  +
\paperwidth 780pt
  +
%\paperheight 138pt
  +
\paperheight 480pt
  +
\textwidth 420pt
  +
\textheight 600pt
  +
\topmargin -108pt
  +
\oddsidemargin -72pt
  +
\newcommand \ds {\displaystyle}
  +
\newcommand \sx {\scalebox}
  +
\newcommand \rme {\mathrm{e}}
  +
\newcommand \rot {\begin{rotate}}
  +
\newcommand \ero {\end{rotate}}
  +
\newcommand \ing {\includegraphics}
  +
\pagestyle{empty}
  +
\parindent 0pt
  +
\begin{document}
  +
\begin{picture}(640,500)
  +
%\put(0,20){\ing{koriplo}}
  +
\put(0,20){\ing{koritestplo}}
  +
%\put(384,468){\sx{2.8}{$y\!=\!\mathrm{kori}(x)$}}
  +
\put(0,480){\sx{2.8}{$y$}}
  +
\put(0,410){\sx{2.8}{$2$}}
  +
\put(0,310){\sx{2.8}{$1$}}
  +
\put(0,210){\sx{2.8}{$0$}}
  +
%\put(96,0){\sx{2.8}{$0$}}
  +
\put( 14,190){\sx{2.8}{$0$}}
  +
\put(114,190){\sx{2.8}{$1$}}
  +
\put(214,190){\sx{2.8}{$2$}}
  +
\put(314,190){\sx{2.8}{$3$}}
  +
\put(414,190){\sx{2.8}{$4$}}
  +
\put(514,190){\sx{2.8}{$5$}}
  +
\put(536,180){\sx{2.7}{$\frac{L_2^2}{L_1^2}$}}
  +
\put(614,190){\sx{2.8}{$6$}}
  +
\put(714,190){\sx{2.8}{$7$}}
  +
\put(762,190){\sx{2.8}{$x$}}
  +
%\put(806,0){\sx{2.8}{$8$}}
  +
%\put(906,0){\sx{2.8}{$9$}}
  +
%\put(1000,0){\sx{2.9}{$x$}}
  +
\put(170,278){\rot{-12}\sx{2.8}{$y\!=\! \mathrm{kori}(x)$}\ero}
  +
\put(598,468){\rot{0}\sx{2.8}{$y\!=\! 10^{15}\,\mathcal D(x)$}\ero}
  +
%\put(656,148){\rot{65}\sx{2.8}{$y\!=\! 10^{15}\,\mathcal D(x)$}\ero}
  +
\end{picture}\end{document}
  +
</nowiki></nomathjax></poem>
  +
  +
==References==
  +
<references/>
  +
  +
[[Category:Bessel function]]
  +
[[Category:BesselJ0]]
  +
[[Category:Explicit plot]]
  +
[[Category:Kori]]
  +
[[Category:Morinaga function]]

Latest revision as of 08:40, 1 December 2018

Test of implementation of function kori through the built-in C++ function j0.

for various values of the argument $x$, function kori can be expressed with

kori$(x)=\displaystyle \frac{J_0\big(L_1 \sqrt{x}\big)}{1-x}$

where $L_1$ is first zero of the Bessel function function $j_0$. This function, BesselJ0, is denoted in C++ with identifier j0. The straightforward implementation is good for positive values of argument, larger than 2; for smaller, but still positive values of the argument, the representation is not good; for other values (negative and complex) this implementation fails.

The following special expansion is used:

kori$(x)=\,$ $\displaystyle \big(1\!-\!x/M_2\big)\big(1\!-\!x/M_3\big) \sum_{n=0}^{\infty}$ $ c_n x^n$

where $M_m=(L_m/L_1)^2$, and $L_n$ is $n$th zero of function BesselJ0. The coefficients $c$ are evaluated using the Mathematica routine Series; approximations for these coefficients can be extracted from the C++ code below.

The series is truncated at the term with $n=15$, as it is sufficient for the numerical implementation for $-M_2<x\le M_2$ with at least 15 significant figures.

The thick red curve shows the implementation of function kori through the truncated series above. The thin black curve shows the implementation through the built-in function j0; it fully overlaps with the red curve, but can be seen at the strong zooming–in. In order to make the difference visible, deviation

$\mathcal D(x)= (1\!-\!x/M_2)(1\!-\!x/M_3) \sum_{n=0}^{15}$ $ c_n x^n - \, $j0$(L1*\mathrm{sqrt}(x))/(1\!-\!x)$

is shown with another thin black curve. In order to make it visible, it is scaled with factor $10^{15}$, as it is marked in the figure.

Deviation $\mathcal D(x)$ is of order of $10^{-15}$ or even $10^{-14}$, while $0\!<\!x\!<\!2$; this can be attributed to the rounding errors at evaluation of function double j0(double). For $x\!=\!1$, the representation through j0 just fails, and this is the main motivation to built-up the independent implementation of function kori. In order to keep the code transparent and easy for the reading by human, it is preferably to cover the range from $0$ to $M_2$ with single–peace approximation; such a single-peace approximations have advantages also for the numerical evaluation of the contour integrals.

Deviation is of order of $10^{-16}$ for $2\!<\!x\!<6$, and this verifies each of the representations mentioned.

Deviation grows–up at $x\!>\!6$, but still can be used to plot graphics at least until the next zero of function kori, which is $M_2$.

The comparison shows, that in the range $0<x<M_1\approx 5.3$, the implementation kori0 with 15 terms,

kori0$(x)= (1\!-\!x/M_2)(1\!-\!x/M_3) \sum_{n=0}^{15} c_n x^n$

is better than the straightforward implementation through the built-in function j0.

Approximation kori0 is expected to be used for evaluation of integrals expressed through function kori.

C++ generator of curves


#include <stdio.h>
#include <stdlib.h>
#define DB double
#define DO(x,y) for(x=0;x<y;x++)
//using namespace std;
#include <complex>
typedef std::complex<double> z_type;
#define Re(x) x.real()
#define Im(x) x.imag()
#define I z_type(0.,1.)
#define M(x,y) fprintf(o,"%6.4lf %6.4lf M\n",0.+x,0.+y);
#define L(x,y) fprintf(o,"%6.4lf %6.4lf L\n",0.+x,0.+y);

// #include "korifit76.cin"
#include "ado.cin"

long double L= 2.4048255576957727686216318793265 ;

DB LN[9]={ 0., 2.404825557695773, 5.5200781102863115,
 8.653727912911013,11.791534439014281,14.930917708487787,
18.071063967910924,21.21163662987926, 24.352471530749302};

long double kori0(long double x){ int n,m; long double s;
long double c[16]={1,
-0.256005011219941993570670508019,
 0.0281978125768042007057113441199,
-0.00181258393780241662780680613867,
 0.0000775822332475845672943781286224,
-2.38609544935449493722623606533e-6,
 5.5465254159392191445565645766e-8,
-1.01054597146460752951273870914e-9,
 1.48348772753372628782330712313e-11,
-1.79332434839872451210764178471e-13,
 1.81697031752603393402790223753e-15,
-1.56569115453044785387543614077e-17,
 1.16170076283001586118141461154e-19,
-7.50065787627627390138399503073e-22,
 4.25298321710799220467323756624e-24,
-2.13477668517290665487301104539e-26};
m=15;
s=c[m]*x; for(n=m-1;n>0;n--) {s+=c[n]; s*=x;}
return (1+s)*(1-(long double)0.189791479516754136723328431595*x); }

/*
double korid(double x){ int n,m; double s;
double c[16]={1,
-0.256005011219941993570670508019,
 0.0281978125768042007057113441199,
-0.00181258393780241662780680613867,
 0.0000775822332475845672943781286224,
-2.38609544935449493722623606533e-6,
 5.5465254159392191445565645766e-8,
-1.01054597146460752951273870914e-9,
 1.48348772753372628782330712313e-11,
-1.79332434839872451210764178471e-13,
 1.81697031752603393402790223753e-15,
-1.56569115453044785387543614077e-17,
 1.16170076283001586118141461154e-19,
-7.50065787627627390138399503073e-22,
 4.25298321710799220467323756624e-24,
-2.13477668517290665487301104539e-26};
m=15;
s=c[m]*x; for(n=m-1;n>0;n--) {s+=c[n]; s*=x;}
return (1+s)*(1-0.189791479516754136723328431595*x); }
*/

int main(){int m,n,i; double x,y,z,t,u; FILE *o;
long double X,Y,Z,T,U;

//o=fopen("61.eps","w");
o=fopen("koritestplo.eps","w");
ado(o,1040,500);
fprintf(o,"20 200 translate 100 100 scale 2 setlinecap 1 setlinejoin\n");

for(m=0;m<11;m+=1) {M(m,0)L(m,2)}
for(n=0;n<3;n+=1){M(0,n)L(10,n)}
fprintf(o,"0.007 W S\n");

for(m=2;m<4;m++){ x=LN[m]/LN[1]; x*=x; M(x,0)L(x,1) }
fprintf(o,"0.004 W 1 0 1 RGB S\n");

for(n=0;n<1060;n+=10){ x=-.3+.01*n; y=kori0(x); if(n==0) M(x,y) else L(x,y);} fprintf(o,"1 0 0 RGB .03 W S\n");

for(n=0;n<1040;n+=10){ x=.01*(n+.5); y=j0(L*sqrt(x))/(1.-x); if(n==0) M(x,y) else L(x,y);} fprintf(o,"0 0 0 RGB .01 W S\n");

//for(n=0;n<720;n+=1){ x=.01*(n+.5); y=korid(x)-j0(L*sqrt(x))/(1.-x); y*=1.e15; if(n==0) M(x,y) else L(x,y);}
//fprintf(o,"0 0 1 RGB .0009 W S\n");

/*
for(n=0;n<720;n+=1){ X=(long double).01*(n+(long double).5); x=(double)X;
                        Y=kori0(X);
                        Z=korid(x);
                        T=Y-Z;
                        U=T*1.e15; u=(double)U; y=u;if(n==0) M(x,y) else L(x,y);}
fprintf(o,"0 0 1 RGB .0008 W S\n");
*/

for(n=0;n<730;n+=1){ X=(long double).01*(n+(long double).5); x=(double)X;
                        Y=kori0(X);
                        Z=(long double)j0(L*sqrt(x))/(1.-x);
                        T=Y-Z;
                        U=T*1.e15; u=(double)U; y=u;if(n==0) M(x,y) else L(x,y);}
fprintf(o,"0 0 0 RGB .0008 W S\n");

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

Latex generator of labels


\documentclass[12pt]{article}
\usepackage{geometry}
\usepackage{graphics}
\usepackage{rotating}
%\paperwidth 570pt
\paperwidth 780pt
%\paperheight 138pt
\paperheight 480pt
\textwidth 420pt
\textheight 600pt
\topmargin -108pt
\oddsidemargin -72pt
\newcommand \ds {\displaystyle}
\newcommand \sx {\scalebox}
\newcommand \rme {\mathrm{e}}
\newcommand \rot {\begin{rotate}}
\newcommand \ero {\end{rotate}}
\newcommand \ing {\includegraphics}
\pagestyle{empty}
\parindent 0pt
\begin{document}
\begin{picture}(640,500)
%\put(0,20){\ing{koriplo}}
\put(0,20){\ing{koritestplo}}
%\put(384,468){\sx{2.8}{$y\!=\!\mathrm{kori}(x)$}}
\put(0,480){\sx{2.8}{$y$}}
\put(0,410){\sx{2.8}{$2$}}
\put(0,310){\sx{2.8}{$1$}}
\put(0,210){\sx{2.8}{$0$}}
%\put(96,0){\sx{2.8}{$0$}}
\put( 14,190){\sx{2.8}{$0$}}
\put(114,190){\sx{2.8}{$1$}}
\put(214,190){\sx{2.8}{$2$}}
\put(314,190){\sx{2.8}{$3$}}
\put(414,190){\sx{2.8}{$4$}}
\put(514,190){\sx{2.8}{$5$}}
\put(536,180){\sx{2.7}{$\frac{L_2^2}{L_1^2}$}}
\put(614,190){\sx{2.8}{$6$}}
\put(714,190){\sx{2.8}{$7$}}
\put(762,190){\sx{2.8}{$x$}}
%\put(806,0){\sx{2.8}{$8$}}
%\put(906,0){\sx{2.8}{$9$}}
%\put(1000,0){\sx{2.9}{$x$}}
\put(170,278){\rot{-12}\sx{2.8}{$y\!=\! \mathrm{kori}(x)$}\ero}
\put(598,468){\rot{0}\sx{2.8}{$y\!=\! 10^{15}\,\mathcal D(x)$}\ero}
%\put(656,148){\rot{65}\sx{2.8}{$y\!=\! 10^{15}\,\mathcal D(x)$}\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:12, 1 December 2018Thumbnail for version as of 06:12, 1 December 20181,618 × 996 (144 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:Koritestplot.jpg&oldid=19728"