реклама на сайте
подробности

 
 
 
Reply to this topicStart new topic
> Вычисление через функцию Бесселя
HardEgor
сообщение Apr 29 2017, 13:56
Сообщение #1


Профессионал
*****

Группа: Свой
Сообщений: 1 697
Регистрация: 3-03-06
Из: Tomsk
Пользователь №: 14 925



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

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

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

Интересно как более простое аналитическое решение, так и численное, типа вбить формулу в эксель/матлаб/маткад/etc и получить результат. Хорошо бы узнать книжку в которой это описывается/решается, наверняка математики всё уже придумали sm.gif
Go to the top of the page
 
+Quote Post
Guest_TSlicer_*
сообщение Apr 29 2017, 15:56
Сообщение #2





Guests






https://www.wolframalpha.com
Прикрепленное изображение


Графики для первых двух положительных корней:
(Wolfram Mathematica v.11.1 2017 г.)
Прикрепленное изображение


Код для вычисления позитивных корней:
Код
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 - Apr 29 2017, 15:54
Go to the top of the page
 
+Quote Post
HardEgor
сообщение Apr 29 2017, 16:06
Сообщение #3


Профессионал
*****

Группа: Свой
Сообщений: 1 697
Регистрация: 3-03-06
Из: Tomsk
Пользователь №: 14 925



Цитата(TSlicer @ Apr 29 2017, 22:50) *
Код для вычисления позитивных корней:
Код
FullSimplify[1 - (x BesselJ[0, x])/BesselJ[1, x] == 0, x > 0]

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

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

Фортран?
Go to the top of the page
 
+Quote Post
andrew_b
сообщение Apr 29 2017, 16:18
Сообщение #4


Профессионал
*****

Группа: Свой
Сообщений: 1 800
Регистрация: 30-12-04
Из: Воронеж
Пользователь №: 1 757



Цитата(HardEgor @ Apr 29 2017, 19:06) *
Фортран?
Видимо, вы никогда код на Фортране не видели, ибо такое не забывается.
И похоже, что Паскаль вы не видели тоже, ибо это он.
Go to the top of the page
 
+Quote Post
Guest_TSlicer_*
сообщение Apr 29 2017, 16:25
Сообщение #5





Guests






Цитата(andrew_b @ Apr 29 2017, 16:18) *
Видимо, вы никогда код на Фортране не видели, ибо такое не забывается.
И похоже, что Паскаль вы не видели тоже, ибо это он.

Да ладно, не все же программисты, тем более универсалы sm.gif
Go to the top of the page
 
+Quote Post
HardEgor
сообщение Apr 29 2017, 17:16
Сообщение #6


Профессионал
*****

Группа: Свой
Сообщений: 1 697
Регистрация: 3-03-06
Из: Tomsk
Пользователь №: 14 925



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

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

Меня смутили функции:
AP_FP_Less()
AP_FP_Greater()
BesselAsympt0()
AP_Sqr(X)
Думал что это стандартные фортрановские функции...
Go to the top of the page
 
+Quote Post
Guest_TSlicer_*
сообщение Apr 29 2017, 17:28
Сообщение #7





Guests






Цитата(HardEgor @ Apr 29 2017, 18:16) *
Меня смутили функции:
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-программ.
Go to the top of the page
 
+Quote Post

Reply to this topicStart new topic
2 чел. читают эту тему (гостей: 2, скрытых пользователей: 0)
Пользователей: 0

 


RSS Текстовая версия Сейчас: 22nd September 2017 - 06:05
Рейтинг@Mail.ru


Страница сгенерированна за 0.02197 секунд с 7
ELECTRONIX ©2004-2016