Перейти к содержанию

Вычисление через функцию Бесселя

Надо найти первый положительный корень уравнения:

x*J0(x)/J1(x) = 1-y , где y известен, Jx() - функции Бесселя, найти x

Прямой путь - разложение J0 и J1 в ряд и вычисление, как-то сложновато выходит.
Нутром чую, что есть более простые способы, подскажите пожалуйста, а то как-то всё подзабылось....

Интересно как более простое аналитическое решение, так и численное, типа вбить формулу в эксель/матлаб/маткад/etc и получить результат. Хорошо бы узнать книжку в которой это описывается/решается, наверняка математики всё уже придумали sm.gif

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты
Гость TSlicer
https://www.wolframalpha.com
[attachment=106913:roots.png]

Графики для первых двух положительных корней:
(Wolfram Mathematica v.11.1 2017 г.)
[attachment=106915:roots_2.png]

Код для вычисления позитивных корней:
Код
FullSimplify[1 - (x BesselJ[0, x])/BesselJ[1, x] == 0, x > 0]


Можно и самому через ряды Чебышева sm.gif

Код
function BesselJ0(X : Double):Double;
var
    XSq : Double;
    NN : Double;
    PZero : Double;
    QZero : Double;
    P1 : Double;
    Q1 : Double;
begin
    if AP_FP_Less(X,0) then
    begin
        X := -X;
    end;
    if AP_FP_Greater(X,8.0) then
    begin
        BesselAsympt0(X, PZero, QZero);
        NN := X-Pi/4;
        Result := Sqrt(2/Pi/X)*(PZero*Cos(NN)-QZero*Sin(NN));
        Exit;
    end;
    XSq := AP_Sqr(X);
    P1 := 26857.86856980014981415848441;
    P1 := -40504123.71833132706360663322+XSq*P1;
    P1 := 25071582855.36881945555156435+XSq*P1;
    P1 := -8085222034853.793871199468171+XSq*P1;
    P1 := 1434354939140344.111664316553+XSq*P1;
    P1 := -136762035308817138.6865416609+XSq*P1;
    P1 := 6382059341072356562.289432465+XSq*P1;
    P1 := -117915762910761053603.8440800+XSq*P1;
    P1 := 493378725179413356181.6813446+XSq*P1;
    Q1 := 1.0;
    Q1 := 1363.063652328970604442810507+XSq*Q1;
    Q1 := 1114636.098462985378182402543+XSq*Q1;
    Q1 := 669998767.2982239671814028660+XSq*Q1;
    Q1 := 312304311494.1213172572469442+XSq*Q1;
    Q1 := 112775673967979.8507056031594+XSq*Q1;
    Q1 := 30246356167094626.98627330784+XSq*Q1;
    Q1 := 5428918384092285160.200195092+XSq*Q1;
    Q1 := 493378725179413356211.3278438+XSq*Q1;
    Result := P1/Q1;
end;
Изменено пользователем TSlicer

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты
Цитата(TSlicer @ Apr 29 2017, 22:50) <{POST_SNAPBACK}>
Код для вычисления позитивных корней:
Код
FullSimplify[1 - (x BesselJ[0, x])/BesselJ[1, x] == 0, x > 0]

Спасибо большое, а то я вписал общую формулу, а он мне "Standard computation time exceeded..." и денежку запросил.

Цитата(TSlicer @ Apr 29 2017, 22:56) <{POST_SNAPBACK}>
Можно и самому через ряды Чебышева sm.gif

Фортран?

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты
Цитата(HardEgor @ Apr 29 2017, 19:06) <{POST_SNAPBACK}>
Фортран?
Видимо, вы никогда код на Фортране не видели, ибо такое не забывается.
И похоже, что Паскаль вы не видели тоже, ибо это он.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты
Гость TSlicer
Цитата(andrew_b @ Apr 29 2017, 16:18) <{POST_SNAPBACK}>
Видимо, вы никогда код на Фортране не видели, ибо такое не забывается.
И похоже, что Паскаль вы не видели тоже, ибо это он.

Да ладно, не все же программисты, тем более универсалы sm.gif

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты
Цитата(andrew_b @ Apr 29 2017, 23:18) <{POST_SNAPBACK}>
Видимо, вы никогда код на Фортране не видели, ибо такое не забывается.
И похоже, что Паскаль вы не видели тоже, ибо это он.

На паскале писал, просто я писал на всех языках за последние 30 лет sm.gif и они уже как-то перепутались у меня в голове wink.gif

Меня смутили функции:
AP_FP_Less()
AP_FP_Greater()
BesselAsympt0()
AP_Sqr(X)
Думал что это стандартные фортрановские функции...

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты
Гость TSlicer
Цитата(HardEgor @ Apr 29 2017, 18:16) <{POST_SNAPBACK}>
Меня смутили функции:
AP_FP_Less()
AP_FP_Greater()
BesselAsympt0()
AP_Sqr(X)
Думал что это стандартные фортрановские функции...

Это частные локальные функции, эквивалентные стандартным в Pascal (Delphi);
<
>
sqrt

Вариант вычисления Бесселя и работает на асимптоте:
procedure BesselAsympt0(X : Double; var PZero : Double; var QZero : Double);

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Для публикации сообщений создайте учётную запись или авторизуйтесь

Вы должны быть пользователем, чтобы оставить комментарий

Создать учетную запись

Зарегистрируйте новую учётную запись в нашем сообществе. Это очень просто!

Регистрация нового пользователя

Войти

Уже есть аккаунт? Войти в систему.

Войти
Авторизация