4

I would like SymPy to evaluate an expression like the following:

formula

How would I define the symbols and the expression so that SymPy could handle it nicely? I would like to keep N as just a symbol, i.e. not make an actual finite list of x's. I have tried various combinations of IndexedBase and Sum /Product, but didn't get it working right.

Bendik
  • 1,097
  • 1
  • 8
  • 27
  • It is not difficult to compute that analytically. Why even use SymPy? – Rodrigo de Azevedo Jan 23 '18 at 11:09
  • By the way, the desired result is twice the sum of the `x`'s. – Rodrigo de Azevedo Jan 23 '18 at 11:20
  • 1
    @RodrigodeAzevedo The expression I gave would be part a larger one, which I cut out in the interest of keeping the question about SymPy, rather than my specific problem. I have done it all analytically, but I would like to make sure it is correct. – Bendik Jan 23 '18 at 11:31

1 Answers1

3

Ideally it would be this:

x = IndexedBase("x")
i, j, N = symbols("i j N")
expr = Sum(Product(exp(-x[j]**2), (j, 1, N)).diff(x[i]), (i, 1, N))

So far this is unevaluated, expr is 

Sum(Derivative(Product(exp(-x[j]**2), (j, 1, N)), x[i]), (i, 1, N))

The method doit can be used to evaluate it. Unfortunately the differentiation of a product doesn't quite work yet: expr.doit() returns

N*Derivative(Product(exp(-x[j]**2), (j, 1, N)), x[i])

Rewriting the product as the sum prior to differentiation helps: 

expr = Sum(Product(exp(-x[j]**2), (j, 1, N)).rewrite(Sum).diff(x[i]), (i, 1, N))
expr.doit()

returns 

Sum(Piecewise((-2*exp(Sum(log(exp(-x[j]**2)), (j, 1, N)))*x[i], (1 <= i) & (i <= N)), (0, True)), (i, 1, N))

which is the correct result of differentiation. Sadly we have that extraneous condition in Piecewise, and also log(exp(...)) that should have been simplified. SymPy doesn't infer that (1 <= i) & (i <= N) is True from the context of the outer sum, and it also hesitates to simplify log(exp thinking x[j] might be complex. So I resort to surgical procedure with Piecewise, replacing it by the first piece, and to forcefully expanding logs:

e = expr.doit()
p = next(iter(e.atoms(Piecewise)))
e = expand_log(e.xreplace({p: p.args[0][0]}), force=True)

Now e is 

Sum(-2*exp(Sum(-x[j]**2, (j, 1, N)))*x[i], (i, 1, N))

Couldn't get exp(Sum(..)) to become a Product again, unfortunately.

  • Also, `factor(e)` recognizes that the double sum can be factored as a product of two sums. –  Jan 22 '18 at 15:59