Функции Бесселя в программе символьной математики SymPy / Хабр

Введение:
Большое число самых разнообразных задач, относящихся практически ко всем важнейшим разделам математической физики и призванных ответить на актуальные технические вопросы, связано с применением функций Бесселя.

Функции Бесселя широко используются при решении задач акустики, радиофизики, гидродинамики, задач атомной и ядерной физики. Многочисленные приложения функций Бесселя к теории теплопроводности и теории упругости (задачи о колебаниях пластинок, задачи теории оболочек, задачи определения концентрации напряжения вблизи трещин).

Такая популярность функций Бесселя объясняется тем, что решение уравнений математической физики, содержащих оператор Лапласа в цилиндрических координатах, классическим методом разделения переменных приводит к обыкновенному дифференциальному уравнению, служащему для определения этих функций[1].

Функции Бесселя названы по имени немецкого астронома Фридриха Бесселя, который в работе 1824 года, изучая движение планет вокруг солнца, вывел рекуррентные соотношения для функций Бесселя $J_{v}(x)$ , получил для целых $inline$ интегральное представление функции $J_{v}(x)$ , доказал наличие бесчисленного множества нулей функции $J_{0}(x)$ и составил первые таблицы для функций $J_{1}(x)$ и $J_{2}(x)$ .

Однако, впервые одна из функций Бесселя $J_{0}(x)$ была рассмотрена еще в 1732 году Даниилом Бернулли в работе, посвященной колебанию тяжелых цепей. Д. Бернулли нашел выражение функции $J_{0}(x)$ в виде степенного ряда и заметил (без доказательства), что уравнение $J_{0}(x)=0$ имеет бесчисленное множество действительных корней.

Следующей работой, в которой встречаются функции Бесселя, была работа Леонардо Эйлера 1738 года, посвященная изучению колебаний круглой мембраны. В этой работе Л. Эйлер нашел для целых $inline$ выражение функции Бесселя $J_{v}(x)$ в виде ряда по степеням $inline$ , а в последующих работах распространил это выражение на случай произвольных значений индекса $inline$ . Кроме того, Л. Эйлер доказал, что для $inline$ , равного целому числу с половиной, функции $J_{v}(x)$ выражаются через элементарные функции.

Он заметил (без доказательства), что при действительных $inline$ функции $J_{v}(x)$ имеют бесчисленное множество действительных нулей и дал интегральное представление для $J_{v}(x)$ . Некоторые исследователи считают, что основные результаты, связанные с функциями Бесселя и их приложениями в математической физике, связаны с именем Л. Эйлера.

Изучить свойство функций Бесселя и одновременно освоить методы решения уравнений, сводящихся к функциям Бесселя, позволяет свободно распространяемая программа символьной математики SymPy — библиотеки Python.

В программе символьной математики SymPy графики функций Бесселя первого рода целых порядков можно построить, пользуясь соотношением для суммы ряда:

$J_{p}(x)=\sum_{m=0}^{\infty}\frac{x^{2m+p}(-1)^{m}}{2^{2m+p}m!\Gamma (p+m+1)}$

Функции Бесселя первого рода

from sympy import*
from sympy.plotting import plot
x,n, p=var('x,n, p')
def besselj(p,x):
    return summation(((-1)**n*x**(2*n+p))/(factorial(n)*gamma(n+p+1)*2**(2*n+p)),[n,0,oo])
st="J_{p}(x)"
p1=plot(besselj(0,x),(x,-20,20),line_color='b',title=' $'+st+ '$',show=False)    
p2=plot(besselj(1,x),(x,-20,20),line_color='g',show=False)  
p3=plot(besselj(2,x),(x,-20,20),line_color='r',show=False)
p4=plot(besselj(3,x),(x,-20,20),line_color='c',show=False)
p1.extend(p2)
p1.extend(p3)
p1.extend(p4)
p1.show()

При помощи соотношения для суммы ряда можно доказать свойство этих функций для целых порядков

$J_{1}(x)=-J_{-1}(x):$

Свойство функции Бесселя первого рода

from sympy import*
from sympy.plotting import plot
x,n, p=var('x,n, p')
def besselj(p,x):
    return summation(((-1)**n*x**(2*n+p))/(factorial(n)*gamma(n+p+1)*2**(2*n+p)),[n,0,oo])
