robert_dodier

diff_sum -- derivative of summation wrt an indexed variable

Nov 13th, 2017
1,528
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. /* diff_sum.mac -- simplify 'diff('sum(f(x[i]), i,. ..), x[j])
  2.  * copyright 2013 by Robert Dodier
  3.  * I release this work under terms of the GNU General Public License
  4.  */
  5.  
  6. put ('diff_sum, true, 'present);
  7.  
  8. matchdeclare (mm, lambda ([e], not atom(e) and op(e)="*"));
  9. matchdeclare (xx, all);
  10. foo : 'diff(mm,xx);
  11. tellsimpafter (''foo,
  12.                block([aa:args(mm), nn:length(mm)],
  13.                      sum (product (if jj=ii then 'diff(aa[ii], xx) else aa[jj], jj, 1, nn), ii, 1, nn))) ;
  14.  
  15. find_subscripts (x, e) := block ([L : []], find_subscripts1 (x, e), L);
  16. find_subscripts1 (x, e) := block ([L1 : match_subscript (e, x)],
  17.   if L1 = false
  18.     then
  19.       if not mapatom (e)
  20.         then map (lambda ([ee], find_subscripts1 (x, ee)),  e)
  21.         else e
  22.     else (push (rhs (first (L1)), L), e));
  23. matchdeclare (xx, symbolp);
  24. matchdeclare (bb, all);
  25. defmatch (match_subscript, xx[bb], xx);
  26. push (a, b) ::= buildq ([a, b], b : cons (a, b));
  27.  
  28. matchdeclare (aa, all);
  29. matchdeclare (ii1, "<" (minf));
  30. matchdeclare (nn, ">" (inf));
  31. matchdeclare ([ii, xx], symbolp);
  32. matchdeclare (jj, all);
  33.  
  34. simp : false $
  35.  
  36. tellsimpafter ('diff ('sum (aa, ii, ii1, nn), xx[jj], 1),
  37.   block ([L, diff_summand],
  38.     L : find_subscripts (xx, aa),
  39.     /* diff_summand : apply ("+", makelist ('diff (aa, xx[j1], 1) * kron_delta (jj, j1), j1, L)), ?? */
  40.     diff_summand : apply ("+", makelist (diff (aa, xx[j1], 1) * kron_delta (jj, j1), j1, L)),
  41.     buildq ([diff_summand, ii, ii1, nn],
  42.       'sum (diff_summand, ii, ii1, nn))));
  43.  
  44. simp : true $
  45.  
  46. /* Find subscripts of expressions such as x[i](t).
  47.  * Dunno what to call it; "find_subscripts_with_function_calls" is clumsy and still not very precise.
  48.  * "find_subscripts_2" it is, then.
  49.  */
  50. find_subscripts_2 (x, t, e) := block ([L : []], find_subscripts1_2 (x, t, e), L);
  51. find_subscripts1_2 (x, t, e) := block ([L1 : match_subscript_2 (e, x, t)],
  52.   if L1 = false
  53.     then
  54.       if not mapatom (e)
  55.         then map (lambda ([ee], find_subscripts1_2 (x, t, ee)),  e)
  56.         else e
  57.     else (push (rhs (first (L1)), L), e));
  58. matchdeclare ([xx, tt], symbolp);
  59. matchdeclare (bb, all);
  60. defmatch (match_subscript_2, xx[bb](tt), xx, tt);
  61.  
  62. simp : false $
  63.  
  64. tellsimpafter ('diff ('sum (aa, ii, ii1, nn), xx[jj](tt), 1),
  65.   block ([L, diff_summand],
  66.     L : find_subscripts_2 (xx, tt, aa),
  67.     /* diff_summand : apply ("+", makelist ('diff (aa, xx[j1](tt), 1) * kron_delta (jj, j1), j1, L)), ?? */
  68.     diff_summand : apply ("+", makelist (diff (aa, xx[j1](tt), 1) * kron_delta (jj, j1), j1, L)),
  69.     buildq ([diff_summand, ii, ii1, nn],
  70.       'sum (diff_summand, ii, ii1, nn))));
  71.  
  72. simp : true $
RAW Paste Data