Continuity error with the function mathieu_modcem1

[brumann]

When looking at the shape of the derivative of the modified Mathieu function of order 1, it appears that there is a continuity problem.It sees to be a sign problem, I thought it would come from a square root branch mistake, but I cannot find that in the Mathieu function defined in Abramowitz and stegun.

Here are the parameters to use to obtain the graph given:

• q vary from 5 to 30 (200 points)
• m = 6
• x = 0.856
• function ussed sp.mathieu_modcem1(m, q, x)[1]

Thank you in advance for your response.

Sincerely,

[rlucas7]

Thanks for reporting the issue. Can you post the code that you used to produce the plot please? I think I see what is happening, would like to confirm. Sincerely,

…

-Lucas Roberts

On May 29, 2020, at 5:21 AM, brumann ***@***.***> wrote:  When looking at the shape of the derivative of the modified Mathieu function of order 1, it appears that there is a continuity problem.It sees to be a sign problem, I thought it would come from a square root branch mistake, but I cannot find that in the Mathieu function defined in Abramowitz and staging. Here are the parameters to use to obtain the graph given: q vary from 5 to 30 (200 points) m = 6 x = 0.856 function ussed sp.mathieu_modcem1(m, q, x)[1] Thank you in advance for your response. Sincerely, — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.

[brumann]

Thank you for your answer
Here is the code to obtain the graph given above.

import matplotlib.pyplot as plt
import numpy as np
import scipy.special as sp

q = np.linspace(5,30,200)
m = 6
x = 0.856
f_x = sp.mathieu_modcem1(m, q, x)[1]
plt.plot(q,f_x)
plt.xlabel('q')
plt.title('Derivative $Ce_m(x,q)$')
plt.show()

In fact for my code I was just using the mathieu fonction inside a Newton's method and this one was finding an erroneous value.

Sincerely

[rlucas7]

On my side I can confirm that I see the same thing in master w/python 3.6

>>> scipy.__version__
'1.6.0.dev0+f57e5bf'

unfortunately it isn't what I thought it was and looks like it goes into the fortran codes.

[rlucas7]

Some notes that may prove helpful to those who speak fortran:

mathieu_modcem1 is a wrapper here:

scipy/scipy/special/functions.json

Line 948 in f57e5bf

"mathieu_modcem1": {

which leads to the source .c file:

scipy/scipy/special/specfun_wrappers.c

Line 525 in f57e5bf

int mcm1_wrap(double m, double q, double x, double *f1r, double *d1r)

and that calls the fortran:

scipy/scipy/special/specfun_wrappers.c

Line 537 in f57e5bf

F_FUNC(mtu12,MTU12)(&kf,&kc,&int_m, &q, &x, f1r, d1r, &f2r, &d2r);

and the fortran source starts here:

scipy/scipy/special/specfun/specfun.f

Line 11872 in 42e50b4

SUBROUTINE MTU12(KF,KC,M,Q,X,F1R,D1R,F2R,D2R)

not immediately clear to me why an q swapping from 25.1xxx to 25.2xxx would do this but I'm not familiar with fortran or with the function in particular.

>>> q[160]
25.100502512562812
>>> q[161]
25.22613065326633

[ianthomas23]

The explanation of this problem is that the algorithm requires finding the root of a function that has multiple roots, and within a certain parameter region the initial estimate for the location of the root is poor such that the secant method root finding converges to a different root.

The Fortran function MTU12 calls CVA2

scipy/scipy/special/specfun/specfun.f

Lines 1402 to 1406 in 7815f95

	SUBROUTINE CVA2(KD,M,Q,A)
	C
	C ======================================================
	C Purpose: Calculate a specific characteristic value of
	C Mathieu functions

with KD=1, M=6 and Q=25.0 works fine but Q=25.25 shows the problem. CVA2 calls CV0 to obtain the initial estimate for A

scipy/scipy/special/specfun/specfun.f

Lines 764 to 769 in 7815f95

	SUBROUTINE CV0(KD,M,Q,A0)
	C
	C =====================================================
	C Purpose: Compute the initial characteristic value of
	C Mathieu functions for m ≤ 12 or q ≤ 300 or
	C q ≥ m*m

and then calls REFINE to improve it

scipy/scipy/special/specfun/specfun.f

Lines 1797 to 1801 in 7815f95

	SUBROUTINE REFINE(KD,M,Q,A)
	C
	C =====================================================
	C Purpose: calculate the accurate characteristic value
	C by the secant method

It is best illustrated by the following plot.

The code is trying to find a root of the blue line, i.e. where it crosses the orange line, and the initial estimate is indicated by the green line. For Q=25.0 the initial estimate is just within the central region of the plot and the secant method successfully finds the nearest root at about A=49.0. For Q=25.25the initial estimate is not within the central region of the plot so the secant method locates the root to the left at about A=27 rather than the required root at about A=49.25.

According to the file header, all this Fortran code was taken from the book "Computation of Special Functions" by Zhang and Jin, 1996, Wiley and Sons. To proceed any further with this probably requires the assistance of someone who has a copy of the book.