st="J_{1}(x)=-J_{-1}(x)"
p1=plot(besselj(1,x),(x,-10,10),line_color='b',title=' $'+st+ '$',show=False)    
p2=plot(besselj(-1,x),(x,-10,10),line_color='r',show=False)  
p1.extend(p2)
p1.show()

Для демонстрации условий Коши, построим функцию $J_{1/3}(x)$ и её производную $\frac{dJ_{1/3}(x)}{dx}:$ :

Функция дробного порядка и её производная

from sympy import*
from sympy.plotting import plot
x,n, p=var('x,n, p')
def besselj(p,x):
    return summation(((-1)**n*x**(2*n+p))/(factorial(n)*gamma(n+p+1)*2**(2*n+p)),[n,0,oo])
st="J_{1/3}(x),J{}'_{1/3}(x)"
p1=plot(besselj(1/3,x),(x,-1,10),line_color='b',title=' $'+st+ '$',ylim=(-1,2),show=False)
def dbesselj(p,x):
    return diff(summation(((-1)**n*x**(2*n+p))/(factorial(n)*gamma(n+p+1)*2**(2*n+p)),[n,0,oo]),x)
p2=plot(dbesselj(1/3,x),(x,-1,10),line_color='g',show=False)  
p1.extend(p2)
p1.show()

Однако, для практических расчётов используется замечательный модуль mpmath, позволяющий численно не только решать уравнения с функциями Бесселя первого и второго рода в том числе и модифицированные всех допустимых порядков, но и строить графики с автоматическим масштабированием.

Кроме того, модуль mpmath не требует специальных средств для совместного использования символьной и численной математики. Историю создания этого модуля и возможности его использования для обратного преобразования Лапласа я уже рассматривал в публикации [2]. Теперь продолжим рассмотрение mpmath для работы с функциями Бесселя [3].

Функция Бесселя первого рода $J_{N}(x)$
mpmath.besselj(n, x, derivative=0) — дает функцию Бесселя первого рода $inline$ . Функции $J_{N}(x)$ является решением следующего дифференциального уравнения:

$x^{2}{y}''+x{y}'+(x^{2}-n^{2})y=0$

Для целых положительных $inline$ ведёт себя как синус или косинус, умноженный на коэффициент, который медленно убывает при $x\rightarrow \pm \infty$
Функция Бесселя первого рода $J_{N}(x)$ является частным случаем гипергеометрической функции ${o}F_{1}$ :

$J_{n}(x)=\frac{x^{n}}{2^{n}\Gamma (n+1)} {o}F_{1} (n+1,\frac{x^{2}}{4})$

Функцию Бесселя можно дифференцировать $inline$ раз при условии, что m-я производная не равна нулю:

$\frac{d^{m}}{dx^{m}}J_{n}(x)$

Функция Бесселя первого рода $J_{N}(x)$ для положительных целых порядков n = 0,1,2,3 — решение уравнения Бесселя:


