Блог

У меня есть сложная функция ниже:

  import numpy as np
 import scipy as sp
 from scipy.special import jv, hankel1, jvp, h1vp, h2vp, gamma

 nr = 2
 m = 20

 def Dm(x):
     return nr * jvp(m,nr*x,1) * hankel1(m,x) - jv(m,nr*x) * h1vp(m,x,1)
  

Я ищу, чтобы найти столько комплексных корней Dm (x), которые лежат в 4-м квадранте комплексной плоскости, сколько смогу, используя newton() из scipy.optimize, а затем сохранить их в одномерный массив. Лучший метод, который я могу придумать, — это перебирать его с помощью newton() с регулярными интервалами через конечную часть 4-го квадранта, проверить, является ли корень дубликатом предыдущего корня, проверить, действительно ли корень является корнем, а затем сохранить его в массив. Как только этот алгоритм завершится, я хочу отсортировать массив путем увеличения реальных компонентов. Мои вопросы:

(i) Могу ли я создать массив неопределенной длины, к которому я могу продолжать добавлять значения по мере их нахождения?

(ii) Могу ли я построить функцию таким образом, чтобы я мог визуализировать корни? Математика говорит, что все они находятся на одном листе комплексной плоскости.

(iii) Есть ли лучший метод поиска корней? Я чувствую, что с моим методом я пропущу много корней в домене.

Некоторые ответы:

(i) использовать список. Массивы имеют фиксированный размер. Добавление в список — очень дешевый вариант. Когда вы добавляете новый корень в список, убедитесь, что предыдущий корень отсутствует в списке, например, вычисляя np.amin(np.abs(np.array(a)-b)) , где a находится список существующих корней и b является новым корнем. Если это значение очень мало, вы попали в существующий корень. (Насколько мал, зависит от функции. Оно не может быть 0.0, так как тогда вы обычно не распознаете один и тот же корень из-за неточностей с плавающей запятой и итерации.)

Если у вас очень большое количество корней (тысячи), вы можете захотеть отсортировать их, как только вы их получите. Это ускоряет поиск соответствующих корней. С другой стороны, скорее всего, > 90% времени тратится на итерацию корней, и вам не нужно беспокоиться о других проблемах с производительностью. Затем вы просто компилируете список, сортируете его (сортировка списка проста и быстра) и преобразуете в массив, если вам нужно.

(ii) Да. Два примера ниже: (За countour материал спасибо Уоррену Векессеру и его очень хороший ответ!)

 import numpy as np
from scipy.special import jv, hankel1, jvp, h1vp
import matplotlib.pyplot as plt


nr = 2
m = 20

# create a 2000 x 2000 sample complex plane between -2-2i ..  2 2i
x = np.linspace(-2, 2, 2000)
y = np.linspace(-2, 2, 2000)
X, Y = np.meshgrid(x, y)
C = X   1j * Y 
z = 1-C**2

# draw a contour image of the imaginary part (red) and real part (blue)
csr = plt.contour(x, y, z.real, 5, colors='b')
plt.clabel(csr)
csi = plt.contour(x, y, z.imag, 5, colors='r')
plt.clabel(csi)
plt.axis('equal')
plt.savefig('contours.png')

# draw an image of the absolute value of the function, black representing zeros
plt.figure()
plt.imshow(abs(z), extent=[-2,2,-2,2], cmap=plt.cm.gray)
plt.axis('equal')
plt.savefig('absval.png')
  

Это дает countours.png :

и absval.png :

Обратите внимание, что если вы хотите увеличить изображения, обычно необходимо изменить пределы и пересчитать комплексные значения z , чтобы избежать пропущенных деталей. Изображения, конечно, могут быть нанесены друг на друга, цветовая палитра изображения может быть изменена, countour имеет миллион вариантов. Если вы хотите рисовать только нули, замените число 5 (количество контуров) на [0] (отображать только указанный контур) в countour вызовах.

Конечно, вы замените my (1-C ^ 2) своей собственной функцией. Единственное, о чем нужно позаботиться, это то, что если функция получает массив комплексных чисел, она возвращает массив результатов в той же форме, вычисленной точечно. Imshow должен получить массив скаляров. Для получения дополнительной информации, пожалуйста, обратитесь к imshow документации.

(iii) Могут быть, но нет общих методов для нахождения всех минимумов / максимумов / нулей произвольных функций. (Функция может даже иметь бесконечное количество корней.) Ваша идея первого построения функции — хорошая. Тогда вам будет легче понять ее поведение.

ответ @DrV выглядит хорошо, поэтому здесь я предложу только другой подход к части (ii) вашего вопроса. Полезным способом визуализации корней сложной функции является построение 0 контуров действительной и мнимой частей. То есть вычислите z = Dm(...) на достаточно плотной сетке, а затем используйте matplotlib contour функцию ‘s для построения контуров, где z.real 0, а где z.imag ноль. Корни функции — это точки, в которых эти контуры пересекаются.

Например, вот некоторый код, который генерирует график контуров.

 import numpy as np
from scipy.special import jv, hankel1, jvp, h1vp
import matplotlib.pyplot as plt


nr = 2
m = 20

def Dm(x):
    return nr * jvp(m,nr*x,1) * hankel1(m,x) - jv(m,nr*x) * h1vp(m,x,1)


x = np.linspace(0, 40, 2500)
y = np.linspace(-15, 5, 2000)
X, Y = np.meshgrid(x, y)
z = Dm(X   1j*Y)

plt.contour(x, y, z.real, [0], colors='b', lw=0.25)
plt.contour(x, y, z.imag, [0], colors='r', lw=0.25)
plt.savefig('contours.png')
  

Вот сюжет. Каждое пересечение красной и синей линий является корнем.

График показывает, что вы не должны ожидать нахождения корней с большими отрицательными мнимыми частями. Вероятно, вы могли бы изучить асимптотическое поведение функций Бесселя и Ханкеля с большими аргументами для идей, чтобы доказать это.