Интеграл функции Intensity в python

Проблема, похоже, связана с поведением функции вблизи нуля. Если функция построена, она выглядит гладко:

Однако scipy.integrate.quad жалуется на ошибки округления, что очень странно с этой красивой кривой. Однако функция не определена в 0 (разумеется, вы делите на ноль!), Поэтому интеграция идет не так хорошо.

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

Однако, я думаю, у меня есть догадка о вашей проблеме. Насколько я помню, интеграл, который вы показали, представляет собой фактически интенсивность (мощность / площадь) дифракции Фраунгофера как функцию расстояния от центра. Если вы хотите интегрировать полную мощность в радиусе, вам придется сделать это в двух измерениях.

По простым правилам интеграции области вы должны умножить свою функцию на 2 pi r перед тем, как интегрировать (или x вместо r в вашем случае). Тогда это становится:

 f = lambda(r): r*(sp.j1(r)/r)**2 

или

 f = lambda(r): sp.j1(r)**2/r 

или даже лучше:

 f = lambda(r): r * (sp.j0(r) + sp.jn(2,r)) 

Последняя форма лучше, поскольку она не страдает от каких-либо особенностей. Он основан на комментарии Хайме к первоначальному ответу (см. Комментарий ниже этого ответа!).

(Заметим, что я опустил пару констант.) Теперь вы можете интегрировать его с нуля на бесконечность (без отрицательных радиусов):

 fullpower = quad(f, 1e-9, np.inf)[0] 

Затем вы можете интегрироваться из какого-то другого радиуса и нормализовать по полной интенсивности:

 pwr = quad(f, 1e-9, 3.8317)[0] / fullpower 

И вы получаете 0,839 (что близко к 84%). Если вы попробуете более дальний радиус (13.33):

 pwr = quad(f, 1e-9, 13.33) 

что дает 0,954.

Следует отметить, что мы вводим небольшую ошибку, начиная интеграцию из 1e-9 вместо 0. Величину ошибки можно оценить, попробовав разные значения для начальной точки. Результат интеграции очень мало меняется между 1е-9 и 1е-12, поэтому они кажутся безопасными. Конечно, вы могли бы использовать, например, 1е-30, но тогда в делении может быть численная нестабильность. (В этом случае нет, но в целом особенности являются численно злыми).

Давайте сделаем еще одно:

 import matplotlib.pyplot as plt import numpy as np x = linspace(0.01, 20, 1000) intg = np.array([ quad(f, 1e-9, xx)[0] for xx in x]) plt.plot(x, intg/fullpower) plt.grid('on') plt.show() 

И это то, что мы получаем:

По крайней мере, это выглядит правильно, темные полосы диска Airy хорошо видны.

Что касается последней части вопроса: I0 определяет максимальную интенсивность (единицы измерения могут быть, например, W / m2), тогда как интеграл дает полную мощность (если интенсивность находится в Вт / м2, полная мощность в W) , Установка максимальной интенсивности до 100 не гарантирует ничего об общей мощности. Вот почему важно рассчитать общую мощность.

Фактически существует замкнутое уравнение формы для полной мощности, излучаемой на круговой области:

P ( x ) = P0 (1 – J0 ( x ) ^ 2 – J1 ( x ) ^ 2),

где P0 – полная мощность.

Обратите внимание, что вы также можете получить закрытое решение для вашей интеграции с помощью Sympy :

 import sympy as sy sy.init_printing() # LaTeX like pretty printing in IPython x,d = sy.symbols("x,d", real=True) I0=100 dist=3.8317 f = I0*((2*sy.besselj(1,x)/x)**2) # the integrand F = f.integrate((x, -d, d)) # symbolic integration print(F.evalf(subs={d:dist})) # numeric evalution 

F оценивает:

 1600*d*besselj(0, Abs(d))**2/3 + 1600*d*besselj(1, Abs(d))**2/3 - 800*besselj(1, Abs(d))**2/(3*d) 

с besselj(0,r) соответствующей sp.j0(r) .

Они могут быть сингулярностью в алгоритме интеграции при выполнении jacobian при x = 0. Вы можете исключить эти точки из интеграции с « точками »:

 f = lambda x:( I0*((2*sp.j1(x)/x)**2)) I = quad(f, -dist, dist, points = [0]) 

Я получаю следующий результат (это ваш желаемый результат?)

 331.4990321315221