Интеграл функции интенсивности в питоне

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

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

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

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

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

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

or

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. Величину ошибки можно оценить, попробовав разные значения начальной точки. Результат интеграции очень мало меняется между 1e-9 и 1e-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()

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

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

Что касается последней части вопроса: I0 определяет максимальную интенсивность (единицы могут быть, например, Вт/м2), тогда как интеграл дает общую мощность (если интенсивность в Вт/м2, общая мощность в Вт) . Установка максимальной интенсивности на 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).

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

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

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

331.4990321315221