Sage functions arenot accepted by scipy fsolve

[Ne reknu]

Hello. 

I need a root finder of a real function which work with a initial guess. That can be done by scipy fsolve, but if I try to aply it on a gamma (or Bessel) function in Sage, I recieve following error message:

NotImplementedError: The Function gamma does not support numpy arrays as arguments

I have encoutered this before - then the solution was replacing  x[i]  by  x.item(i). I am however not skilled enough to do this chenge inside scipy library. Any sugestions?

Here is my code:

def P(l,m,x):
    return (-1)^m*gamma(l+m+1)*(1-x^2)^(m/2)/2^m/gamma(l-m+1)/gamma(m+1)*hypergeometric([m-l,m+l+1],[1+m],(1-x)/2)
dx=1e-10
def dP(l,m,x):
    return (P(l,m,x+dx) - P(l,m,x-dx))/(2*dx)
def K(nu,x):
    return bessel_K(nu,x)
def dK(nu,x):
    return (K(nu,x+dx)-K(nu,x-dx))/(2*dx)

 

import matplotlib.pyplot as plt
import numpy as np
import sympy
from scipy import optimize as opt

def Rce(KR):
    a=3
    m=0
    nu=sqrt(abs((1+a^2)*m^2-(a^2/4)))*I
    xD=a*KR
    ccD=1/sqrt(1+a^2)
    l=-0.5+sqrt(abs((1/4)-KR^2))
    return sqrt(1+a^2)*KR*P(l,m,ccD)*dK(nu,xD)+a*K(nu,xD)*dP(l,m,ccD)
    
    

ktip=0.2
k0=opt.fsolve(Rce,ktip)

Thanks for any advice.

[Nils Bruin]

On Friday, April 10, 2020 at 1:02:12 AM UTC-7, Ne reknu wrote:

Hello. 

I need a root finder of a real function which work with a initial guess. That can be done by scipy fsolve, but if I try to aply it on a gamma (or Bessel) function in Sage, I recieve following error message:

NotImplementedError: The Function gamma does not support numpy arrays as arguments

I have encoutered this before - then the solution was replacing  x[i]  by  x.item(i). I am however not skilled enough to do this chenge inside scipy library. Any sugestions?

Change the function you pass on instead, so that it does follow numpy's broadcasting rules: opt.fsolve(numpy.vectorize(Rce),ktip) probably does work.

While numpy's broadcasting rules are useful for numpy, they would probably not play nice with the other uses of functions in sage, so I don't think it's something we can accommodate by default. Users will have to rely on numpy's convenient wrapper.

[Ne reknu]

Thank you for your answer. Unfortunatelly the  numpy.vectorize(Rce) generates the very same error (I tried to plot it). Or pehaps I misunderstood your intention?

Dne pátek 10. dubna 2020 19:40:38 UTC+2 Nils Bruin napsal(a):

[Nils Bruin]

On Friday, April 10, 2020 at 1:05:15 PM UTC-7, Ne reknu wrote:

Thank you for your answer. Unfortunatelly the  numpy.vectorize(Rce) generates the very same error (I tried to plot it). Or pehaps I misunderstood your intention?

No, my guess would have been that would work, but then apparently it doesn't. It does give you a route to arrive at a solution, though. The problem here isn't sage-specific: it's just that normally, python functions don't broadcast over numpy arrays. There should be plenty of examples around the web that show how to get broadcast behaviour in place. It looked to me like numpy.vectorize would do the trick, but there should be other solutions, or perhaps you should use the vectorized function differently.

By realizing that the question is not sage-specific, I think it will be easy to find a ready-made solution elsewhere or find another forum with more specific expertise. Or perhaps someone else here knows more about numpy.