By Mark van Hoeij

Posted: Feb. 29, 2024

The recursion relation inside Doron's program K(n, x) (included at the end of this file) has a polynomial solution:
    P := (n-x+1)*(n-x)*(3*n^2+2*x*n+3*n+x^2+x);
(Several computer algebra systems have code to find polynomial/rational/hypergeometric solutions.)

One can construct a new recurrence L such that the solutions of L are:
    the solutions of the original recurrence, divided by P.
Since the original recurrence had P among its solutions, the new recurrence L will have 1 among its solutions.
We can write L in Q(x, n)[Sn] where Sn is the shift operator, which sends f(n,x) to f(n+1,x).

Saying that 1 is a solution of L is equivalent to saying that L factors like this:
    L = L1 * (Sn - 1)
for some L1 in Q(x, n)[Sn] where * = composition of operators.

If u is a solution of L1, then sum(u) is a solution of L, and P * sum(u) is a solution of the original recurrence (the recurrence for K(n,x)). Here `sum(u)' refers to an expression such that if you apply Sn - 1 to it, you get u. In Maple, `sum(u)' would be something like this:
    ui := subs(n=i, u); add( ui, i=0..n-1 );
because if you apply Sn - 1 to that you get:
    add(ui, i=0..n) - add(ui, i=0..n-1) = subs(i = n, ui) = u.
After computing P and u, and putting them into a program, I get a program MyK(n, x) for which MyK(n, x) = K(n, x) for n = 1, 2. But since MyK(n, x) satisfies the same recurrence, we get MyK(n, x) = K(n, x) for all positive integers n.

However, if I try to compute MyK(n, n+1) then I get a division by zero error. The reason is because u, the solution of L1, which turned out to be a rational function, it has a pole in the summation interval when x = n+1.

To solve this problem, suppose that v is some rational function for which sum(v) is also rational. Then sum(u) is the same as sum(v) + sum(u - v).

By choosing v in a good way we we can do a couple of things:
    (1) choose v to cancel out that pole when x = n+1.
    (2) pick v in such a way that u - v is as small as possible.
That produces the following improved MyK(n, x) program:

MyK := proc(n,x) local i; factor(
    4*(x-n-1)*(x-n)*(3*n^2+2*n*x+x^2+3*n+x) * 'add'(i/(i^2+i*x+x^2-1)/(x-i), i=1..n)
    + n/45*(n+1)*(250*n^4-456*n^3*x-198*n^2*x^2+494*n*x^3-45*x^4+500*n^3-684*n^2*x
      -198*n*x^2+247*x^3+92*n^2-2*n*x+45*x^2-158*n+113*x)/(x^3-x) )
end:

One can verify that the formula inside this improved program MyK(n,x) still satisfies the 2nd order recurrence, and that this improved MyK(n, x) still equals K(n, x) for n = 1, 2, and hence for all positive integers n. But this time, if x = n+1 then MyK(n, x) does not produce a division by zero error. And MyK(n-1, n) returns n*(n+1)^2.

For completeness: here is the original K(n,x) program that contains the initial values K(1,x) and K(2,x) as well as the order-2 recurrence relation:

