0

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)
dejanualex
  • 3,872
  • 6
  • 22
  • 37
bgabor
  • 1

2 Answers2

0

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
Oscar Benjamin
  • 12,649
  • 1
  • 12
  • 14
  • Thanks for the quick reply! It seams to be working! If I understand the number is a close value to the answer. But what if for example I want to use this for mora data, and don't really know the range or a close value to the answer, what can a I do? Use maybe the library checksol and make a loop till its true? – bgabor May 20 '20 at 20:30
  • Sympy does not yet have a global equation solver: https://github.com/sympy/sympy/issues/19164#issuecomment-632022578 – Oscar Benjamin May 21 '20 at 11:00
0

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.

smichr
  • 16,948
  • 2
  • 27
  • 34