2

I'm having trouble getting my head around the following numerical evaluation. I have a function with two variables, r and gamma. I now wish to plot the root of that function as a function of gamma from 0 to 500. I know how to use lambdify to evaluate a function over a 1d array of a variable. However, I do not understand how to make sympy find the root of the function for a given variable gamma and then step over many gamma values from 0 to 500. 

import numpy as np
import sympy as sm
import sympy.mpmath as mpmath
from sympy.interactive import printing
printing.init_printing(use_latex=True)

r, gamma = sm.symbols("r gamma")

function = -0.5+(sm.exp(-2*r**2)+sm.exp(-1-r**2*(4+sm.exp(-2*r**2*gamma)))*r**2*gamma)/(sm.exp(2*r**2)+r**2*gamma)  #the function with two variables
root_gammazero = sm.nsolve(function.subs(gamma,0),0.1) # this gives r = 0.416277

x_vector = np.linspace(0,500,1000)    #step through parameter gamma on x-axis

functionlambdify = sm.lambdify(gamma, sm.nsolve(function.subs(gamma,gamma),r, 0.1), "numpy") #find root for a given gamma and then return that root in the y_vector
y_vector = functionlambdify(x_vector)

How do I pass the changing variable gamma into the lambdify function. In the code example I tried this by using function.subs(gamma,gamma),r, 0.1) in functionlambdify = sm.lambdify(gamma, sm.nsolve(function.subs(gamma,gamma),r, 0.1), "numpy") but it doesnt work.

Any help is greatly appreciated.

Thanks a lot!

Kokomoking
  • 363
  • 1
  • 4
  • 5

1 Answers1

2

nsolve is essentially an algorithm implemented in Python so you cannot "lambdify" it. Instead, just write an explicit loop calling nsolve for each of your elements. e.g:

y_vector = np.empty(len(x_vector))
for i in range(len(x_vector)):
    y_vector[i] = sm.nsolve(function.subs(gamma,x_vector[i]), 0.1)

If this solution is too slow for you, you could create a callback of your expression using lambdify and use a non-linear solver from e.g. scipy.

Bjoern Dahlgren
  • 931
  • 8
  • 18