Solving Integral with changing range with Python

256 Views Asked by At

I'm new to python but I been working on a code which can solve an integral equation which range is also changing according the unknown parameter. I tried to use Sympy solve function but it does not return any result. I solved the problem with a for loop, but its really slow, inefficient. I'm sure there must be a better solution, a solver. Maybe there is an other approach? I'am messing up something? I attach also the code.


import sympy as sy
from sympy.solvers import solve

alphasum = 1.707
Lky = 3.078
g = 8
Ep = 195
sigp1 = 1401.927
sigp0 = 1476
e = 2.718282
u = 0.05
k = 0.007

lsl = sy.Symbol('lsl')


def lefthand(g, Ep):
    return g * Ep


def rigthhand(sigp1, sigp0, e, u, alphasum, Lky, k, lsl):
    return (sigp1 - (-sigp1 + 2 * sigp0 - 2 * sigp0 * (1 - e ** (-u * (((alphasum / Lky) * lsl) + k * lsl)))))


equr = (sy.integrate(rigthhand(sigp1, sigp0, e, u, alphasum, Lky, k, lsl), (lsl, 0, lsl)))
equl = lefthand(g, Ep)

print(equr)
print (equl)
print (equr-equl)

result = solve(equr-equl, lsl,warn =True,check=False,minimal=True,quick=True,simplify=True,quartics=True)

print(result)
2

There are 2 best solutions below

2
On

You can use nsolve to find numerical solutions:

In [11]: sympy.nsolve(equr-equl, lsl, 8)                                                                                          
Out[11]: 8.60275245116315

In [12]: sympy.nsolve(equr-equl, lsl, -4)                                                                                         
Out[12]: -4.53215114758428
0
On

When u is 0 your equation is linear and you can solve if for lsl; let this be an initial guess for the value of lsl at a larger value of u. Increase your value of u to the target value. Let's use capital U for the parameter we will control, replacing it with values closer and closer to u:

>>> U = Symbol('U')
>>> equr = (sy.integrate(rigthhand(sigp1, sigp0, e, U, alphasum, Lky, k, lsl), (lsl, 0, lsl)))
>>> equl = lefthand(g, Ep)
>>> z = equr-equl
>>> u0 = solve(z.subs(U,0),lsl)[0]
>>> for i in range(1,10):  # get to u in 9 steps
...   u0 = nsolve(z.subs(U,i*u/10.), u0)
>>> print(u0)
-4.71178322070344

Now define a larger value of u

>>> u = 0.3
>>> for i in range(1,10):  # get to u in 9 steps
...    u0 = nsolve(z.subs(U,i*u/10.), u0)
...
>>> u0
-2.21489271112540

Our intitial guess will (usually) be pretty good since it is exact when u is 0; if it fails then you might need to take more than 9 steps to reach the target value.