0
$\begingroup$

I have defined the Jacobi polynomials in terms of hypergeometric functions in two ways. While calculating the zeros of Jacobi polynomials, Mathematica shows different values of first few zeros. I am confused with the correct one. Any help regarding this will be highly appreciated. Thanks

sfact[a, n] === Pochhammer[a, n]

jacobi[n_, a_, b_, z_] :=  Pochhammer[a + 1, n]/n! Sum[Pochhammer[-n, k] Pochhammer[1 + n + a + b, k] (1 - z)^k/(Pochhammer[a + 1, k] k! 2^k), {k, 0, n}]

jacobii[n_, a_, b_, z_] = Pochhammer[1 + a, n])/n! * Hypergeometric2F1[-n, 1 + a + b + n, 1 + a, (1 - z)/2]

NSolve[jacobi[15, 8.3, -1.55, z] == 0, z]

{{z -> -1.00197}, {z -> -0.989872}, {z -> -0.946571}, {z ->-0.876951}, {z -> -0.781069}, {z -> -0.66255}, {z -> -0.524204}, {z -> -0.37}, {z -> -0.204121}, {z -> -0.0311152}, {z -> 0.144319}, {z -> 0.317478}, {z -> 0.483868}, {z -> 0.639596}, {z -> 0.78275}}

NSolve[jacobii[15, 8.3, -1.55, z] == 0, z]

{{z -> -1.00214}, {z -> -0.989586}, {z -> -0.946753}, {z -> -0.876849}, {z -> -0.781115}, {z -> -0.662534}, {z -> -0.524208}, {z -> -0.369999}, {z -> -0.204121}, {z -> -0.0311152}, {z -> 0.144319}, {z -> 0.317478}, {z -> 0.483868}, {z -> 0.639596}, {z -> 0.78275}}
Improve this question
$\endgroup$
8
  • $\begingroup$ Hypergeometric function is given by an infinite sum, but your first definition has only a finite number of terms. $\endgroup$
    – yarchik
    Jun 30, 2022 at 8:02
  • 3
    $\begingroup$ It's hard to give any help if you don't supply code other people can easily copy; no one is going to bother retyping stuff from a screenshot. $\endgroup$
    – J. M.'s persistent exhaustion ♦
    Jun 30, 2022 at 8:33
  • 3
    $\begingroup$ Please include code in copyable form. Please include symbolic results for a simple case, such as jacobi[4,2,7,z] and jacobii[4,2,7,z] and the built-in JacobiP[4,2,7,z] to see if, at the symbolic level, your expressions are probably correct. If symbolically everything is good, then welcome to floating point world. For instance z/.NSolve[JacobiP[15,5.3,-0.55,z]==0,z] and w+1/.NSolve[JacobiP[15,5.3,-0.55,w+1]==0,w] should be identical in theory but are not. These are polynomials with large coefficients and so on... $\endgroup$
    – user293787
    Jun 30, 2022 at 10:17
  • $\begingroup$ @yarchik The hypergeometric representation of the Jacobi polynomial is $$\frac{(a+1)_n}{n!} \sum_{k=0}^{n}\frac{(-n)_k (n+a+b+1)_k}{(a+1)_k k!}\left(\frac{1-x}{2}\right)^k$$ $\endgroup$
    – Pinaki Prasad Kar
    Jun 30, 2022 at 12:21
  • $\begingroup$ @user293787 Thanks a lot. Using JacobiP, I got the zeros as in the published paper. But still have confusions with the zeros of the polynomial defined manually. $\endgroup$
    – Pinaki Prasad Kar
    Jun 30, 2022 at 12:31

1 Answer 1

Reset to default
4
$\begingroup$ 
$Version

(* "13.1.0 for Mac OS X x86 (64-bit) (June 16, 2022)" *)

Clear["Global`*"]

sfact[a_, n_] := Product[(a + k - 1), {k, 1, n}]

Your sfact is just the function Pochhammer

sfact[a, n] === Pochhammer[a, n]

(* True *)

jacobi[n_, a_, b_, z_] := 
 Pochhammer[a + 1, n]/n! Sum[
   Pochhammer[-n, k] Pochhammer[1 + n + a + b, 
     k] (1 - z)^k/(Pochhammer[a + 1, k] k! 2^k), {k, 0, n}]

Your definition for jacobii is what you get for jacobi if you use Set rather than SetDelayed

jacobii[n_, a_, b_, z_] = jacobi[n, a, b, z]

