Естественный логарифм функции Бесселя, переполнение

Я пытаюсь вычислить логарифм модифицированной функции Бесселя второго типа в MATLAB, т.е. что-то вроде этого:

log(besselk(nu, Z))

где, например,

nu = 750;
Z = 1;

У меня проблема, потому что значение log(besselk(nu, Z)) переходит в бесконечность, потому что besselk(nu, Z) - бесконечность. Однако log(besselk(nu, Z)) должно быть действительно маленьким.

Я пытаюсь написать что-то вроде

f = ******(sym('ln(besselk(******(nu), ******(Z)))'));

Однако я получаю следующую ошибку:

Ошибка при использовании mupadmex Ошибка в команде MuPAD: DOUBLE не может преобразовать входное выражение в двойной массив. Если входное выражение содержит символическую переменную, вместо этого используйте функцию VPA.

Ошибка в sym/****** (строка 514) Xstr = mupadmex ('symobj:: ******', S.s, 0) `;

Как избежать этой ошибки?

3 ответа

Как указывал ньюффа, DLMF дает асимптотические разложения K_nu (z) при больших nu. Из 10.41.2 находим для вещественных положительных аргументов z:

besselk(nu,z) ~ sqrt(pi/(2nu)) (e z/(2nu))^-nu

который дает после некоторого упрощения

log( besselk(nu,z) ) ~ 1/2*log(pi) + (nu-1/2)*log(2nu) - nu(1 + log(z))

Итак, это O (nu log (nu)). Неудивительно, что прямой расчет не выполняется при nu > 750.

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

EDIT: Я просто пробовал для nu = 750 и z = 1: приведенное выше приближение дает 4.7318e + 03, а с результатом horchler получаем log (1.02 * 10 ^ 2055) = 2055 * log (10) + log (1.02) = 4.7318e + 03. Поэтому это правильно, по крайней мере, 5 значащих цифр, для nu >= 750 и z = 1! Если это будет достаточно для вас, это будет намного быстрее, чем символическая математика.


Вы делаете несколько вещей неправильно. Нет смысла использовать ****** для ваших двух аргументов в besselk и преобразовать вывод в символический. Вы также должны избегать ввода старой строки на sym. Вместо этого вы хотите условно оценить besselk (который вернет около 1.02 × 10 2055 намного больше realmax), берем журнал результата символически, а затем преобразуем обратно в двойную точность.

Достаточно следующего: если один или несколько входных аргументов являются символьными, будет использоваться символическая версия besselk:

f = ******(log(besselk(sym(750), sym(1))))

или в старой строковой форме:

f = ******(sym('log(besselk(750, 1))'))

Если вы хотите сохранить свои параметры символическими и оценить в более позднее время:

syms nu Z;
f = log(besselk(nu, Z))
******(subs(f, {nu, Z}, {750, 1}))

Убедитесь, что вы не перевернули значения nu и Z в своей математике, поскольку большие ордера (nu) не очень распространены.


Вы пробовали интегральное представление?

Журнал [Интеграция [Cosh [Nu t]/E ^ (Z Cosh [t]), {t, 0, бесконечность}]]

licensed under cc by-sa 3.0 with attribution.