1

I have a system of equations expressed by sympy:

def test1a(A):
    t, tt = sym.symbols('t tt')
    return sym.cos(t+tt+A)*A

def test1b(B):
    t, tt = sym.symbols('t tt')
    return sym.sin(t-tt+B)*B

That I want to convert to a numpy expression before evaluating the result with fsolve():

def testprep(variables, A, B):
    
    t, tt = sym.symbols('t tt')
    
    fA = lambdify([(t, tt)], test1a(A))
    fB = lambdify([(t, tt)], test1b(B))
    
    t, tt = variables
    
    return [fA,fB]

def testsolve(A, B):
    
    print(so.fsolve(testprep, [-1, 1], args = (A, B)))
    
    return

When ran, I get the following result:

import scipy.optimize as so
import sympy as sym
from sympy.utilities.lambdify import lambdify as lambdify
import numpy as np

def main():
    A = 1
    B = 2
    testsolve(A,B)
    return

if __name__ == '__main__':
    main()

Output:

error: Result from function call is not a proper array of floats.

As a sanity check, I drafted the same system in terms of numpy expressions and solved it:

def standard(variables, A, B):
    
    t, tt = variables
    
    eq1 = np.cos(t+tt+A)*A
    eq2 = np.sin(t-tt+B)*B
    
    return [eq1, eq2]

def solvestandard(A, B):
    
    print(so.fsolve(standard, np.array([-1, 1]), args = (A,B)))
    
    return

With output:

[-0.71460184  1.28539816]

I'm new to lambdify( ) and am not too familiar with the process of converting from sympy to numpy. What would I need to change in order to make the test-case work?

t.o.
  • 832
  • 6
  • 20

1 Answers1

2

You've presented your code in a very convoluted way for an SO question. You don't need so many functions just to show what is basically 5 lines of code!

Please reduce your examples to something as simple as possible that is a single block of code that can be copy-pasted complete with all imports. Then please test it by copy-pasting and running in a fresh python process.

The way to use lambdify with fsolve is something like:

import sympy as sym
import scipy.optimize as so

t, tt = sym.symbols('t tt')
A, B = 1, 2
eq1 = sym.cos(t+tt+A)*A
eq2 = sym.sin(t-tt+B)*B
f = lambdify([(t, tt)], [eq1, eq2])
print(so.fsolve(f, [-1, 1]))

The idea is that lambdify makes an efficient function that can be computed many times (e.g. fsolve will call it iteratively). Then you pass that efficient function to fsolve.

The function that you pass to fsolve should not call lambdify itself (as your testprep does) because lambdify is a lot slower than evaluating the function:

In [22]: %time f = lambdify([(t, tt)], [eq1, eq2])
CPU times: user 4.74 ms, sys: 77 µs, total: 4.82 ms
Wall time: 4.96 ms

In [23]: %time f([1, 2])
CPU times: user 36 µs, sys: 1e+03 ns, total: 37 µs
Wall time: 41 µs
Out[23]: [-0.6536436208636119, 1.682941969615793]

Maybe your real problem is more complicated but I would use sympy's nsolve for something like this:

In [16]: nsolve([eq1, eq2], [t, tt], [-1, 1])
Out[16]: 
⎡-0.714601836602552⎤
⎢                  ⎥
⎣ 1.28539816339745 ⎦

In [17]: nsolve([eq1, eq2], [t, tt], [-1, 1], prec=100)
Out[17]: 
⎡-0.714601836602551690384339154180124278950707650156223544756263851923045898428447750342991293664470733⎤
⎢                                                                                                      ⎥
⎣1.285398163397448309615660845819875721049292349843776455243736148076954101571552249657008706335529267 ⎦

SymPy's nsolve takes care of calling lambdify for you. It is somewhat slower than fsolve because it works with arbitrary precision floating point but that means it can also compute more accurate results.

Oscar Benjamin
  • 12,649
  • 1
  • 12
  • 14
  • `help(f)` should show the `numpy` code that was created. Sometimes `lambdify` produces a bad result because it does a rather naive lexical translation (without a 'deep' knowledge of `numpy`). – hpaulj Jul 07 '21 at 20:29
  • Thanks for the pointers, this is a great answer. Coding in a more pythonic way is something I struggle with, but your explanation addresses my main problems. If nsolve can accept additional arguments in the form of an array of constants it may be something I look into. – t.o. Jul 07 '21 at 20:29
  • Rather than passing in additional arguments just substitute those values into `eq1` and `eq2` before calling `nsolve`. – Oscar Benjamin Jul 08 '21 at 10:08
  • I.e. if `A` and `B` are symbols then you can do `eq1.subs({A:5, B:4})`. – Oscar Benjamin Jul 08 '21 at 10:09