Алгоритм быстрого интегрирования функции Планка

[sciko]

Нужен сабж. Кто не знает, что это такое, может почитать в рукипедии.

Интегрирование не обязательно должно быть точным. Если есть, то возьму и обратный алгоритм. Если есть быстрые алгоритмы интегрирования уравнений Сакумы-Хаттори, то сойдут и они.

[Crazy]

Квадратурная формула Гаусса?

[sciko]

Квадратурная формула Гаусса?

Проблема в том, что для неё функцию раскладывают в ряд полиномов. Но в формуле Планка стоит экспонента, что потребует высоких степеней полиномов. А на высоких степенях квадратурная формула Гаусса будет не только увеличится время работы (O(N^2)), но и падать точность.

[Crazy]

В квадратурной формуле Гаусса узлами интерполяции являются корни полинома Лежандра.
Там же нужны только узлы и веса, которые известны.

[Фантом]

В принципе, ответ зависит от ожидаемого количества интегрирований и платформы. Если прогонов будет много, а существенных ограничений памяти нет, то проще всего сделать так: обезразмерить функцию, перейдя к переменной h\nu/kT, потом один раз теми же гауссовыми квадратурами (никакого разложения в ряд там не будет, так что это не проблема) насчитать таблицу интегралов, например, от 0 до разных верхних пределов в интересующей Вас спектральной области (по-видимому, Вам нужны окрестности максимума, иначе можно просто воспользоваться приближениями Вина или Рэлея-Джинса, они интегрируются аналитически, так что таблица будет не слишком большой), а потом просто считать разности двух чисел из таблицы по мере надобности и домножать их на размерный коэффициент.

Второй вариант, если память жалко: ту же обезразмеренную функцию разложить в ряд Тейлора в окрестности 3 (это близко к максимуму, хотя можно взять и точнее, где-то примерно 2.8), оставить любое разумное число членов (хоть пару десятков) и проинтегрировать сумму аналитически (что проблем не составляет). Дальше все очевидно.

[sciko]

В принципе, ответ зависит от ожидаемого количества интегрирований и платформы.

Около миллиона вычислений уравнения вида:

Код: Выделить всё

R(T[0])=tay*(eps*R(T[1])+(1-eps)*R(T[2]))+(1-tay)*R(T[3])

Причём R(T[3]) тоже вычисляется не напрямую.

Платформа i686 (мне это нужно для демонстрации).

Если прогонов будет много, а существенных ограничений памяти нет, то проще всего сделать так: обезразмерить функцию, перейдя к переменной h\nu/kT, потом один раз теми же гауссовыми квадратурами (никакого разложения в ряд там не будет, так что это не проблема) насчитать таблицу интегралов, например, от 0 до разных верхних пределов в интересующей Вас спектральной области (по-видимому, Вам нужны окрестности максимума, иначе можно просто воспользоваться приближениями Вина или Рэлея-Джинса, они интегрируются аналитически, так что таблица будет не слишком большой), а потом просто считать разности двух чисел из таблицы по мере надобности и домножать их на размерный коэффициент.

Не пойдёт. Пределы интегрирования не меняются: Вы не поверите, но у регистрирующей аппаратуры спектр чувствительности всегда один и тот же.

ту же обезразмеренную функцию разложить в ряд Тейлора в окрестности 3 (это близко к максимуму, хотя можно взять и точнее, где-то примерно 2.8), оставить любое разумное число членов (хоть пару десятков) и проинтегрировать сумму аналитически (что проблем не составляет).

После подстановки

Код: Выделить всё

z=c2/(k*T)

получается выражение вида

Код: Выделить всё

R(T)=c0*(z^3)/(exp(z)-1)

, которая весьма непросто интегрируется даже будучи разложена в ряд Тейлора.

[sciko]

Нашёл в интернете вариант аппроксимации функции вида

Код: Выделить всё

F(z)=1/(exp(z)-1)

как

Код: Выделить всё

\sum\limits_{m=1}^{t-1} \exp(-mz) + \frac{\exp(-(t-0.5)z)}{z}

[Фантом]

Не пойдёт. Пределы интегрирования не меняются: Вы не поверите, но у регистрирующей аппаратуры спектр чувствительности всегда один и тот же.

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

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

После подстановки
получается выражение вида
, которая весьма непросто интегрируется даже будучи разложена в ряд Тейлора.

А Вы попробовали? Похоже, что нет (напоминаю, что речь идет о разложении в ряд Тейлора в окрестности z \approx 2.8, а не в ряд Маклорена). Попробуйте - обнаружите, что абсолютные значения коэффициентов разложения первообразной уменьшаются примерно на порядок при увеличении номера члена на единицу (да и ряд не знакопостоянный). Если работать приходится вдалеке от максимума, то сходимость только улучшится (впрочем, для высокоэнергетичной области проще сразу взять приближение Вина)ю

[sciko]

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

Очень даже ограничивается. Любой численный расчёт (кроме может экзотического моделирования) основан на данных регистрирующей аппаратуры.

при интегрировании безразмерной функции при одном и том же диапазоне частот (или длин волн) пределы интегрирования будут меняться - они зависят от температуры.

Согласен, ступил.

