QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

gsl_sf_bessel_jl_array and gsl_sf_bessel_jl_steed_array are inprecise. 

Must do a different implementation.


Former-commit-id: 37da4cf29b5d8706c639a9a3af61653d90ccc37e
This commit is contained in:
Marek Nečada 2018-03-14 14:10:11 +02:00
parent 3b407ab40a
commit fb9bbff4b8
3 changed files with 57 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
notes/tftests.lyx
View File

@ -343,7 +343,7 @@ key "xu_efficient_1998"
\begin_inset Formula $B_{mn\mu\nu}$
\end_inset
se liší v několika ohledech:
se lišejí v několika ohledech:
\end_layout
\begin_layout Enumerate
@ -389,9 +389,18 @@ Ve starší práci je poslední člen sumy
\begin_inset Formula $\left(m,n,\mu,\nu\right)=\left(0,1,-1,1\right)$
\end_inset
!!! (Opraveno, použit Cruzanův tvar v
!!! (Při numerickém srovnání Xuových vzorců
\begin_inset CommandInset citation
LatexCommand cite
after "(77–80)"
key "xu_calculation_1996"
\end_inset
ve staré práci a cruzanovských vzorců
\begin_inset CommandInset citation
LatexCommand cite
after "(59–61, ...)"
key "xu_efficient_1998"
\end_inset
@ -430,7 +439,9 @@ Legendreovy polynomy,
\end_layout
\begin_layout Itemize
Besselovy funkce.
Besselovy funkce – nepřesné jak sviňa zejména u derivací besselových funkcí
prvního druhu.
Nutno zvolit jinou implementaci.
\end_layout
\end_deeper

41
tests/bessel_precision/besseltest.c Normal file
View File

@ -0,0 +1,41 @@
#include <gsl/gsl_sf_bessel.h>
#include <stdio.h>
#include <math.h>
int main() {
int lMax;
while (1 == scanf("%d", &lMax)) {
double x;
if (1 != scanf("%lf", &x))
abort();
double gsl[lMax+2], relerr[lMax+1], orig[lMax+1];
for (int l = 0; l <= lMax; l++)
if (1 != scanf("%lf", orig+l))
abort();
#if defined JTEST || defined DJTEST
if (gsl_sf_bessel_jl_array(lMax+1, x, gsl))
#elif defined YTEST || defined DYTEST
if (gsl_sf_bessel_yl_array(lMax+1, x, gsl))
#else
if (gsl_sf_bessel_jl_steed_array(lMax+1, x, gsl))
#endif
abort();
#if defined DJTEST || defined DYTEST || defined DJTEST_STEED
for (int l = 0; l <= lMax; l++)
gsl[l] = -gsl[l+1] + (l/x) * gsl[l];
#endif
printf("x = %.16g, lMax = %d:\nsage: ", x, lMax);
for (int l = 0; l <= lMax; l++)
printf("%.16g ", orig[l]);
printf("\ngsl: ");
for (int l = 0; l <= lMax; l++)
printf("%.16g ", gsl[l]);
printf("\nrell: ");
for (int l = 0; l <= lMax; l++)
printf("%.16g ", 2 * fabs(gsl[l] - orig[l]) / (fabs(gsl[l]) + fabs(gsl[l])));
putchar('\n');
}
return 0;
}

3
tests/bessel_precision/compileall.sh
View File

@ -2,7 +2,8 @@
c99 -ggdb -o jtest -DJTEST besseltest.c -lgsl -lblas -lm
c99 -ggdb -o ytest -DYTEST besseltest.c -lgsl -lblas -lm
c99 -ggdb -o jtest_steed -DJTEST_STEED besseltest.c -lgsl -lblas -lm
c99 -ggdb -o djtest -DDJTEST besseltest.c -lgsl -lblas -lm
c99 -ggdb -o dytest -DDYTEST besseltest.c -lgsl -lblas -lm
c99 -ggdb -o djtest_steed -DJTEST_STEED besseltest.c -lgsl -lblas -lm
c99 -ggdb -o djtest_steed -DDJTEST_STEED besseltest.c -lgsl -lblas -lm