from mpmath import*
j0 = lambda x: besselj(0,x)
j1 = lambda x: besselj(1,x)
j2 = lambda x: besselj(2,x)
j3 = lambda x: besselj(3,x)
plot([j0,j1,j2,j3],[0,14]

Функция Бесселя первого рода $J_{N}(x)$ в комплексной плоскости:

from sympy import*
from mpmath import*
cplot(lambda z: besselj(1,z), [-8,8], [-8,8], points=50000)

Примеры:
Функция $inline$ обеспечивает результат с заданным числом цифр $inline$ после запятой:


from mpmath import*
mp.dps = 15; mp.pretty = True
print(besselj(2, 1000))
nprint(besselj(4, 0.75))
nprint(besselj(2, 1000j))
mp.dps = 25
nprint( besselj(0.75j, 3+4j))
mp.dps = 50
nprint( besselj(1, pi))

Аргумент функции может быть большим числом:


from mpmath import*
mp.dps = 25
nprint( besselj(0, 10000))
nprint(besselj(0, 10**10))
nprint(besselj(2, 10**100))
nprint( besselj(2, 10**5*j))

Функции Бесселя первого рода удовлетворяют простым симметриям относительно $inline$ :

from sympy import*
from mpmath import*
mp.dps = 15
nprint([besselj(n,0) for n in range(5)])
nprint([besselj(n,pi) for n in range(5)])
nprint([besselj(n,-pi) for n in range(5)])

Корни не периодические, но расстояние между последовательными корнями асимптотически приближается к $inline$ . Функции Бесселя первого рода имеют следующий код:


from mpmath import*
print(quadosc(j0, [0, inf], period=2*pi))
print(quadosc(j1, [0, inf], period=2*pi))

Для $inline$ или $inline$ функция Бесселя сводится к тригонометрической функции:

from sympy import*
from mpmath import*
x = 10
print(besselj(0.5, x))
print(sqrt(2/(pi*x))*sin(x))
print(besselj(-0.5, x))
print(sqrt(2/(pi*x))*cos(x))

Могут быть вычислены производные любого порядка, отрицательные порядки соответствуют интегрированию:


from mpmath import*
mp.dps = 25
print(besselj(0, 7.5, 1))
print(diff(lambda x: besselj(0,x), 7.5))
print(besselj(0, 7.5, 10))
print(diff(lambda x: besselj(0,x), 7.5, 10))
print(besselj(0,7.5,-1) - besselj(0,3.5,-1))
print(quad(j0, [3.5, 7.5]))

Дифференцирование с нецелым порядком дает дробную производную в смысле дифференциального интеграла Римана-Лиувилля, вычисляемую с помощью функции $inline$ :


from mpmath import*
mp.dps = 15
print(besselj(1, 3.5, 0.75))
print(differint(lambda x: besselj(1, x), 3.5, 0.75))

Другие способы вызова функции Бесселя первого рода нулевого и первого порядков
mpmath.j0(x) — Вычисляет функцию Бесселя $J_{0}(x)$ ;
mpmath.j1(x) — Вычисляет функцию Бесселя $J_{1}(x)$ ;

Функции Бесселя второго рода
bessely(n, x, derivative=0) Вычисляет функцию Бесселя второго рода по соотношению:

$Y_{n}(x)=\frac{J_{n}(x)cos(\pi\cdot n)-J_{-n}(x) }{sin(\pi\cdot n)}$

Для целого числа $inline$ следующую формулу следует понимать как предел. Функцию Бесселя можно дифференцировать $inline$ раз при условии, что m-я производная не равна нулю:
$\frac{d^{m}}{dx^{m}} Y_{n}(x)$
Функция Бесселя второго рода $Y_{n}(x)$ для целых положительных порядков $inline$ .

from sympy import*
from mpmath import*
y0 = lambda x: bessely(0,x)
y1 = lambda x: bessely(1,x)
y2 = lambda x: bessely(2,x)
y3 = lambda x: bessely(3,x)
plot([y0,y1,y2,y3],[0,10],[-4,1])

Функция Бесселя 2-го рода $Y_{n}(x)$ в комплексной плоскости

from sympy import*
from mpmath import*
cplot(lambda z: bessely(1,z), [-8,8], [-8,8], points=50000)

Примеры:
Некоторые значения функции $Y_{n}(x)$ :

from sympy import*
from mpmath import*
mp.dps = 25; mp.pretty = True
print(bessely(0,0))
print(bessely(1,0))
print(bessely(2,0))
print(bessely(1, pi))
print(bessely(0.5, 3+4j))

Аргументы могут быть большими:

from sympy import*
from mpmath import*
mp.dps = 25; mp.pretty = True
print(bessely(0, 10000))
print(bessely(2.5, 10**50))
print(bessely(2.5, -10**50))

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

from sympy import*
from mpmath import*
mp.dps = 25; mp.pretty = True
print(bessely(2, 3.5, 1))
print(diff(lambda x: bessely(2, x), 3.5))
print(bessely(0.5, 3.5, 1))
print(diff(lambda x: bessely(0.5, x), 3.5))
print(diff(lambda x: bessely(2, x), 0.5, 10))
print(bessely(2, 0.5, 10))
print(bessely(2, 100.5, 100))
print(quad(lambda x: bessely(2,x), [1,3]))
print(bessely(2,3,-1) - bessely(2,1,-1))

Модифицированная функция Бесселя первого рода

mpmath.besseli(n, x, derivative=0)

besseli(n, x, derivative=0) модифицированная функция Бесселя первого рода

$I_{n}(x)= \mathit{i}^{-n}J_{n}(ix)$

$\frac{d^{m}}{dx^{m}}I_{n}(x)$

Модифицированная функция Бесселя $inline$ для вещественных порядков $inline$ :

from mpmath import*
i0 = lambda x: besseli(0,x)
i1 = lambda x: besseli(1,x)
i2 = lambda x: besseli(2,x)
i3 = lambda x: besseli(3,x)
plot([i0,i1,i2,i3],[0,5],[0,5])

Модифицированная функция Бесселя $I_{n}(x)$ в комплексной плоскости

from mpmath import*
cplot(lambda z: besseli(1,z), [-8,8], [-8,8], points=50000)

Примеры:
Некоторые значения $I_{n}(x)$

from mpmath import*
mp.dps = 25; mp.pretty = True
print(besseli(0,0))
print(besseli(1,0))
print(besseli(0,1))
print(besseli(3.5, 2+3j))

Аргументы могут быть большими:

from mpmath import*
mp.dps = 25; mp.pretty = True
print(besseli(2, 1000))
print(besseli(2, 10**10))
print(besseli(2, 6000+10000j))

Для целых чисел n выполняется следующее интегральное представление:

from mpmath import*
mp.dps = 15; mp.pretty = True
n = 3
x = 2.3
print(quad(lambda t: exp(x*cos(t))*cos(n*t), [0,pi])/pi)
print(besseli(n,x))

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

from mpmath import*
mp.dps = 25; mp.pretty = True
print(besseli(2, 7.5, 1))
print(diff(lambda x: besseli(2,x), 7.5))
print(besseli(2, 7.5, 10))
print(diff(lambda x: besseli(2,x), 7.5, 10))
print(besseli(2,7.5,-1) - besseli(2,3.5,-1))
print(quad(lambda x: besseli(2,x), [3.5, 7.5]))

Модифицированные функции Бесселя второго рода,

mpmath.besselk(n, x)

besselk(n, x) модифицированные функции Бесселя второго рода

$K_{n}(x)=\frac{\pi }{4}\frac{I_{-n}(x)-I_{n}(x)}{sin(\pi\cdot n)}$

Для целого числа $inline$ эту формулу следует понимать как предел.
Модифицированная функция Бесселя 2-го рода $K_{n}(x)$ для вещественных $inline$ :

from mpmath import*
k0 = lambda x: besselk(0,x)
k1 = lambda x: besselk(1,x)
k2 = lambda x: besselk(2,x)
k3 = lambda x: besselk(3,x)
plot([k0,k1,k2,k3],[0,8],[0,5])

Модифицированная функция Бесселя 2-го рода $K_{n}(x))$ в комплексной плоскости

