Полные эллиптические интегралы в MATLAB

[Pumpov]

Добрый день!

Для расчета магнитного поля диаметрально намагниченного цилиндра возникла необходимость вычислить эллиптические интегралы 1-го и 2-го рода: $K[M]=\int\limits_{0}^{\pi/2}\frac{d\theta}{\sqrt{1-M{\sin^2}\theta}}$ и $E[M]=\int\limits_{0}^{\pi/2}\sqrt{1-M{\sin^2}\thetad\theta}$, соответственно. Я использую MATLAB.
Проблема заключается в том, что в MATLAB (и, насколько я понял, по определению этих интегралов) предусмотрен расчёт с помощью функции "ellipke()" только для значений из промежутка $[0;1]$. Но в формулах для магнитного поля требуется иногда подставлять в аргумент эллиптических интегралов и отрицательные значения. Собственно, как легко видеть из формул, эти интегралы существуют и при отрицательных $M$. Однако программа выдаёт ошибку, что "M must be in the range 0 <= M <= 1.".
Интересно, что данные эллиптические интегралы можно вычислить в MATLAB, если в качестве аргумента "x" написать, например, не "ellipke(x)", а "ellipke(sym(x))". И программа считает интеграл. Но это не подходит, если нужно создать функцию, чтобы потом подставлять в неё массив "X", созданный через "meshgrid", т.е. "ellipke(sym(X))" не срабатывает для поэлементного вычисления.
Какие могут быть варианты решения моей проблемы?

[Pumpov]

Точнее, мой вопрос такой: можно ли заставить MATLAB вычислить $K[M]$ и $E[M]$ при отрицательных значениях $M$ с помощью функции "ellipke()" сразу при подстановке массива аргументов, чтобы не прописывать самому цикл для последовательной постановки каждого значения?

[novichok2018]

Это гипергеометрическая функция Гаусса. Для неё есть формулы продолжения на все значения, формулы явные. Выразите одно через другое, продолжайте по известным формулам, и получите ответ.

[Pumpov]

novichok2018 в сообщении #1547723 писал(а):

Это гипергеометрическая функция Гаусса. Для неё есть формулы продолжения на все значения, формулы явные. Выразите одно через другое, продолжайте по известным формулам, и получите ответ.

Спасибо, я нашёл формулы на Википедии. В Матлабе, правда, нет встроенной гипергеометрической функции, но нашлись программы, которые выложены в виде скриптов на mathworks.com. Буду пробовать.

[GAA]

Pumpov в сообщении #1547722 писал(а):

можно ли заставить MATLAB вычислить $K[M]$ и $E[M]$ при отрицательных значениях $M$ с помощью функции "ellipke()" сразу при подстановке массива аргументов

Да, можно (стандартный синтаксис).

>> M= sym([-0.5 -0.1 0.3]);
>> vpa(ellipke(M))
ans = 
[ 1.4157372084259561988921659654236, 1.5335928197134568814982146172928, 1.7138894481787910620389348450438]

Или

>> vpa(ellipke(sym([-0.5 -0.1 0.3])))
ans =
[ 1.4157372084259561988921659654236, 1.5335928197134568814982146172928, 1.7138894481787910620389348450438]

[novichok2018]

МАТЕМАТИКА считает:
In[2]:= N[EllipticK[-5]]
Out[2]= 0.955504

[GAA]

Pumpov в сообщении #1547708 писал(а):

Но это не подходит, если нужно создать функцию, чтобы потом подставлять в неё массив "X", созданный через "meshgrid", т.е. "ellipke(sym(X))" не срабатывает для поэлементного вычисления.

В Matlab 2013b "срабатывает"

>> M = meshgrid(-0.1:0.1:0.1);
>> vpa(ellipke(sym(M)))
ans =
[ 1.5335928197134568814982146172928, 1.5707963267948966192313216916398, 1.6124413487202193982299163630854]
[ 1.5335928197134568814982146172928, 1.5707963267948966192313216916398, 1.6124413487202193982299163630854]
[ 1.5335928197134568814982146172928, 1.5707963267948966192313216916398, 1.6124413487202193982299163630854]

Укажите, пожалуйста, когда не "срабатывает". [Естественно символьный toolbox должен быть установлен]

[Pumpov]

GAA в сообщении #1547740 писал(а):

Pumpov в сообщении #1547722 писал(а):

можно ли заставить MATLAB вычислить $K[M]$ и $E[M]$ при отрицательных значениях $M$ с помощью функции "ellipke()" сразу при подстановке массива аргументов

Да, можно (стандартный синтаксис).

>> M= sym([-0.5 -0.1 0.3]);
>> vpa(ellipke(M))
ans = 
[ 1.4157372084259561988921659654236, 1.5335928197134568814982146172928, 1.7138894481787910620389348450438]

Или

>> vpa(ellipke(sym([-0.5 -0.1 0.3])))
ans =
[ 1.4157372084259561988921659654236, 1.5335928197134568814982146172928, 1.7138894481787910620389348450438]

Спасибо! Я вижу, через команду "vpa" таким образом вычисляется первая часть - интеграл 1-го рода от каждого значения массива. А можно сделать, чтобы и второго рода интеграл тоже выводился по принципу

[K,E] = ellipke(M)

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

[K0,E0] = ellipke(sym(M))

, а потом

K=double(K0)
E=double(E0)

-- 02.02.2022, 19:23 --

Цитата:

Укажите, пожалуйста, когда не "срабатывает". [Естественно символьный toolbox должен быть установлен]

А Вы знаете, срабатывает, перепутал - не срабатывало дальнейшее использование чисел, т.к. я забыл их конвертировать в "double".

[GAA]

Используется синтаксис Matlab M

>> M= sym([-0.5 -0.1 0.3]);
>> [K0, E0] = ellipke(M)
K0 = [ ellipticK(-1/2), ellipticK(-1/10), ellipticK(3/10)]
E0 = [ ellipticE(-1/2), ellipticE(-1/10), ellipticE(3/10)]
>> vpa(K0), vpa(E0)
ans = 
[ 1.4157372084259561988921659654236, 1.5335928197134568814982146172928, 1.7138894481787910620389348450438]
ans =
[ 1.7517712756948178620265018389845, 1.6093590249375295937017355096028, 1.4453630644126652620116176014803]

[Функция] double — это конструктор массива чисел с плавающей точкой [двойной точности].
Если нужны результаты в виде массива с плавающей точкой [двойной точности], то, конечно, можно использовать double для преобразования из sym в double.

[Pumpov]

Понятно. Большое спасибо!