# solve equation with double sum This post is a wiki. Anyone with karma >750 is welcome to improve it.

Hi!

1) f(x)==h(x)

2) g(x)+S_{i,j,k}(x) == 0

I know I can solve (numerically) eq.(1) doing:

x=var('x')
find_root(f(x)==h(x),x,x_min,x_max)


In eq.(2) S_{i,j,k}(x) is a triple sum function of 'x' and i,j and k are the index of the sum.

How can I solve (numerically) eq.(2)?

Update:

If I run the next code:

import sympy.mpmath

N=20
A=0.7
G_0 = 37.7
B = 0.36

x = sympy.symbols('x')

def S(x_):
return sympy.mpmath.nsum(lambda i, j, k: (12*A**4*x_**6*i**4-30*A**2*x_**3*i**2*(j**2+k**2)+3*(j**2+k**2)**2)/(2*(A**2*x_**3*i**2+j**2+k**2)**(7/2)),[1,N],[1,N],[1,N])

def F(x_):
return G_0 * (x_ - 1/(x_**2))

print(sympy.mpmath.findroot(F(x) + B*A*sqrt(x)*S(x), [0.85,1]) )


I get the next error:

TypeError: unsupported operand parent(s) for '*': 'Symbolic Ring' and '<class 'sympy.mpmath.ctx_mp_python.mpf'>'


What am I doing wrong?

Best regards!

edit retag close merge delete

Can you post more details on the sum? We will need more details to be able to give much help.

Sort by » oldest newest most voted

The following code running with Sage rather than mpmath runs for me, but it says there are no zeros on the interval.

var('x')

N=20
A=0.7
G_0 = 37.7
B = 0.36

def S(x_):
return add([(12*A**4*x_**6*i**4-30*A**2*x_**3*i**2*(j**2+k**2)+3*(j**2+k**2)**2)/(2*(A**2*x_**3*i**2+j**2+k**2)**(7/2)) for i in range(1,N+1) for j in range(1,N+1) for k in range(1,N+1)])

def F(x_):
return G_0 * (x_ - 1/(x_**2))

print(find_root(F(x) + B*A*sqrt(x)*S(x), 0.85,1) )

more

Thanks!! I fix it defining the function G(x)=F(x) + B*A*sqrt(x)*S(x) becausesympy.mpmath.findroot need that. One more question about the summation: I wanna use the limits -N to N but i want to exclude the point (i,j,k)=(0,0,0) (this is the tree index zero simultaneously) . How can I do that? Best regards. Thanks again

Try:

add([(12*A**4*x_**6*i**4-30*A**2*x_**3*i**2*(j**2+k**2)+3*(j**2+k**2)**2)/(2*(A**2*x_**3*i**2+j**2+k**2)**(7/2)) for i in range(-N,N+1) for j in range(-N,N+1) for k in range(-N,N+1) if (i!=0 or j!=0 or k!=0) ] )