scipy icon indicating copy to clipboard operation
scipy copied to clipboard

ENH: Jacobi elliptic functions with negative second argument

Open miflynn opened this issue 1 year ago • 3 comments

Is your feature request related to a problem? Please describe.

I am working with a class of differential equations whose solutions are described by Jacobi elliptic functions. A representative example of such a differential equation is $x''(t) = -ax(t) - bx^{3}(t)$ with $a, b$ positive real constants. Equations of this type describe an oscillator with a quartic nonlinearity in physics.

Describe the solution you'd like.

Solutions to this class of differential equations can be written in terms of the Jacobi sine function $\text{sn}(u,m)$ with $m<0$. Currently, SciPy supports these functions with the restriction that $0<m<1$. My ideal solution to this problem would be for SciPy to support computation of the Jacobi elliptic functions for $m<1$, as is available in other libraries such as mpmath.

Describe alternatives you've considered.

No response

Additional context (e.g. screenshots, GIFs)

Currently, evaluating the SciPy Jacobi functions evaluate to nan when $m<0$.

miflynn avatar Aug 19 '24 16:08 miflynn

Hi @miflynn, thanks for the request. I did a quick check and extending the Jacobi elliptic functions to $m < 0$ would actually be pretty straightforward. For negative $m$ we want to take a look at DLMF 22.17, in particular DLMF 22.17 (5-8). In their notation $m = k^2$, and $\mathrm{sd} = \frac{\mathrm{sn}}{\mathrm{dn}}$, $\mathrm{cd} = \frac{\mathrm{cn}}{\mathrm{dn}}$, $\mathrm{nd} = \frac{1}{\mathrm{dn}}$.

I've checked that with 22.17.5 and 22.17.6 we can compute values of $\mathrm{sn}(u, m)$ for negative $m$ that are reasonably close to mpmath, and I'm pretty sure the formulas for $\mathrm{dn}$ and $\mathrm{cn}$ will also work as expected. There are formulas for $m > 1$ in DLMF 22.17, but they do not agree with mpmath, perhaps because they give different branches.

I unfortunately don't have bandwidth to work on this at the moment but could probably get to it within a month or two. Let me you know if you'd be interested in submitting a PR, if so I can help get you pointed in the right direction.

steppi avatar Aug 21 '24 11:08 steppi

hi @steppi, can i give this a go ( if you are not looking into it right now ) ? I can go through a bit of the code flow and get back if i have any queries this week.

dannyi96 avatar Aug 26 '24 17:08 dannyi96

hi @steppi, can i give this a go ( if you are not looking into it right now ) ? I can go through a bit of the code flow and get back if i have any queries this week.

That would be great @dannyi96. I’d be happy to help with any questions you have and review your PR when it’s ready.

steppi avatar Aug 27 '24 00:08 steppi