How to solve an equation that contains sigma sum to find the upper limit of sigma in Python

Question

I want to find an integer value from a sigma contained equation with two variables like this post

where x (a real decimal value) range is between two limits e.g. known_xmin_value <= x < known_xmax_value. 1 is the lower limit of k (which is the integer) but don't know the upper limit (which is the goal to be derived from the solution); Perhaps just could be guess that it will be lower than 2000000. Is it possible to find the integer for k? How?

from sympy import Sum, solve
from sympy.abc import k, x

known_value = 742.231       # just for example
known_xmin_value = 3.652    # just for example

solve(-known_value + Sum((x + (k - 1) * (x - known_xmin_value)) ** 3, (k, 1, integer_unknown_limit)), x)

score 1 · Answer 1 · answered Nov 16 '22 at 09:19

Since you are looking for a numerical solution, a naive approach would be a brute force over a set of integer_unknown_limit and use numerical root finding algorithms. For example:

from sympy import *
from scipy.optimize import root
import matplotlib.pyplot as plt

x, k, integer_unknown_limit = symbols("x, k, u")
known_value = 742.231       # just for example
known_xmin_value = 3.652    # just for example
expr = -known_value + Sum((x + (k - 1) * (x - known_xmin_value)) ** 3, (k, 1, integer_unknown_limit))

f = lambdify([x, integer_unknown_limit], expr)
res = {}

# NOTE: this loop may take a few minutes, depending on the machine.
# I suggest to lower the upper limit
for ul in range(1, 1000):
    print(ul)
    try:
        r = root(f, 2, args=(ul,))
        if r.success:
            res[ul] = r.x[0]
    except:
        pass

plt.figure()
plt.plot(list(res.keys()), res.values())
plt.xlabel("u (upper limit)")
plt.ylabel("x where f(x)=0")
plt.title("f(x) = $%s$" % latex(expr))
plt.show()

As you can see, for each upper limit value in the considered subset, the expression is satisfied at some value of x.

EDIT: note that the same can be done with sympy's nsolve. Just replace the for loop with:

for ul in range(1, 1000):
    print(ul)
    try:
        res[ul] = nsolve(expr.subs(integer_unknown_limit, ul), x, 1)
    except:
        pass

The `1` in `nsolve` is the initial guess. To improve performance you could break the for loop when the first solution is found. — Davide_sd, Nov 16 '22 at 13:45

score 1 · Accepted Answer · answered Nov 16 '22 at 11:51

SymPy can compute the sum symbolically:

from sympy import *

from sympy import Sum, solve
from sympy.abc import k, x, y

# Use exact rational numbers
known_value = Rational('742.231')       # just for example
known_xmin_value = Rational('3.652')    # just for example

sum_value = Sum((x + (k - 1) * (x - known_xmin_value)) ** 3, (k, 1, y))

This is the summation and its result computed using doit:

In [22]: sum_value
Out[22]: 
  y                           
 ____                         
 ╲                            
  ╲                          3
   ╲  ⎛            ⎛    913⎞⎞ 
   ╱  ⎜x + (k - 1)⋅⎜x - ───⎟⎟ 
  ╱   ⎝            ⎝    250⎠⎠ 
 ╱                            
 ‾‾‾‾                         
k = 1                         

In [23]: sum_value.doit()
Out[23]: 
                                       ⎛ 2    ⎞   ⎛      2                         ⎞ ⎛ 3    2    ⎞   ⎛ 4    3    2⎞ ⎛           2                        ⎞
761048497⋅y   ⎛2500707⋅x   2283145491⎞ ⎜y    y⎟   ⎜2739⋅x    2500707⋅x   2283145491⎟ ⎜y    y    y⎟   ⎜y    y    y ⎟ ⎜ 3   2739⋅x    2500707⋅x   761048497⎟
─────────── + ⎜───────── - ──────────⎟⋅⎜── + ─⎟ + ⎜─────── - ───────── + ──────────⎟⋅⎜── + ── + ─⎟ + ⎜── + ── + ──⎟⋅⎜x  - ─────── + ───────── - ─────────⎟
  15625000    ⎝  62500      15625000 ⎠ ⎝2    2⎠   ⎝  250       31250      15625000 ⎠ ⎝3    2    6⎠   ⎝4    2    4 ⎠ ⎝       250       62500      15625000⎠

This is a quartic polynomial in y so the equation can be solved explicitly in radicals giving 4 symbolic expressions for the roots:

In [24]: s1, s2, s3, s4 = solve(sum_value.doit() - known_value, y)

In [25]: print(s1)
(-250*x*sqrt(62500*x**2/(62500*x**2 - 456500*x + 833569) - 456500*x/(62500*x**2 - 456500*x + 833569) - 1000*sqrt(833569*x**2 + 185557750*x - 677656903)/(62500*x**2 - 456500*x + 833569) + 833569/(62500*x**2 - 456500*x + 833569)) - 250*x + 913*sqrt(62500*x**2/(62500*x**2 - 456500*x + 833569) - 456500*x/(62500*x**2 - 456500*x + 833569) - 1000*sqrt(833569*x**2 + 185557750*x - 677656903)/(62500*x**2 - 456500*x + 833569) + 833569/(62500*x**2 - 456500*x + 833569)) - 913)/(500*x - 1826)

For any particular value of x probably none of these 4 roots will be equal to an integer though.

How to solve an equation that contains sigma sum to find the upper limit of sigma in Python

2 Answers2