K:=proc(n,x) option remember: if n=1 then 2/15*(15*x^3+262*x^2+115*x-588)/(x+1)/x: elif n=2 then 2/15*(45*x^6+1325*x^5+1372*x^4-10076*x^3-8557*x^2+1599*x+33876)/(x^2+2*x+3)/x/(x+1)/(x-1): else normal( (500*n^13-456*n^12*x-196*n^11*x^2-1706*n^10*x^3+1726*n^9*x^4+946*n^8 *x^5+1402*n^7*x^6-2084*n^6*x^7-1304*n^5*x^8-226*n^4*x^9+814*n^3*x^10+554*n^2*x^ 11+30*n*x^12-3250*n^12+2736*n^11*x+1078*n^10*x^2+8530*n^9*x^3-7767*n^8*x^4-3784 *n^7*x^5-4907*n^6*x^6+6252*n^5*x^7+3260*n^4*x^8+452*n^3*x^9-1221*n^2*x^10-554*n *x^11-15*x^12+6000*n^11-3996*n^10*x-4480*n^9*x^2-12635*n^8*x^3+9184*n^7*x^4+ 10758*n^6*x^5+7519*n^5*x^6-6275*n^4*x^7-6810*n^3*x^8-4124*n^2*x^9+547*n*x^10+ 262*x^11+2750*n^10-5100*n^9*x+12075*n^8*x^2-640*n^7*x^3+4102*n^6*x^4-19030*n^5* x^5-6530*n^4*x^6+2130*n^3*x^7+6955*n^2*x^8+3898*n*x^9-70*x^10-19000*n^9+15517*n ^8*x-13586*n^7*x^2+10746*n^6*x^3-14138*n^5*x^4+16200*n^4*x^5+12799*n^3*x^6+7957 *n^2*x^7-2701*n*x^8-1374*x^9+11250*n^8-1372*n^7*x-1253*n^6*x^2+5828*n^5*x^3+ 6967*n^4*x^4-5098*n^3*x^5-15122*n^2*x^6-7980*n*x^7+300*x^8+16000*n^7-17938*n^6* x+13177*n^5*x^2-16440*n^4*x^3-9334*n^3*x^4-6474*n^2*x^5+5499*n*x^6+2550*x^7-\ 17750*n^6+7100*n^5*x-7025*n^4*x^2+5360*n^3*x^3+14263*n^2*x^4+6482*n*x^5-330*x^6 -1500*n^5+7589*n^4*x+3227*n^3*x^2+2803*n^2*x^3-5233*n*x^4-2026*x^5+7000*n^4-\ 3364*n^3*x-4875*n^2*x^2-1846*n*x^3+115*x^4-2000*n^3-716*n^2*x+1858*n*x^2+588*x^ 3)/(n-1)/(n-1-x)/(n^2+n*x+x^2-1)/(250*n^9-228*n^8*x-98*n^7*x^2-603*n^6*x^3+635* n^5*x^4+375*n^4*x^5+98*n^3*x^6-407*n^2*x^7-277*n*x^8-15*x^9-2250*n^8+1824*n^7*x +686*n^6*x^2+3618*n^5*x^3-3175*n^4*x^4-1500*n^3*x^5-294*n^2*x^6+814*n*x^7+277*x ^8+8250*n^7-5521*n^6*x-2693*n^5*x^2-8410*n^4*x^3+5715*n^3*x^4+2885*n^2*x^5+946* n*x^6-362*x^7-15750*n^6+7590*n^5*x+6605*n^4*x^2+9520*n^3*x^3-4445*n^2*x^4-2770* n*x^5-750*x^6+16500*n^5-4057*n^4*x-9145*n^3*x^2-5870*n^2*x^3+797*n*x^4+965*x^5-\ 9000*n^4-324*n^3*x+6503*n^2*x^2+2348*n*x^3+473*x^4+2000*n^3+716*n^2*x-1858*n*x^ 2-588*x^3)*K(n-1,x) -n*(n-x)*(250*n^9-228*n^8*x-98*n^7*x^2-603*n^6*x^3+635*n^5*x^4+375*n ^4*x^5+98*n^3*x^6-407*n^2*x^7-277*n*x^8-15*x^9-750*n^7+863*n^6*x-635*n^5*x^2+ 635*n^4*x^3-635*n^3*x^4+635*n^2*x^5+652*n*x^6+45*x^7+750*n^5-1042*n^4*x+635*n^3 *x^2-635*n^2*x^3-473*n*x^4-45*x^5-250*n^3+407*n^2*x+98*n*x^2+15*x^3)*(n^2+n*x+x ^2-2*n-x)/(n-1)/(n-1-x)/(n^2+n*x+x^2-1)/(250*n^9-228*n^8*x-98*n^7*x^2-603*n^6*x ^3+635*n^5*x^4+375*n^4*x^5+98*n^3*x^6-407*n^2*x^7-277*n*x^8-15*x^9-2250*n^8+ 1824*n^7*x+686*n^6*x^2+3618*n^5*x^3-3175*n^4*x^4-1500*n^3*x^5-294*n^2*x^6+814*n *x^7+277*x^8+8250*n^7-5521*n^6*x-2693*n^5*x^2-8410*n^4*x^3+5715*n^3*x^4+2885*n^ 2*x^5+946*n*x^6-362*x^7-15750*n^6+7590*n^5*x+6605*n^4*x^2+9520*n^3*x^3-4445*n^2 *x^4-2770*n*x^5-750*x^6+16500*n^5-4057*n^4*x-9145*n^3*x^2-5870*n^2*x^3+797*n*x^ 4+965*x^5-9000*n^4-324*n^3*x+6503*n^2*x^2+2348*n*x^3+473*x^4+2000*n^3+716*n^2*x -1858*n*x^2-588*x^3)*K(n-2,x)): fi: end:

ShouldBeZeros := seq(normal(K(n,x) - MyK(n,x)), n=1..5); # (It suffices to check just n = 1..2)

(Note that one could shorten K(n,x) by tossing the initival value K(2,x) and replacing it with the initial value K(0,x) = 0.)

back