Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/besselj.cc @ 8278:ab0674a8b345

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fix scaling factor for negative alpha in zbesi,cbesi add bessel function tests add all scaling factors to bessel documentation
author Brian Gough <bjg@gnu.org>
date Tue, 28 Oct 2008 10:47:26 -0400
parents 9d080df0c843
children eb63fbe60fab
comparison
equal deleted inserted replaced
8277:b4a1a5472191 8278:ab0674a8b345
386 @deftypefnx {Loadable Function} {[@var{h}, @var{ierr}] =} besselh (@var{alpha}, @var{k}, @var{x}, @var{opt})\n\ 386 @deftypefnx {Loadable Function} {[@var{h}, @var{ierr}] =} besselh (@var{alpha}, @var{k}, @var{x}, @var{opt})\n\
387 Compute Bessel or Hankel functions of various kinds:\n\ 387 Compute Bessel or Hankel functions of various kinds:\n\
388 \n\ 388 \n\
389 @table @code\n\ 389 @table @code\n\
390 @item besselj\n\ 390 @item besselj\n\
391 Bessel functions of the first kind.\n\ 391 Bessel functions of the first kind. If the argument @var{opt} is supplied, \n\
392 the result is multiplied by @code{exp(-abs(imag(x)))}.\n\
392 @item bessely\n\ 393 @item bessely\n\
393 Bessel functions of the second kind.\n\ 394 Bessel functions of the second kind. If the argument @var{opt} is supplied,\n\
395 the result is multiplied by @code{exp(-abs(imag(x)))}.\n\
394 @item besseli\n\ 396 @item besseli\n\
395 Modified Bessel functions of the first kind.\n\ 397 Modified Bessel functions of the first kind. If the argument @var{opt} is supplied,\n\
398 the result is multiplied by @code{exp(-abs(real(x)))}.\n\
396 @item besselk\n\ 399 @item besselk\n\
397 Modified Bessel functions of the second kind.\n\ 400 Modified Bessel functions of the second kind. If the argument @var{opt} is supplied,\n\
401 the result is multiplied by @code{exp(x)}.\n\
398 @item besselh\n\ 402 @item besselh\n\
399 Compute Hankel functions of the first (@var{k} = 1) or second (@var{k}\n\ 403 Compute Hankel functions of the first (@var{k} = 1) or second (@var{k}\n\
400 = 2) kind.\n\ 404 = 2) kind. If the argument @var{opt} is supplied, the result is multiplied by\n\
405 @code{exp (-I*@var{x})} for @var{k} = 1 or @code{exp (I*@var{x})} for\n\
406 @var{k} = 2.\n\
401 @end table\n\ 407 @end table\n\
402 \n\
403 If the argument @var{opt} is supplied, the result is scaled by the\n\
404 @code{exp (-I*@var{x})} for @var{k} = 1 or @code{exp (I*@var{x})} for\n\
405 @var{k} = 2.\n\
406 \n\ 408 \n\
407 If @var{alpha} is a scalar, the result is the same size as @var{x}.\n\ 409 If @var{alpha} is a scalar, the result is the same size as @var{x}.\n\
408 If @var{x} is a scalar, the result is the same size as @var{alpha}.\n\ 410 If @var{x} is a scalar, the result is the same size as @var{alpha}.\n\
409 If @var{alpha} is a row vector and @var{x} is a column vector, the\n\ 411 If @var{alpha} is a row vector and @var{x} is a column vector, the\n\
410 result is a matrix with @code{length (@var{x})} rows and\n\ 412 result is a matrix with @code{length (@var{x})} rows and\n\
628 630
629 return retval; 631 return retval;
630 } 632 }
631 633
632 /* 634 /*
635 %! # Test values computed with GP/PARI version 2.3.3
636 %!
637 %!shared alpha, x, jx, yx, ix, kx, nix
638 %!
639 %! # Bessel functions, even order, positive and negative x
640 %! alpha = 2; x = 1.25;
641 %! jx = 0.1710911312405234823613091417;
642 %! yx = -1.193199310178553861283790424;
643 %! ix = 0.2220184483766341752692212604;
644 %! kx = 0.9410016167388185767085460540;
645 %!
646 %!assert(besselj(alpha,x), jx, 100*eps)
647 %!assert(bessely(alpha,x), yx, 100*eps)
648 %!assert(besseli(alpha,x), ix, 100*eps)
649 %!assert(besselk(alpha,x), kx, 100*eps)
650 %!assert(besselh(alpha,1,x), jx + I*yx, 100*eps)
651 %!assert(besselh(alpha,2,x), jx - I*yx, 100*eps)
652 %!
653 %!assert(besselj(alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
654 %!assert(bessely(alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
655 %!assert(besseli(alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
656 %!assert(besselk(alpha,x,1), kx*exp(x), 100*eps)
657 %!assert(besselh(alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
658 %!assert(besselh(alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
659 %!
660 %!assert(besselj(-alpha,x), jx, 100*eps)
661 %!assert(bessely(-alpha,x), yx, 100*eps)
662 %!assert(besseli(-alpha,x), ix, 100*eps)
663 %!assert(besselk(-alpha,x), kx, 100*eps)
664 %!assert(besselh(-alpha,1,x), jx + I*yx, 100*eps)
665 %!assert(besselh(-alpha,2,x), jx - I*yx, 100*eps)
666 %!
667 %!assert(besselj(-alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
668 %!assert(bessely(-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
669 %!assert(besseli(-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
670 %!assert(besselk(-alpha,x,1), kx*exp(x), 100*eps)
671 %!assert(besselh(-alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
672 %!assert(besselh(-alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
673 %!
674 %! x *= -1;
675 %! yx = -1.193199310178553861283790424 + 0.3421822624810469647226182835*I;
676 %! kx = 0.9410016167388185767085460540 - 0.6974915263814386815610060884*I;
677 %!
678 %!assert(besselj(alpha,x), jx, 100*eps)
679 %!assert(bessely(alpha,x), yx, 100*eps)
680 %!assert(besseli(alpha,x), ix, 100*eps)
681 %!assert(besselk(alpha,x), kx, 100*eps)
682 %!assert(besselh(alpha,1,x), jx + I*yx, 100*eps)
683 %!assert(besselh(alpha,2,x), jx - I*yx, 100*eps)
684 %!
685 %!assert(besselj(alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
686 %!assert(bessely(alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
687 %!assert(besseli(alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
688 %!assert(besselk(alpha,x,1), kx*exp(x), 100*eps)
689 %!assert(besselh(alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
690 %!assert(besselh(alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
691 %!
692 %! # Bessel functions, odd order, positive and negative x
693 %! alpha = 3; x = 2.5;
694 %! jx = 0.2166003910391135247666890035;
695 %! yx = -0.7560554967536709968379029772;
696 %! ix = 0.4743704087780355895548240179;
697 %! kx = 0.2682271463934492027663765197;
698 %!
699 %!assert(besselj(alpha,x), jx, 100*eps)
700 %!assert(bessely(alpha,x), yx, 100*eps)
701 %!assert(besseli(alpha,x), ix, 100*eps)
702 %!assert(besselk(alpha,x), kx, 100*eps)
703 %!assert(besselh(alpha,1,x), jx + I*yx, 100*eps)
704 %!assert(besselh(alpha,2,x), jx - I*yx, 100*eps)
705 %!
706 %!assert(besselj(alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
707 %!assert(bessely(alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
708 %!assert(besseli(alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
709 %!assert(besselk(alpha,x,1), kx*exp(x), 100*eps)
710 %!assert(besselh(alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
711 %!assert(besselh(alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
712 %!
713 %!assert(besselj(-alpha,x), -jx, 100*eps)
714 %!assert(bessely(-alpha,x), -yx, 100*eps)
715 %!assert(besseli(-alpha,x), ix, 100*eps)
716 %!assert(besselk(-alpha,x), kx, 100*eps)
717 %!assert(besselh(-alpha,1,x), -(jx + I*yx), 100*eps)
718 %!assert(besselh(-alpha,2,x), -(jx - I*yx), 100*eps)
719 %!
720 %!assert(besselj(-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
721 %!assert(bessely(-alpha,x,1), -yx*exp(-abs(imag(x))), 100*eps)
722 %!assert(besseli(-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
723 %!assert(besselk(-alpha,x,1), kx*exp(x), 100*eps)
724 %!assert(besselh(-alpha,1,x,1), -(jx + I*yx)*exp(-I*x), 100*eps)
725 %!assert(besselh(-alpha,2,x,1), -(jx - I*yx)*exp(I*x), 100*eps)
726 %!
727 %! x *= -1;
728 %! jx = -jx;
729 %! yx = 0.7560554967536709968379029772 - 0.4332007820782270495333780070*I;
730 %! ix = -ix;
731 %! kx = -0.2682271463934492027663765197 - 1.490278591297463775542004240*I;
732 %!
733 %!assert(besselj(alpha,x), jx, 100*eps)
734 %!assert(bessely(alpha,x), yx, 100*eps)
735 %!assert(besseli(alpha,x), ix, 100*eps)
736 %!assert(besselk(alpha,x), kx, 100*eps)
737 %!assert(besselh(alpha,1,x), jx + I*yx, 100*eps)
738 %!assert(besselh(alpha,2,x), jx - I*yx, 100*eps)
739 %!
740 %!assert(besselj(alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
741 %!assert(bessely(alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
742 %!assert(besseli(alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
743 %!assert(besselk(alpha,x,1), kx*exp(x), 100*eps)
744 %!assert(besselh(alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
745 %!assert(besselh(alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
746 %!
747 %! # Bessel functions, fractional order, positive and negative x
748 %!
749 %! alpha = 3.5; x = 2.75;
750 %! jx = 0.1691636439842384154644784389;
751 %! yx = -0.8301381935499356070267953387;
752 %! ix = 0.3930540878794826310979363668;
753 %! kx = 0.2844099013460621170288192503;
754 %!
755 %!assert(besselj(alpha,x), jx, 100*eps)
756 %!assert(bessely(alpha,x), yx, 100*eps)
757 %!assert(besseli(alpha,x), ix, 100*eps)
758 %!assert(besselk(alpha,x), kx, 100*eps)
759 %!assert(besselh(alpha,1,x), jx + I*yx, 100*eps)
760 %!assert(besselh(alpha,2,x), jx - I*yx, 100*eps)
761 %!
762 %!assert(besselj(alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
763 %!assert(bessely(alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
764 %!assert(besseli(alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
765 %!assert(besselk(alpha,x,1), kx*exp(x), 100*eps)
766 %!assert(besselh(alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
767 %!assert(besselh(alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
768 %!
769 %! nix = 0.2119931212254662995364461998;
770 %!
771 %!assert(besselj(-alpha,x), yx, 100*eps)
772 %!assert(bessely(-alpha,x), -jx, 100*eps)
773 %!assert(besseli(-alpha,x), nix, 100*eps)
774 %!assert(besselk(-alpha,x), kx, 100*eps)
775 %!assert(besselh(-alpha,1,x), -I*(jx + I*yx), 100*eps)
776 %!assert(besselh(-alpha,2,x), I*(jx - I*yx), 100*eps)
777 %!
778 %!assert(besselj(-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
779 %!assert(bessely(-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
780 %!assert(besseli(-alpha,x,1), nix*exp(-abs(real(x))), 100*eps)
781 %!assert(besselk(-alpha,x,1), kx*exp(x), 100*eps)
782 %!assert(besselh(-alpha,1,x,1), -I*(jx + I*yx)*exp(-I*x), 100*eps)
783 %!assert(besselh(-alpha,2,x,1), I*(jx - I*yx)*exp(I*x), 100*eps)
784 %!
785 %! x *= -1;
786 %! jx *= -I;
787 %! yx = -0.8301381935499356070267953387*I;
788 %! ix *= -I;
789 %! kx = -0.9504059335995575096509874508*I;
790 %!
791 %!assert(besselj(alpha,x), jx, 100*eps)
792 %!assert(bessely(alpha,x), yx, 100*eps)
793 %!assert(besseli(alpha,x), ix, 100*eps)
794 %!assert(besselk(alpha,x), kx, 100*eps)
795 %!assert(besselh(alpha,1,x), jx + I*yx, 100*eps)
796 %!assert(besselh(alpha,2,x), jx - I*yx, 100*eps)
797 %!
798 %!assert(besselj(alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
799 %!assert(bessely(alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
800 %!assert(besseli(alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
801 %!assert(besselk(alpha,x,1), kx*exp(x), 100*eps)
802 %!assert(besselh(alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
803 %!assert(besselh(alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
804 %!
805 %! # Bessel functions, even order, complex x
806 %!
807 %! alpha = 2; x = 1.25 + 3.625 * I;
808 %! jx = -1.299533366810794494030065917 + 4.370833116012278943267479589*I;
809 %! yx = -4.370357232383223896393056727 - 1.283083391453582032688834041*I;
810 %! ix = -0.6717801680341515541002273932 - 0.2314623443930774099910228553*I;
811 %! kx = -0.01108009888623253515463783379 + 0.2245218229358191588208084197*I;
812 %!
813 %!assert(besselj(alpha,x), jx, 100*eps)
814 %!assert(bessely(alpha,x), yx, 100*eps)
815 %!assert(besseli(alpha,x), ix, 100*eps)
816 %!assert(besselk(alpha,x), kx, 100*eps)
817 %!assert(besselh(alpha,1,x), jx + I*yx, 100*eps)
818 %!assert(besselh(alpha,2,x), jx - I*yx, 100*eps)
819 %!
820 %!assert(besselj(alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
821 %!assert(bessely(alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
822 %!assert(besseli(alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
823 %!assert(besselk(alpha,x,1), kx*exp(x), 100*eps)
824 %!assert(besselh(alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
825 %!assert(besselh(alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
826 %!
827 %!assert(besselj(-alpha,x), jx, 100*eps)
828 %!assert(bessely(-alpha,x), yx, 100*eps)
829 %!assert(besseli(-alpha,x), ix, 100*eps)
830 %!assert(besselk(-alpha,x), kx, 100*eps)
831 %!assert(besselh(-alpha,1,x), jx + I*yx, 100*eps)
832 %!assert(besselh(-alpha,2,x), jx - I*yx, 100*eps)
833 %!
834 %!assert(besselj(-alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
835 %!assert(bessely(-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
836 %!assert(besseli(-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
837 %!assert(besselk(-alpha,x,1), kx*exp(x), 100*eps)
838 %!assert(besselh(-alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
839 %!assert(besselh(-alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
840 %!
841 %! # Bessel functions, odd order, complex x
842 %!
843 %! alpha = 3; x = 2.5 + 1.875 * I;
844 %! jx = 0.1330721523048277493333458596 + 0.5386295217249660078754395597*I;
845 %! yx = -0.6485072392105829901122401551 + 0.2608129289785456797046996987*I;
846 %! ix = -0.6182064685486998097516365709 + 0.4677561094683470065767989920*I;
847 %! kx = -0.1568585587733540007867882337 - 0.05185853709490846050505141321*I;
848 %!
849 %!assert(besselj(alpha,x), jx, 100*eps)
850 %!assert(bessely(alpha,x), yx, 100*eps)
851 %!assert(besseli(alpha,x), ix, 100*eps)
852 %!assert(besselk(alpha,x), kx, 100*eps)
853 %!assert(besselh(alpha,1,x), jx + I*yx, 100*eps)
854 %!assert(besselh(alpha,2,x), jx - I*yx, 100*eps)
855 %!
856 %!assert(besselj(alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
857 %!assert(bessely(alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
858 %!assert(besseli(alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
859 %!assert(besselk(alpha,x,1), kx*exp(x), 100*eps)
860 %!assert(besselh(alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
861 %!assert(besselh(alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
862 %!
863 %!assert(besselj(-alpha,x), -jx, 100*eps)
864 %!assert(bessely(-alpha,x), -yx, 100*eps)
865 %!assert(besseli(-alpha,x), ix, 100*eps)
866 %!assert(besselk(-alpha,x), kx, 100*eps)
867 %!assert(besselh(-alpha,1,x), -(jx + I*yx), 100*eps)
868 %!assert(besselh(-alpha,2,x), -(jx - I*yx), 100*eps)
869 %!
870 %!assert(besselj(-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
871 %!assert(bessely(-alpha,x,1), -yx*exp(-abs(imag(x))), 100*eps)
872 %!assert(besseli(-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
873 %!assert(besselk(-alpha,x,1), kx*exp(x), 100*eps)
874 %!assert(besselh(-alpha,1,x,1), -(jx + I*yx)*exp(-I*x), 100*eps)
875 %!assert(besselh(-alpha,2,x,1), -(jx - I*yx)*exp(I*x), 100*eps)
876 %!
877 %! # Bessel functions, fractional order, complex x
878 %!
879 %! alpha = 3.5; x = 1.75 + 4.125 * I;
880 %! jx = -3.018566131370455929707009100 - 0.7585648436793900607704057611*I;
881 %! yx = 0.7772278839106298215614791107 - 3.018518722313849782683792010*I;
882 %! ix = 0.2100873577220057189038160913 - 0.6551765604618246531254970926*I;
883 %! kx = 0.1757147290513239935341488069 + 0.08772348296883849205562558311*I;
884 %!
885 %!assert(besselj(alpha,x), jx, 100*eps)
886 %!assert(bessely(alpha,x), yx, 100*eps)
887 %!assert(besseli(alpha,x), ix, 100*eps)
888 %!assert(besselk(alpha,x), kx, 100*eps)
889 %!assert(besselh(alpha,1,x), jx + I*yx, 100*eps)
890 %!assert(besselh(alpha,2,x), jx - I*yx, 100*eps)
891 %!
892 %!assert(besselj(alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
893 %!assert(bessely(alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
894 %!assert(besseli(alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
895 %!assert(besselk(alpha,x,1), kx*exp(x), 100*eps)
896 %!assert(besselh(alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
897 %!assert(besselh(alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
898 %!
899 %! nix = 0.09822388691172060573913739253 - 0.7110230642207380127317227407*I;
900 %!
901 %!assert(besselj(-alpha,x), yx, 100*eps)
902 %!assert(bessely(-alpha,x), -jx, 100*eps)
903 %!assert(besseli(-alpha,x), nix, 100*eps)
904 %!assert(besselk(-alpha,x), kx, 100*eps)
905 %!assert(besselh(-alpha,1,x), -I*(jx + I*yx), 100*eps)
906 %!assert(besselh(-alpha,2,x), I*(jx - I*yx), 100*eps)
907 %!
908 %!assert(besselj(-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
909 %!assert(bessely(-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
910 %!assert(besseli(-alpha,x,1), nix*exp(-abs(real(x))), 100*eps)
911 %!assert(besselk(-alpha,x,1), kx*exp(x), 100*eps)
912 %!assert(besselh(-alpha,1,x,1), -I*(jx + I*yx)*exp(-I*x), 100*eps)
913 %!assert(besselh(-alpha,2,x,1), I*(jx - I*yx)*exp(I*x), 100*eps)
914 */
915
916 /*
633 ;;; Local Variables: *** 917 ;;; Local Variables: ***
634 ;;; mode: C++ *** 918 ;;; mode: C++ ***
635 ;;; End: *** 919 ;;; End: ***
636 */ 920 */