from mpmath import*
cplot(lambda z: besselk(1,z), [-8,8], [-8,8], points=50000)

Примеры:
Сложные и комплексные аргументы:

from mpmath import *
mp.dps = 25; mp.pretty = True
print(besselk(0,1))
print(besselk(0, -1))
print(besselk(3.5, 2+3j))
print(besselk(2+3j, 0.5))

Аргументы – большие числа

from mpmath import *
mp.dps = 25; mp.pretty = True
print(besselk(0, 100))
print(besselk(1, 10**6))
print(besselk(1, 10**6*j))
print(besselk(4.5, fmul(10**50, j, exact=True)))

Особенности поведения функции в точке $inline$ :

from mpmath import *
print(besselk(0,0))
print(besselk(1,0))
for n in range(-4, 5):
    print(besselk(n, '1e-1000'))

Нули функции Бесселя

besseljzero()
mpmath.besseljzero(v, m, derivative=0)

Для реального порядка $\mathit{\nu}\geq 0$ и положительного целого числа $inline$ возвращает $j_{\nu,m}$ , m-й положительный нуль функции Бесселя первого рода $J_{\nu}(z)$ (см. besselj ()). Альтернативно, с $inline$ , дает первый неотрицательный простой ноль ${j}'_{\nu,m}$ из ${J}'_{\nu}(z)$ . Обозначения по соглашению об индексации с использованием Abramowitz & Stegun и DLMF. Обратите внимание на особый случай ${j}'_{0,1}= 0$ , в то время как все остальные нули положительны.