(* (Hypergeometric2F1[-n, 1 + a + b + n, 1 + a, (1 - z)/2] * 
  Pochhammer[1 + a, n])/n!

The differences that you observed are precision issues

soli1 = NSolveValues[jacobi[15, 8.3, -1.55, z] == 0, z]

(* {-1.00164, -0.990395, -0.946253, -0.877127, -0.780995, -0.662574, -0.524198, \
-0.370001, -0.204121, -0.0311152, 0.144319, 0.317478, 0.483868, 0.639596, \
0.78275} *)

Note that the results for jacobi change with exact input

soli2 = NSolveValues[jacobi[15, 83/10, -155/100, z] == 0, z]

(* {-1.00223, -0.989436, -0.946852, -0.876792, -0.78114, -0.662525, -0.52421, \
-0.369999, -0.204122, -0.0311152, 0.144319, 0.317478, 0.483868, 0.639596, \
0.78275} *)

Similarly with jacobii

solii1 = NSolveValues[jacobii[15, 8.3, -1.55, z] == 0, z]

(* {-1.00214, -0.989586, -0.946753, -0.876849, -0.781115, -0.662534, -0.524208, \
-0.369999, -0.204121, -0.0311152, 0.144319, 0.317478, 0.483868, 0.639596, \
0.78275} *)

solii2 = NSolveValues[jacobii[15, 83/10, -155/100, z] == 0, z]

(* {-1.00223, -0.989436, -0.946852, -0.876792, -0.78114, -0.662525, -0.52421, \
-0.369999, -0.204122, -0.0311152, 0.144319, 0.317478, 0.483868, 0.639596, \
0.78275} *)

However, comparing the results from both methods using exact input

soli2 - solii2

(* {-9.66565*10^-11, -1.32849*10^-10, 2.83897*10^-11, 
 1.39112*10^-11, -5.60085*10^-12, 1.06459*10^-12, 
 3.10862*10^-14, -5.55112*10^-16, 8.32667*10^-17, -6.93889*10^-18, 0., 
 1.11022*10^-16, 3.33067*10^-16, -3.33067*10^-16, 2.22045*10^-16} *)

EDIT: The built-in JacobiP gives equivalent results with exact input.

jacobi[n, a, b, z]/JacobiP[n, a, b, z] == 
  jacobii[n, a, b, z]/JacobiP[n, a, b, z] == 1 // 
    FullSimplify

(* True *)

solP = NSolveValues[JacobiP[15, 83/10, -155/100, z] == 0, z];

solP - solii2

(* {-4.89355*10^-11, -1.84909*10^-10, 
 3.09808*10^-11, -3.90776*10^-12, -3.9404*10^-12, -6.37379*10^-13, 
 4.40759*10^-14, -8.88178*10^-16, 5.55112*10^-17, 1.04083*10^-17, 0., 0., 
 1.11022*10^-16, -9.99201*10^-16, 5.55112*10^-16} *)

To get better results use FullSimplify with the exact inputs.

solPfs = NSolveValues[FullSimplify[JacobiP[15, 83/10, -155/100, z]] == 0, z];

solii3 = NSolveValues[FullSimplify[jacobii[15, 83/10, -155/100, z]] == 0, z];

solPfs - solii3

(* {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.} *)
Improve this answer
$\endgroup$
4
  • $\begingroup$ Thank you for this wonderful clarification. $\endgroup$
    – Pinaki Prasad Kar
    Jun 30, 2022 at 12:54
  • $\begingroup$ But as suggested above in the comments, when we take the in-built syntax JacobiP it gives different value, what I exactly want. How does it works? $\endgroup$
    – Pinaki Prasad Kar
    Jul 1, 2022 at 0:42
  • $\begingroup$ Dear @PinakiPrasadKar, NSolve tries to numerically solve whatever equation you feed into it. To understand what is going on, one should not just compare the output of NSolve, but compare directly say jacobi[15,5.3,-0.55,z] and JacobiP[15,5.3,-0.55,z]. Perhaps look at something like the Wilkinson polynomial for numerical subtleties of numerical root finding. In any case, in Mathematica something like N[SolveValues[JacobiP[15,53/10,-55/100,z]==0,z],100] should give you $100$ correct digits and you should be fine. $\endgroup$
    – user293787
    Jul 1, 2022 at 5:10
  • $\begingroup$ @user293787 Thanks a lot for the explanation. $\endgroup$
    – Pinaki Prasad Kar
    Jul 2, 2022 at 7:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.