Software: Freud machine by Alphonse P. Magnus math. , Univ. of Louvain (Belgium) magnus@anma.ucl.ac.be http://www.math.ucl.ac.be/~magnus/ Aug. 1996 See here FORTRAN programs associated to the numerical approximate computation of the recurrence coefficients a1, a2, ... , aN of the orthonormal polynomials p0, p1, ... : a_{n+1} p_{n+1}(x) = x p_n(x) -a_n p_{n-1}(x) , n=0,1,... (a0=0), related to the (even) weight function |x|^\rho \exp(-P(x^2)) on the whole real line. Elementary recurrence transformations also show that a_{2n+2} a_{2n+1} q_{n+1}(y) = (y- a_{2n+1}^2 -a_{2n}^2 ) q_n(y) -a_{2n} a_{2n-1} q_{n-1}(y) is the recurrence relation of the orthonormal polynomials q_n(y) related to the weight y^{(\rho -1)/2} \exp(-P(y)) on the positive real line (take y=x^2). Let x_n = a_n^2. The programs build and use nonlinear recurrence relations F_n(x) = n + \rho (1-(-1)^n)/2 (FREUD'S EQUATIONS [1-6]) for the unknowns. F_n is linear in the coefficients of P. If m is the degree of P, F_n contains sums of products of up to m factors taken in the sequence x_{n-m+1}, ... , x_{n+m-1}. For instance, if P(y)=ay^2+by, F_n(x)= 4a x_n(x_{n-1}+x_n+x_{n+1}) +2bx_n. Program freud1.f asks for an highest power m and puts relevant data in a condensed way in the file freud00m.dat . You don't need to use freud1 if you already have a file freud00m.dat with m high enough. If you need m>10, you will have to experiment with the dimensions inside freud1.f ... freud2.f prints the data from a file freud00m.dat in a readable (???!!!) way. It gives all the terms of F_n, something like (-3) (4) must be understood as x_{n-3} x_{n+4}. It then gives all the partial derivatives (needed in freud3) of F_n with respected to x_{n-m+1}, ... , x_n. Actually, the given derivative of F_n with respect to x_n is in fact the derivative of (F_n/x_n). Well, try with a low value of m. freud3.f tries to find the an's from \rho and the coefficients of P (\rho must be > -1 and the main coefficient of P must be positive). The coefficients are computed between n=0 and a given maximum index N. Try several values of N, as the an's with n close to N are estimated from a crude asymptotic formula. The method is close to Lew & Quarles [2] one: a trial sequence x_0, ..., x_N is corrected according to Newton-Raphson method (whence the need for the derivatives of F_n). Successive error norms are printed. THE PROGRAM freud3.f IS BASED ON THE CONJECTURE THAT THE MATRIX OF DERIVATIVES IS A POSITIVE DEFINITE SYMMETRIC MATRIX near the solution. The matrix is known to be OK if P has positive coefficients. If Cholesky factorization attempt fails at an intermediate step, the program gives up. Don't worry, restart with a less offensive polynomial P, and ask again for P after completion, as each new run starts with the solution of the preceding one. For instance, P(y)=y^3-10y was solved by solving successively y^3-y, y^3-2y, etc. (see below). When all is well, output is, from n=0 to n=N: x_n = a_n^2 a_n if n is odd: a_{n-1}^2+a_n^2 a_{n-2}^2 a_{n-2}^2 [1] G. Freud , On the coefficients in the recursion formulae of orthogonal polynomials, Proc. Roy. Irish Acad. Sect.A (1) vol. 76 (1976),1-6. [2] J.S. Lew, D.A. Quarles, Nonnegative solutions of a nonlinear recurrence, J. Approx. Theory 38 (1983) 357-379. [3] A.P. Magnus, A proof of Freud's conjecture about orthogonal polynomials related to $|x|^\rho \exp(-x^{2m})$, for integer $m$, pp. 362-372 in 'Polynomes orthogonaux et applications, Proceedings, Bar-le-Duc 1984', edited by C.Brezinski, A.Draux, A.P.Magnus, P.Maroni, and A.Ronveaux, Springer Lect. Notes Math. 1171, Springer, Berlin, 1985. [4] A.P. Magnus, On Freud's equations for exponential weights, J. Approx. Theory 46 (1986) 65-99. [5] A.P. Magnus, Refined asymptotics for Freud's recurrence coefficients, pp.196-200 in 'Linear and Complex Analysis Problem Book 3, Part II' edited by V.P. Havin and N.K. Nikolski, Springer Lect. Notes Math. 1574, Springer, Berlin, 1994. [6] A.P. Magnus, Freud's equations as discrete Painleve equations. See entry in http://www.math.ucl.ac.be/~magnus/cvtimes.htm Example: runs of freud2 and freud3: Script started on Fri Aug 30 10:05:21 1996 /u18/grpanma/magnus/freud @ ux12 [2] %freud2 Freud's equations for x^(2m), 1<=m<=10 . Give m. (stop if m outside 1..10) 3 freud010.dat freud008.dat verbose output? (1=yes, 0=no) 1 1 3 (terms) F_n(x^ 4)= 4 x_n times + 1( -1)( + 1( 0)( + 1( 1)( 1 (terms) d F_n/d log x_(n -1) = 4 x_n x_(n -1) times + 1( 1 (terms) d F_n/d log x_(n 0) = 4 x_n x_(n 0) times + 1( 13 8 (terms) F_n(x^ 6)= 6 x_n times + 1( -2)( -1)( + 1( -1)( -1)( + 2( -1)( 0)( + 1( -1)( 1)( + 1( 0)( 0)( + 2( 0)( 1)( + 1( 1)( 1)( + 1( 1)( 2)( 1 (terms) d F_n/d log x_(n -2) = 6 x_n x_(n -2) times + 1( -1)( 4 (terms) d F_n/d log x_(n -1) = 6 x_n x_(n -1) times + 1( -2)( + 2( -1)( + 2( 0)( + 1( 1)( 3 (terms) d F_n/d log x_(n 0) = 6 x_n x_(n 0) times + 2( -1)( + 2( 0)( + 2( 1)( STOP: /u18/grpanma/magnus/freud @ ux12 [3] %freud3 max degree ? 3 freud010.dat freud008.dat 2 3 1 -1 1 0 1 1 1 .. 1 1 3 2 -1 2 0 2 1 degree (end if <0), number of recurrence coef.? 3 50 give rho and polynomial coef. (starting with largest degree) 0.3 1 0 -10 0.30 60.0 ; 1.0000 x 3 0.0000 x 2 -10.0000 x 1 0.5054 0.5673 0.6070 0.6368 0.6609 0.6813 0.6990 0.7148 0.7289 0.7418 -4.4163 -4.6030 -4.4295 -4.1534 -3.7702 -3.3116 -2.7968 -2.2385 -1.6453 -1.0232 -3.9643 -3.4276 -2.6450 -1.7648 -0.7792 0.2810 1.3968 2.5559 3.7498 4.9723 err Cholesky n= 1 -3.96430896485237 sorry, try another intermediate polynomial STOP: /u18/grpanma/magnus/freud @ ux12 [4] %freud3 max degree ? 3 freud010.dat freud008.dat 2 3 1 -1 1 0 1 1 1 .. 1 1 3 2 -1 2 0 2 1 degree (end if <0), number of recurrence coef.? 3 25 give rho and polynomial coef. (starting with largest degree) 0.3 1 0 -1 0.30 60.0 ; 1.0000 x 3 0.0000 x 2 -1.0000 x 1 0.5054 0.5673 0.6070 0.6368 0.6609 0.6813 0.6990 0.7148 0.7289 0.7418 0.1816 1.1899 2.2017 3.1452 4.0920 5.0433 5.9986 6.9572 7.9187 8.8826 0.6335 2.3653 3.9863 5.5338 7.0830 8.6358 10.1922 11.7516 13.3137 14.8781 1.7379 0.4004 0.7499 0.6241 0.7123 0.6842 0.7364 0.7190 0.7595 0.7467 178.9354 11.7344 5.1385 4.2951 5.5407 6.1027 7.4603 8.0843 9.4329 10.0747 527.1357 12.8891 9.3589 6.9508 9.7567 10.0934 12.8488 13.3319 16.0285 16.5753 2 iter. err 1: 2.47016407117108 1.4750 0.3574 0.7209 0.6286 0.7033 0.6860 0.7315 0.7198 0.7564 0.7471 65.7723 5.1433 3.6691 4.0199 5.3073 6.0016 7.3022 8.0007 9.3011 10.0004 196.6053 5.6960 7.0488 6.6611 9.3006 9.9882 12.5393 13.2368 15.7758 16.4870 3 iter. err 1: 0.328060799308702 1.2537 0.3460 0.7220 0.6278 0.7033 0.6860 0.7315 0.7198 0.7564 0.7471 24.4270 2.8369 3.3587 4.0004 5.3000 6.0000 7.3000 8.0000 9.3000 10.0000 74.5709 3.2176 6.7334 6.6278 9.2899 9.9872 12.5345 13.2365 15.7735 16.4867 4 iter. err 1: 0.325129801636409 1.0678 0.3796 0.7164 0.6286 0.7033 0.6859 0.7315 0.7198 0.7564 0.7471 9.5105 2.1753 3.3182 4.0027 5.3001 6.0000 7.3000 8.0000 9.3000 10.0000 29.5478 2.6233 6.6460 6.6320 9.2942 9.9853 12.5352 13.2364 15.7735 16.4867 5 iter. err 1: 0.320968572406916 0.9149 0.4436 0.7010 0.6323 0.7029 0.6858 0.7316 0.7198 0.7564 0.7471 4.1786 2.0994 3.3431 4.0103 5.3005 6.0000 7.3000 8.0000 9.3000 10.0000 12.8693 2.8081 6.4956 6.6668 9.2972 9.9815 12.5370 13.2359 15.7736 16.4867 6 iter. err 2: 0.311551165896484 0.8006 0.4972 0.6856 0.6356 0.7027 0.6857 0.7316 0.7198 0.7564 0.7471 2.1988 2.0765 3.3323 4.0074 5.3005 6.0000 7.3000 8.0000 9.3000 10.0000 6.5777 3.0725 6.3042 6.6854 9.3016 9.9775 12.5388 13.2355 15.7736 16.4867 7 iter. err 1: 0.266827202802609 0.7356 0.5243 0.6772 0.6371 0.7026 0.6856 0.7317 0.7198 0.7564 0.7471 1.4886 2.0278 3.3096 4.0019 5.3001 6.0000 7.3000 8.0000 9.3000 10.0000 4.3551 3.1837 6.1862 6.6875 9.3048 9.9751 12.5398 13.2353 15.7736 16.4867 8 iter. err 1: 0.169460872532181 0.7160 0.5313 0.6751 0.6375 0.7026 0.6856 0.7317 0.7198 0.7564 0.7471 1.3137 2.0026 3.3008 4.0001 5.3000 6.0000 7.3000 8.0000 9.3000 10.0000 3.8198 3.1984 6.1536 6.6871 9.3057 9.9744 12.5400 13.2353 15.7736 16.4867 9 iter. err 1: 5.405943684023323E-002 0.7144 0.5318 0.6750 0.6375 0.7026 0.6856 0.7317 0.7198 0.7564 0.7471 1.3001 2.0000 3.3000 4.0000 5.3000 6.0000 7.3000 8.0000 9.3000 10.0000 3.7788 3.1986 6.1513 6.6870 9.3058 9.9744 12.5401 13.2353 15.7736 16.4867 10 iter. err 1: 4.411331647301447E-003 0.7144 0.5318 0.6750 0.6375 0.7026 0.6856 0.7317 0.7198 0.7564 0.7471 1.3000 2.0000 3.3000 4.0000 5.3000 6.0000 7.3000 8.0000 9.3000 10.0000 3.7786 3.1986 6.1512 6.6870 9.3058 9.9744 12.5401 13.2353 15.7736 16.4867 11 iter. err 1: 2.710121906404358E-005 0.7144 0.5318 0.6750 0.6375 0.7026 0.6856 0.7317 0.7198 0.7564 0.7471 1.3000 2.0000 3.3000 4.0000 5.3000 6.0000 7.3000 8.0000 9.3000 10.0000 3.7786 3.1986 6.1512 6.6870 9.3058 9.9744 12.5401 13.2353 15.7736 16.4867 12 iter. err 1: 1.019885273038833E-009 0.7144 0.5318 0.6750 0.6375 0.7026 0.6856 0.7317 0.7198 0.7564 0.7471 1.3000 2.0000 3.3000 4.0000 5.3000 6.0000 7.3000 8.0000 9.3000 10.0000 3.7786 3.1986 6.1512 6.6870 9.3058 9.9744 12.5401 13.2353 15.7736 16.4867 13 iter. err 24: 2.671080549140649E-016 ----------- result ----------- 3 0.3000000000 1.0000000000 0.0000000000 -1.0000000000 25 0 0.00000000000000 0 0.00000000000000 1 0.51031522934387 1 0.71436351344667 0 0.51031522934387 0.00000000000000 2 0.28281938818713 2 0.53180766089549 3 0.45557338208857 3 0.67496176342706 1 0.73839277027569 0.14432704094561 4 0.40642452299943 4 0.63751433160317 5 0.49360757529984 5 0.70257211394976 2 0.90003209829927 0.18515619450658 6 0.47004458806768 6 0.68559797845945 7 0.53537559121634 7 0.73169364574003 3 1.00542017928403 0.23201756939890 8 0.51806264481447 8 0.71976568743895 9 0.57208578440791 9 0.75636352662454 4 1.09014842922237 0.27735809475465 10 0.55817971546891 10 0.74711425864382 etc. 25 0.77677891851490 25 0.88135062178165 12 1.52039882722557 0.55924978836586 degree (end if <0), number of recurrence coef.? 3 30 give rho and polynomial coef. (starting with largest degree) 0.3 1 0 -2 0.30 60.0 ; 1.0000 x 3 0.0000 x 2 -2.0000 x 1 etc. 8 iter. err 1: 2.372758655702095E-014 ----------- result ----------- 3 0.3000000000 1.0000000000 0.0000000000 -2.0000000000 30 0 0.00000000000000 0 0.00000000000000 1 0.66930197948376 1 0.81810878218227 0 0.66930197948376 0.00000000000000 ... etc. etc. and so on with 1 0 -3 , 1 0 -4 , .... , up to degree (end if <0), number of recurrence coef.? 3 75 give rho and polynomial coef. (starting with largest degree) 0.3 1 0 -10 0.30 60.0 ; 1.0000 x 3 0.0000 x 2 -10.0000 x 1 1.2964 0.2474 1.2450 0.3741 1.1783 0.4985 1.0930 0.6274 1.0043 0.7392 -2.0614 1.8776 0.2001 3.7201 2.5232 5.5030 4.9106 7.2127 7.2829 8.9072 56.9810 2.0255 50.6805 4.4436 43.6246 7.6015 36.3727 12.0420 31.0669 17.6558 1.3402 0.2276 1.3056 0.3291 1.2826 0.4014 1.2823 0.4686 1.2529 0.5722 1.8320 1.9931 4.1844 4.0113 7.5409 6.1935 13.2371 8.7676 17.5352 11.4492 73.3735 2.1076 69.1926 4.4981 69.7130 7.2684 78.9374 10.7546 80.1259 15.5603 2 iter. err 8: 0.583699462128753 1.3344 0.2324 1.2924 0.3437 1.2481 0.4357 1.2103 0.5193 1.1592 0.6231 1.3110 2.0001 3.3534 4.0030 5.6141 6.0169 8.5841 8.0496 11.3258 10.1281 71.1149 2.1227 65.0242 4.5635 59.9486 7.4055 58.1336 10.7369 54.7007 15.3148 3 iter. err 8: 0.205580676791651 1.3342 0.2327 1.2909 0.3469 1.2386 0.4516 1.1769 0.5580 1.1036 0.6709 1.3000 2.0000 3.3010 4.0008 5.3284 6.0102 7.5395 8.0425 9.8238 10.0907 71.0585 2.1232 64.6501 4.5777 57.8199 7.5687 51.2985 11.4339 45.0559 16.6769 4 iter. err 10: 0.147866785858131 1.3342 0.2327 1.2908 0.3473 1.2370 0.4551 1.1682 0.5694 1.0853 0.6858 1.3000 2.0000 3.3000 4.0000 5.3011 6.0008 7.3159 8.0045 9.3482 10.0095 71.0576 2.1232 64.6242 4.5793 57.5049 7.5982 49.6899 11.6211 42.1889 17.0677 5 iter. err 10: 4.376878527495005E-002 1.3342 0.2327 1.2908 0.3474 1.2369 0.4554 1.1675 0.5705 1.0835 0.6873 1.3000 2.0000 3.3000 4.0000 5.3000 6.0000 7.3001 8.0000 9.3005 10.0001 71.0576 2.1232 64.6227 4.5794 57.4844 7.6005 49.5607 11.6383 41.9068 17.1071 6 iter. err 10: 4.363450260503076E-003 1.3342 0.2327 1.2908 0.3474 1.2369 0.4554 1.1674 0.5705 1.0835 0.6873 1.3000 2.0000 3.3000 4.0000 5.3000 6.0000 7.3000 8.0000 9.3000 10.0000 71.0576 2.1232 64.6227 4.5794 57.4842 7.6005 49.5597 11.6385 41.9041 17.1074 7 iter. err 10: 4.132156762264493E-005 1.3342 0.2327 1.2908 0.3474 1.2369 0.4554 1.1674 0.5705 1.0835 0.6873 1.3000 2.0000 3.3000 4.0000 5.3000 6.0000 7.3000 8.0000 9.3000 10.0000 71.0576 2.1232 64.6227 4.5794 57.4842 7.6005 49.5597 11.6385 41.9041 17.1074 8 iter. err 10: 3.777243220638750E-009 1.3342 0.2327 1.2908 0.3474 1.2369 0.4554 1.1674 0.5705 1.0835 0.6873 1.3000 2.0000 3.3000 4.0000 5.3000 6.0000 7.3000 8.0000 9.3000 10.0000 71.0576 2.1232 64.6227 4.5794 57.4842 7.6005 49.5597 11.6385 41.9041 17.1074 9 iter. err 6: 9.197137813722951E-016 ----------- result ----------- 3 0.3000000000 1.0000000000 0.0000000000 -10.0000000000 75 0 0.00000000000000 0 0.00000000000000 1 1.78018428338462 1 1.33423546774347 0 1.78018428338462 0.00000000000000 2 0.05415780697857 2 0.23271829962118 3 1.66610650745614 3 1.29077748177451 1 1.72026431443471 0.09641087680584 4 0.12065400887015 4 0.34735285930902 5 1.52989105387482 5 1.23688764804036 2 1.65054506274497 0.20102242932923 6 0.20741735258589 6 0.45543095259972 7 1.36293063347265 7 1.16744620153249 3 1.57034798605854 0.31732595213955 8 0.32547255692874 8 0.57050202184457 9 1.17387692417975 9 1.08345600934221 4 1.49934948110849 0.44359651819285 10 0.47237557155785 10 0.68729583991019 etc. 74 1.14956997887706 74 1.07218001234730 75 1.24737147513298 75 1.11685785807012 37 2.39694145401004 1.35958367485531 degree (end if <0), number of recurrence coef.? -1 1 STOP: /u18/grpanma/magnus/freud @ ux12 [10] %exit script done on Fri Aug 30 10:16:18 1996