В действительности подсчитываются только простые нули (все нули функций Бесселя простые, за исключением когда $inline$ ), и $j_{\nu,m}$ становится монотонной функцией от $\nu$ и $inline$ . Нули чередуются согласно неравенствам:

${j}'_{\nu,k}< j_{\nu,k}<{j}'_{\nu,k+1}$

$j_{\nu,1}< j_{\nu+1,2}<j_{\nu,2}< j_{\nu+1,2} < j_{\nu,3}\cdots$

Примеры:
Начальные нули функций Бесселя $J_{0}(z)$ , $J_{1}(z)$ , $J_{2}(z)$

from mpmath import *
mp.dps = 25; mp.pretty = True
print(besseljzero(0,1))
print(besseljzero(0,2))
print(besseljzero(0,3))
print(besseljzero(1,1))
print(besseljzero(1,2))
print(besseljzero(1,3))
print(besseljzero(2,1))
print(besseljzero(2,2))
print(besseljzero(2,3))

Начальные нули производных от функций Бесселя ${J}'_{0}(z)$ , ${J}'_{1}(z)$ , ${J}'_{2}(z)$

from mpmath import *
mp.dps = 25; mp.pretty = True
print(besseljzero(0,1,1))
print(besseljzero(0,2,1))
print(besseljzero(0,3,1))
print(besseljzero(1,1,1))
print(besseljzero(1,2,1))
print(besseljzero(1,3,1))
print(besseljzero(2,1,1))
print(besseljzero(2,2,1))
print(besseljzero(2,3,1))

Нули с большим индексом:

from mpmath import *
mp.dps = 25; mp.pretty = True
print(besseljzero(0,100))
print(besseljzero(0,1000))
print(besseljzero(0,10000))
print(besseljzero(5,100))
print(besseljzero(5,1000))
print(besseljzero(5,10000))
print(besseljzero(0,100,1))
print(besseljzero(0,1000,1))
print(besseljzero(0,10000,1))

Нули функций с большим порядком:

from mpmath import *
mp.dps = 25; mp.pretty = True
print(besseljzero(50,1))
print(besseljzero(50,2))
print(besseljzero(50,100))
print(besseljzero(50,1,1))
print(besseljzero(50,2,1))
print(besseljzero(50,100,1))

Нули функций с дробным порядком:

from mpmath import *
mp.dps = 25; mp.pretty = True
print(besseljzero(0.5,1))
print(besseljzero(1.5,1))
print(besseljzero(2.25,4))

И $J_{\nu}(z)$ . и ${J}'_{\nu}(z)$ можно выразить как бесконечные произведения по их нулям:

from mpmath import *
mp.dps = 6; mp.pretty = True
v,z = 2, mpf(1)
nprint((z/2)**v/gamma(v+1) * \
                    nprod(lambda k: 1-(z/besseljzero(v,k))**2, [1,inf]))
print(besselj(v,z))
nprint((z/2)**(v-1)/2/gamma(v) * \
                        nprod(lambda k: 1-(z/besseljzero(v,k,1))**2, [1,inf]))
print(besselj(v,z,1))

besselyzero()
mpmath.besselyzero(v, m, derivative=0)

Для реального порядка $\mathit{\nu}\geq 0$ и положительного целого числа $inline$ возвращает $y_{\nu,m}$ , $inline$ , m-й положительный нуль функции Бесселя второго рода $Y_{\nu}(z)$ (см. Bessely ()). Альтернативно, с $inline$ , дает первый положительный нуль ${y}'_{\nu,m}$ из ${Y}'_{\nu}(z)$ . Нули чередуются согласно неравенствам:

$y_{\nu,k}< {y}'_{\nu,k}<y_{\nu,k+1}$

$y_{\nu,1}< y_{\nu+1,2}<y_{\nu,2}< y_{\nu+1,2} < y_{\nu,3}\cdots$

Примеры:
Начальные нули функций Бесселя $Y_{0}(z)$ , $Y_{1}(z)$ , $Y_{2}(z)$