напоминаю, что речь идет о разложении в ряд Тейлора в окрестности z \approx 2.8, а не в ряд Маклорена

Разложение в ряд Маклорена функции вида \frac{1}{\exp{z}-1} очень нетривиальная задача, т.к. связана с деление на выражение \exp{z}-1, которое при z=0 превращается в 0.

А Вы попробовали?

Вот в этом-то и основная проблема.

Я не поленился и вычислил аналитически 5 производную от функции \frac{z^3}{\exp{z}-1}:

-\frac{120\,{z}^{3}\,{e}^{5\,z}}{{\left( {e}^{z}-1\right) }^{6}}+\frac{240\,{z}^{3}\,{e}^{4\,z}}{{\left( {e}^{z}-1\right) }^{5}}+\frac{360\,{z}^{2}\,{e}^{4\,z}}{{\left( {e}^{z}-1\right) }^{5}}-\frac{150\,{z}^{3}\,{e}^{3\,z}}{{\left( {e}^{z}-1\right) }^{4}}-\frac{540\,{z}^{2}\,{e}^{3\,z}}{{\left( {e}^{z}-1\right) }^{4}}-\frac{360\,z\,{e}^{3\,z}}{{\left( {e}^{z}-1\right) }^{4}}+\frac{30\,{z}^{3}\,{e}^{2\,z}}{{\left( {e}^{z}-1\right) }^{3}}+\frac{210\,{z}^{2}\,{e}^{2\,z}}{{\left( {e}^{z}-1\right) }^{3}}+\frac{360\,z\,{e}^{2\,z}}{{\left( {e}^{z}-1\right) }^{3}}+\frac{120\,{e}^{2\,z}}{{\left( {e}^{z}-1\right) }^{3}}-\frac{{z}^{3}\,{e}^{z}}{{\left( {e}^{z}-1\right) }^{2}}-\frac{15\,{z}^{2}\,{e}^{z}}{{\left( {e}^{z}-1\right) }^{2}}-\frac{60\,z\,{e}^{z}}{{\left( {e}^{z}-1\right) }^{2}}-\frac{60\,{e}^{z}}{{\left( {e}^{z}-1\right) }^{2}}

Его не просто записать без ошибки. Вопрос вычислительной ошибки (хотя бы из-за большого числа вычисления экспоненты) тоже остаётся открытым.

Во-вторых, я посчитал значения коэффициентов для z=2.8:

Код: Выделить всё

F(z)=1.421333898103823+0.0094962419063318*(x-2.8)-0.22291408433486*(x-2.8)^2+0.04481206762846*(x-2.8)^3+0.0062181477546315*(x-2.8)^4-0.0039522996669334*(x-2.8)^5+...

Такое разложение даёт погрешность аппроксимации около 0,1% без учёта погрешностей вычисления.

Кроме того, видно, что

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

малость не соответствует действительности.

[Фантом]

Очень даже ограничивается. Любой численный расчёт (кроме может экзотического моделирования) основан на данных регистрирующей аппаратуры.

Ну, например, для меня такое моделирование экзотическим не является.

Я не поленился и вычислил аналитически 5 производную от функции \frac{z^3}{\exp{z}-1}:

Его не просто записать без ошибки. Вопрос вычислительной ошибки (хотя бы из-за большого числа вычисления экспоненты) тоже остаётся открытым.

Бог мой, зачем так-то?! Берете Maxima, вбиваете туда что-нибудь вроде

Код: Выделить всё

float(taylor(z^3/(exp(z)-1),z,3,20));

(это разложение в окрестности 3 до 20-го порядка) и получаете на выходе

Код: Выделить всё

-9.913780081248967E-16*(z-3.0)^20
+3.411722889745232E-15*(z-3.0)^19
+2.758965562002768E-14*(z-3.0)^18
-3.309048653442929E-13*(z-3.0)^17
+6.478314318465098E-13*(z-3.0)^16
+1.215065260779206E-11*(z-3.0)^15
-1.042689019239971E-10*(z-3.0)^14
+3.700605704183545E-11*(z-3.0)^13
+4.823181390211337E-9*(z-3.0)^12
-3.074667568123652E-8*(z-3.0)^11
-4.762271766256543E-8*(z-3.0)^10
+1.7671934930080643E-6*(z-3.0)^9
-8.533959430475123E-6*(z-3.0)^8
-3.0966952625438154E-5*(z-3.0)^7
+6.162861409090619E-4*(z-3.0)^6
-.003190926420188421*(z-3.0)^5
+.002649806463944675*(z-3.0)^4
+0.0483084608109771*(z-3.0)^3
-.1948352183667711*(z-3.0)^2
-0.0741233432917029*(z-3.0)
+1.414683805263911

Работы на несколько секунд, итоговый результат, при желании, очень легко проверить.

малость не соответствует действительности.

Я вообще-то писал про первообразную (которая, кстати, получается также однострочной командой в той же Maxima). Но даже если брать разложение самой функции Планка, то до 5-го порядка это естественно (или Вы надеялись, что в окрестности максимума функции для первых членов разложения будет иначе?), а дальше картина будет примерно такой, как я описывал (собственно, см.выше). Просто, честно сказать, я не предполагал, что Вы будете считать это вручную (и, соответственно, с таким малым числом членов).