# diff_sum -- derivative of summation wrt an indexed variable

Nov 13th, 2017
3,652
0
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 \$