from mpmath import *
mp.dps = 25; mp.pretty = True
print(besselyzero(0,1))
print(besselyzero(0,2))
print(besselyzero(0,3))
print(besselyzero(1,1))
print(besselyzero(1,2))
print(besselyzero(1,3))
print(besselyzero(2,1))
print(besselyzero(2,2))
print(besselyzero(2,3))

Начальные нули производных от функций Бесселя ${Y}'_{0}(z)$ , ${Y}'_{1}(z)$ , ${Y}'_{2}(z)$

from mpmath import *
mp.dps = 25; mp.pretty = True
print(besselyzero(0,1,1))
print(besselyzero(0,2,1))
print(besselyzero(0,3,1))
print(besselyzero(1,1,1))
print(besselyzero(1,2,1))
print(besselyzero(1,3,1))
print(besselyzero(2,1,1))
print(besselyzero(2,2,1))
print(besselyzero(2,3,1))

Нули с большим индексом:

from mpmath import *
mp.dps = 25; mp.pretty = True
print(besselyzero(0,100))
print(besselyzero(0,1000))
print(besselyzero(0,10000))
print(besselyzero(5,100))
print(besselyzero(5,1000))
print(besselyzero(5,10000))
print(besselyzero(0,100,1))
print(besselyzero(0,1000,1))
print(besselyzero(0,10000,1))

Нули функций с большим порядком:

from mpmath import *
mp.dps = 25; mp.pretty = True
print(besselyzero(50,1))
print(besselyzero(50,2))
print(besselyzero(50,100))
print(besselyzero(50,1,1))
print(besselyzero(50,2,1))
print(besselyzero(50,100,1))

Нули функций с дробным порядком:

from mpmath import *
mp.dps = 25; mp.pretty = True
print(besselyzero(0.5,1))
print(besselyzero(1.5,1))
print(besselyzero(2.25,4))

Приложения функций Бесселя
Важность функций Бесселя обусловлена не только частым появлением уравнения Бесселя в приложениях, но также и тем, что решения многих других линейных дифференциальных уравнений второго порядка могут быть выражены через функции Бесселя. Чтобы увидеть как они появляются, мы начнем с уравнения Бесселя порядка $inline$ в форме:

$z^{2}\frac{d^{2}w}{dz^{2}}+z\frac{dw}{dz}+(z^{2}-p^{2})w= 0, (1)$

и подставим сюда

$w=x^{-\alpha},z=kx^{\beta}, (2)$

Тогда, используя (2) и введя константы $inline$ из уравнения (1), получим:

$x^{2}{y}''+Ax{y}'+(B+Cx^{q})y= 0, (3)$

$A=1-2\alpha,B=\alpha ^{2}-\beta^{2}p^{2},C=\beta ^{2} k^{2},q=2\beta , (4)$

Из уравнения (4) получим:

