robert_dodier

sum_kron_delta -- summation of Kronecker delta in Maxima

Nov 13th, 2017
1,550
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. /* sum_kron_delta.mac -- simplify 'sum(f(i)*kron_delta(i, j), i, ...)
  2.  * copyright 2013 by Robert Dodier
  3.  * I release this work under terms of the GNU General Public License
  4.  */
  5.  
  6. put ('sum_kron_delta, true, 'present);
  7.  
  8. matchdeclare (kk, kron_delta_p);
  9. kron_delta_p (e) := not atom(e) and op(e) = 'kron_delta;
  10. matchdeclare (xx, lambda ([e], not kron_delta_p (e)));
  11. matchdeclare (ii, symbolp);
  12. matchdeclare ([ii1, iin], all);
  13.  
  14. simp : false;
  15. tellsimpafter ('sum (kk, ii, ii1, iin), FOO (kk, 1, ii, ii1, iin));
  16. tellsimpafter ('sum (kk * xx, ii, ii1, iin), FOO (kk, xx, ii, ii1, iin));
  17. simp:true;
  18.  
  19. FOO (kk, xx, ii, ii1, iin) := block ([ii_args],
  20.   kk : kron_delta_product (kk),
  21.   if atom(kk) or not (op(kk) = 'kron_delta)
  22.     then 'sum (kk*xx, ii, ii1, iin)
  23.     else
  24.      (ii_args : sublist (args (kk), lambda ([e], not freeof (ii, e))),
  25.       /* just take the first kron_delta variable other than ii -- this is not unique when there are 3 or more variables!! */
  26.       block ([a],
  27.         a : first (delete_all (ii_args, args (kk))),
  28.         /* oops -- what if ii_args has more (or less) than 1 element?? */
  29.         eq : solve (first (ii_args) = a, ii),
  30.         /* oops -- what if eq has more (or less) than 1 element?? */
  31.         buildq
  32.         ([eq,
  33.           S : index_set (ii1, iin),
  34.           kk_xx : subst (eq, kk * xx),
  35.           rhs_eq : rhs (first (eq))],
  36.           if %elementp (rhs_eq, S) then kk_xx else 0))));
  37.  
  38. kron_delta_product (k) := if op(k) = "*" then apply (append, args (k)) else k;
  39. delete_all (L1, L2) := (for x in L1 do L2 : delete (x, L2), L2);
  40. index_set (i, j) := %intersection (%range (i, j), 'integers);
  41.  
  42. matchdeclare (xx, all);
  43. matchdeclare (%ii, lambda ([e], not atom(e) and op(e) = '%intersection));
  44. tellsimp (%elementp (xx, %ii), apply ("and", map (lambda ([S], %elementp (xx, S)), args (%ii))));
  45.  
  46. matchdeclare (ii, lambda ([e], integerp (e) or featurep (e, 'integer)));
  47. tellsimp (%elementp (ii, 'integers), true);
  48. matchdeclare ([aa, bb], all);
  49. tellsimp (%elementp (xx, %range (aa, bb)), aa <= xx and xx <= bb);
RAW Paste Data