subscripted variables, sum, diff

[rtoy]

Imported from SourceForge on 2024-07-04 23:22:37 Created by zmth on 2020-01-02 10:27:45 Original: https://sourceforge.net/p/maxima/bugs/3604

---

build_info(version="5.43.0",
timestamp="2019-06-01 18:55:31",
host="x86_64-w64-mingw32",
lisp_name="SBCL",
lisp_version="1.4.14",
maxima_userdir="C:/Users/zmth/maxima",
maxima_tempdir="C:/Users/zmth/AppData/Local/Temp",
maxima_objdir="C:/Users/zmth/maxima/binary/5_43_0/sbcl/1_4_14",
maxima_frontend="wxMaxima",
maxima_frontend_version="19.05.7")

(declare([b,d],integer),
assume(b>0,d>1),
sum(diff(sum(x_i^2,i,1,d)^(b/2),x_i),i,1,d))

every thing else being equal ,since the symbol x_i always means exactly the same as x[i] then the above should give same answer as

(declare([b,d],integer),
assume(b>0,d>1),
sum(diff(sum(x[i]^2,i,1,d)^(b/2),x[i]),i,1,d))

but it does NOT and

(declare([b,d],integer),
assume(b>0,d>1),
diff(sum(x_i^2,i,1,d)^(b/2),x_i))

also gives a ridiculous answer and both

(declare([b,d],integer),
assume(b>0,d>1),
t:sum(x[i]^2,i,1,d)^(b/2),
diff(t,x[1]))

(declare([b,d],integer),
assume(b>0,d>1),
diff(sum(x_i^2,i,1,d)^(b/2),x_1))

give WRONG answer of 0

[rtoy]

Imported from SourceForge on 2024-07-04 23:22:38 Created by robert_dodier on 2020-01-06 07:05:23 Original: https://sourceforge.net/p/maxima/bugs/3604/#fc58

---

There is a bug here: Maxima gives an incorrect result for differentiating a summation wrt to a subscripted variable. E.g. diff(sum(f(x[i]), i, 1, n), x[j]) should be diff(f(x[j]), x[j]) but Maxima says it's 0.

I'll post a longer comment showing how to use some add-on packages to get the desired result.

[rtoy]

Imported from SourceForge on 2024-07-04 23:22:42 Created by robert_dodier on 2020-01-06 07:06:36 Original: https://sourceforge.net/p/maxima/bugs/3604/#b508

---

This is an interesting problem. Maxima lacks three parts which are needed for it: chain rule for derivatives, derivative of a sum, and summation of Kronecker delta. I'll help Maxima with that stuff in order to get to the desired result.

Maxima doesn't generally apply the chain rule for differentiation. Here's a simplistic rule to apply the chain rule in a form that's needed for this problem.

(%i2) matchdeclare (ee, lambda ([e], not mapatom(e) and not mapatom(first(e)))) $
(%i3) matchdeclare (xx, mapatom) $
(%i4) tellsimpafter ('diff(ee,xx,1), block([foo:gensym()],
                                      ('at('diff(subst(first(ee) = foo,ee),foo,1),foo = first(ee)))
                                     *'diff(first(ee),xx))) $

Maxima incorrectly computes the derivative of a summation wrt to an arbitrary subscripted variable; to remedy that we'll use the package diff_sum. Also, Maxima can't simplify summations containing the Kronecker delta; for that we'll use the package sum_kron_delta. These packages are part of the project maxima-packages on Github: https://github.com/maxima-project-on-github/maxima-packages/tree/master/robert-dodier

(%i5) load ("sum_kron_delta.mac") $
(%i6) load ("diff_sum.mac") $

Now let's see what we get for summation stated in this bug report. I have written the summation as 'sum(...) and the derivative as 'diff(...) since diff_sum works by simplifying noun expressions.

(%i7) declare ([b, d], integer) $
(%i8) assume (b > 0, d > 1) $
(%i9) 'sum ('diff ('sum (x[i]^2, i, 1, d)^(b/2), x[j]), j, 1, d);
          d
         ====
         \
(%o9) 2 ( >    (if true then x  else 0))
         /                    j
         ====
         j = 1
                                              !
                              d           b/2 !          d
                          (------- (g19020   )!         ====    )
                           dg19020            !         \      2
                                              !g19020 =  >    x
                                                        /      i
                                                        ====
                                                        i = 1

Hmm, looks like we have something promising here, let's encourage Maxima a little by asking for evaluation of the nouns in that expression.

(%i10) %, nouns;
                      d                d
                     ====             ====
                     \      2 b/2 - 1 \
(%o10)            b ( >    x )         >    x
                     /      i         /      j
                     ====             ====
                     i = 1            j = 1

That looks right, let's check it for d equal to 3.

(%i11) %, d = 3;
                      3                3
                     ====             ====
                     \      2 b/2 - 1 \
(%o11)            b ( >    x )         >    x
                     /      i         /      j
                     ====             ====
                     i = 1            j = 1
(%i12) %, nouns;
                              2    2    2 b/2 - 1
(%o12)       (x  + x  + x ) (x  + x  + x )        b
               3    2    1    3    2    1

If we compute the inner summation as a + expression, and then differentiate it, and add up the terms, we should get the same result.

(%i13) sum (x[i]^2, i, 1, d)^(b/2), d = 3;
                          2    2    2 b/2
(%o13)                  (x  + x  + x )
                          3    2    1
(%i14) makelist (diff (%, x[k]), k, 1, 3);
             2    2    2 b/2 - 1         2    2    2 b/2 - 1
(%o14) [x  (x  + x  + x )        b, x  (x  + x  + x )        b, 
         1   3    2    1             2   3    2    1
                                           2    2    2 b/2 - 1
                                      x  (x  + x  + x )        b]
                                       3   3    2    1
(%i15) lsum (u, u, %);
            2    2    2 b/2 - 1          2    2    2 b/2 - 1
(%o15) x  (x  + x  + x )        b + x  (x  + x  + x )        b
        3   3    2    1              2   3    2    1
                                            2    2    2 b/2 - 1
                                     + x  (x  + x  + x )        b
                                        1   3    2    1
(%i16) factor (%);
                              2    2    2 b/2 - 1
(%o16)       (x  + x  + x ) (x  + x  + x )        b
               3    2    1    3    2    1

Is it the same? Yes!

(%i17) is (%o16 = %o12);
(%o17)                        true