$\left\{\begin{matrix} \alpha=\frac{1-A}{2}, \\ \beta =\frac{q}{2},\\ k=\frac{2\sqrt{C}}{q}\\ p=\frac{\sqrt{(1-A^{2}-4B}}{q}\\ \end{matrix}\right. (5)$

Если $inline$ , $q\neq 0$ , $(1-A)^{2}\geqslant 4B$ , то общее решение (для $inline$ ) уравнения (3) имеет вид:

$y(x)=x^{\alpha}\left [c_{1}J_{p}(kx^{\beta})+c_{2}J_{-p}(kx^{\beta })\right] (6)$

где: $\alpha$ , $\beta$ , $inline$ определяются из системы (5). Если $inline$ — целое число, то $J_{p}$ нужно заменить на $Y_{p}$ .

Продольный изгиб вертикальной колонны
Мы теперь рассмотрим задачу, важную для практических приложений. В этой задаче требуется определить, когда однородная вертикальная колонна согнется под ее собственным весом. Мы полагаем $inline$ в свободном верхнем конце колонны и $inline$ в её основании; мы предполагаем, что основание жестко вставлено (т. е. закреплено неподвижно) в основу (в землю), возможно в бетон.

Обозначим угловое отклонение колонны в точке $inline$ через $\theta(x)$ . Из теории эластичности при данных условиях следует, что:

$EI\frac{d^{2}\theta}{dx^{2}}+g\rho x \theta = 0, (7)$

где $inline$ — модуль Юнга материала колонны, $inline$ — момент инерции ее поперечного сечения, $\rho$ — линейная плотность колонны и $inline$ — гравитационное ускорение. Граничные условия имеют вид:

${\theta}'(0)=0,\theta(L)=0, (8)$

Будем решать задачу, используя (7) и (8) при:

$\lambda=\gamma^{2}=\frac{g\rho}{EI} (9)$

Перепишем (7) с учётом (9) при условии (8):

$\frac{d^{2}\theta}{dx^{2}}+\gamma^{2}x\theta=0;\frac{d\theta }{dx}=0,\theta(L)=0 .(10)$

Колонна может деформироваться, только если есть нетривиальное решение задачи (10); иначе колонна останется в не отклоненном от вертикали положении (т. е. физически не сможет отклониться от вертикали).
Дифференциальное уравнение (10) представляет собой уравнение Эйри. Уравнение (10) имеет форму уравнения (3) при $inline$ , $C=\gamma^{2}$ , $inline$ . Из системы уравнений (5) получаем $\alpha=\frac{1}{2}$ , $\beta=\frac{3}{2}$ , $k=\frac{2}{3}\gamma$ , $p=\frac{1}{3}$ .
Поэтому общее решение имеет вид:

$\theta(x)=x^{1/2}\left [c_{1}J_{1/3}(\frac{2}{3}\gamma x^{3/2})+c_{2}J_{-1/3}(\frac{2}{3}\gamma x^{3/2})\right ]. (11)$

Чтобы применить начальные условия, мы подставляем $p=\pm\frac{1}{3}$ в

$J_{p}=\sum_{m=0}^{\infty}\frac{(-1)^{m}}{m!\Gamma(p+m+1)}\left ( \frac{x}{2} \right )^{2m+p} (12)$

После преобразования (12), с учётом решения (11), получим:

$\theta(x)= \frac{c_{1}\gamma^{1/3}}{3^{1/3}\Gamma(4/3)}\left (x-\frac{\gamma^{2}x^{4} }{12}+\frac{\gamma ^{4}x^{7}}{504}-\cdot \cdot \cdot \right )+\\ +\frac{c_{2}3^{1/3}}{\gamma ^{1/3}\Gamma(\frac{2}{3})}\left (1-\frac{\gamma^{2}x^{3} }{6}+\frac{\gamma^{4}x^{6} }{180}-\cdot \cdot \cdot \right ) .(13)$

При условии в начальной точке ${\theta}'(0)=0$ , получим $c_{1}=0$ , тогда (11) примет вид:

$\theta(x)=c_{2}x^{1/2}J_{-1/3}(\frac{2}{3}\gamma x^{3/2}), (14)$

При условии в конечной точке $\theta(L)=0$ , из (14) получим:

$J_{-1/3}(\frac{2}{3}\gamma L^{3/2})=0 (15)$

Следует отметить, что преобразований (13), (14) можно было не делать, если построить графики функций $J_{1/3}(x), J_{-1/3}(x)$ , воспользовавшись рассмотренными возможностями модуля mpmath:

from mpmath import*
mp.dps = 6; mp.pretty = True
f=lambda x: besselj(-1/3,x)
f1=lambda x: besselj(1/3,x)
plot([f,f1], [0, 15])

Из графика следует, что при x=0 функция $J_{1/3}(0)=0$ и с учётом решения (11), мы сразу получаем необходимое уравнение (15), остаётся только найти z, как будет показано далее.

Таким образом, колонна деформируется, только если $z=\frac{2}{3}\gamma L^{3/2}$ — корень уравнения $J_{-1/3}(z)=0$ . Построим функцию $J_{-1/3}(z)$ на отдельном графике:

from mpmath import*
mp.dps = 6; mp.pretty = True
f=lambda x: besselj(-1/3,x)
plot(f, [0, 15])

На графике видно, что первый корень немного меньше 2. Найти корень $inline$ из уравнения $J_{-1/3}(z)=0$ можно, воспользовавшись функцией findroot(f, z0), приняв, согласно графика, точку поиска $inline$ и шесть знаков после запятой mp.dps = 6:

from mpmath import*
mp.dps = 6; mp.pretty = True
f=lambda x: besselj(-1/3,x)
print("z0=%s"%findroot(f, 1)

Получим:
$inline$
Рассчитаем критическую длину, например флагштока, по формуле (15):

Высота флагштока для разных параметров в сечении

from numpy import*
def LRr(R,r):
    E=2.9*10**11#н/м^2
    rou=7900#кг/м^3
    g=9.8#м/с^2
    I=pi*((R-r)**4)/4#м^4
    F=pi*(R-r)**2#м^2
    return 1.086*(E*I/(rou*g*F))**1/3
R=5*10**-3
r=0
L= LRr(R,r)
print(round(L,2),"м")
R=7.5*10**-3
r=2*10**-3
Lr= LRr(R,r)
print(round(Lr,2),"м")

Получим:
8.47 м
10.25 м

Полый флагшток может быть выше, чем сплошной.

Распространение волн в тонкой мембране.

Тонкая мембрана при попадании на неё звуковых волн не только колеблется с частотой волн. Форму колебаний мембраны можно получить в функциях Бесселя по следующему листингу, пользуясь формулами besselj() и besseljzero():

Форма колебаний мембраны

from mpmath import* 
from numpy import*
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
def Membrana(r):
    mp.dps=25
    return cos(0.5) * cos( theta) *float(besselj(1,r*besseljzero(1,1) ,0))
theta =linspace(0,2*pi,50)
radius = linspace(0,1,50)
x = array([r * cos(theta) for r in radius])
y = array([r * sin(theta) for r in radius])
z = array([Membrana(r) for r in radius])
fig = plt.figure("Мембрана")
ax = Axes3D(fig)
ax.plot_surface(x, y, z, rstride=1, cstride=1, cmap=cm.jet)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

Альтернатива модулю mpmath в специальных функциях Бесселя библиотеки SciPy

Не углубляясь в детальное рассмотрение функций Бесселя из библиотеки SciPy[4], приведу только два листинга для построения графиков функций первого и второго рода jv(v, x), yv(v, x):

jv(v, x)

import numpy as np
import pylab as py
import scipy.special as sp
x = np.linspace(0, 15, 500000)
for v in range(0, 6):
    py.plot(x, sp.jv(v, x))
py.xlim((0, 15))
py.ylim((-0.5, 1.1))
py.legend(('$J_{0}(x)$', '$ J_{1}(x)$', '$J_{2}(x)$',
           '$J_{3}(x)$', '$ J_{4}(x)$','$ J_{5}(x)$'),
           loc = 0)
py.xlabel('$x$')
py.ylabel('${J}_n(x)$')                             
py.grid(True)                                     
py.show()

yv(v, x)

import numpy as np
import pylab as py
import scipy.special as sp
x = np.linspace(0, 15, 500000)
for v in range(0, 6):
    py.plot(x, sp.yv(v, x))
py.xlim((0, 15))
py.ylim((-0.5, 1.1))
py.legend(('$Y_{0}(x)$', '$ Y_{1}(x)$', '$Y_{2}(x)$',
           '$Y_{3}(x)$', '$ Y_{4}(x)$','$ Y_{5}(x)$'),
           loc = 0)
py.xlabel('$x$')
py.ylabel('$Y_{n}(x)$')                             
py.grid(True)                                     
py.show()

Выводы:

В статье изложены основы работы с функциями Бесселя при помощи библиотек mpmath, sympy и scipy, приведены примеры применения функций для решения дифференциальных уравнений. Статья может быть полезна при изучении уравнений математической физики.

Ссылки:

1.Функции Бесселя
2.Использование обратного преобразования Лапласа для анализа динамических звеньев систем управления
3.Bessel function related functions
4. Special functions (scipy.special).

Функции Бесселя в программе символьной математики SymPy

Распространение волн в тонкой мембране.

Альтернатива модулю mpmath в специальных функциях Бесселя библиотеки SciPy

Выводы:

Ссылки:

Метод Каллана: я повторял предложения три года

Хайп вокруг ИИ: реальное состояние дел и будущее

Через два месяца первого класса мы забрали сына на семейное обучение

ИИ‑агент внутри КОМПАС-3D: пишет код, строит деталь и проверяет результат

Публикации