Полярный график в Matplotlib путем отображения в декартовых координатах

У меня есть переменная (P), которая является функцией угла (тета):

введите здесь описание изображения

В этом уравнении K — константа, theta_p равна нулю, а I — модифицированная функция Бесселя первого рода (порядок 0). который определяется как:

введите здесь описание изображения

Теперь я хочу построить график зависимости P от тета для разных значений константы K. Сначала я рассчитал параметр I, а затем подставил его в первое уравнение, чтобы рассчитать P для различных углов тета. Я сопоставил его с декартовой координатой, поставив:

х = P * cos (тета)

у = P * грех (тета)

Вот моя реализация Python с использованием matplotlib и scipy, когда константа k = 2.0:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

def integrand(x, a, k):
   return a*np.exp(k*np.cos(x))

theta = (np.arange(0, 362, 2))
theta_p = 0.0

X = []
Y = []

for i in range(len(theta)):
    a = (1 / np.pi)
    k = 2.0
    Bessel = quad(integrand, 0, np.pi, args=(a, k))
    I = list(Bessel)[0]
    P = (1 / (np.pi * I)) * np.exp(k * np.cos(2 * (theta[i]*np.pi/180. - theta_p)))
    x = P*np.cos(theta[i]*np.pi/180.)
    y = P*np.sin(theta[i]*np.pi/180.)
    X.append(x)
    Y.append(y)

plt.plot(X,Y,  linestyle='-', linewidth=3, color='red')
axes = plt.gca()
plt.show()

Я должен получить набор графиков, как показано на рисунке ниже, для разных значений K: введите здесь описание изображения

(Обратите внимание, что распределения были нанесены на круг единицы 1 для облегчения визуализации)

Однако кажется, что графики, созданные приведенным выше кодом, не похожи на приведенный выше рисунок. Есть идеи, в чем проблема с приведенной выше реализацией? Заранее спасибо за помощь.

Вот как это выглядит (для k=2): введите здесь описание изображения

Эталоном для этих формул являются уравнения 5 и 6, которые вы можете найти здесь< /а>


person Leo    schedule 24.12.2019    source источник
comment
графики, созданные приведенным выше кодом, не похожи на приведенный выше рисунок - как они выглядят?   -  person ForceBru    schedule 24.12.2019
comment
Почему бы не использовать функцию Бесселя, определенную в scipy?: from scipy.special import i0   -  person s.ouchene    schedule 24.12.2019
comment
Не могли бы вы добавить ссылку, откуда вы взяли формулы?   -  person s.ouchene    schedule 24.12.2019
comment
Я отредактировал свой пост и ответил на ваши вопросы   -  person Leo    schedule 25.12.2019


Ответы (1)


У вас ошибка в формуле.

Ваша формула дает дельту вашей функции над единичным кругом. Итак, в вашей функции, чтобы получить нужный сюжет, просто добавьте к нему 1.

Вот то, что вы хотите, с некоторым приведенным в порядок питоном. ... обратите внимание, что вы можете выполнить весь расчет значений «P» в виде пустой векторной линии, вам не нужно перебирать индексы. ... также вы можете просто сделать полярный график прямо в matplotlib - вам не нужно преобразовывать его в декартов.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

theta = np.arange(0, 2*np.pi+0.1, 2*np.pi/100)

def integrand(x, a, k):
   return a*np.exp(k*np.cos(x))

for k in np.arange(0, 5, 0.5):
    a = (1 / np.pi)
    Bessel = quad(integrand, 0, np.pi, args=(a, k))
    I = Bessel[0]    
    P = 1 + (1/(np.pi * I)) * np.exp(k * np.cos(2 * theta))
    plt.polar(theta, P)

plt.show()

Полярный график

person Richard    schedule 25.12.2019