Введение:
Большое число самых разнообразных задач, относящихся практически ко всем важнейшим разделам математической физики и призванных ответить на актуальные технические вопросы, связано с применением функций Бесселя.
Функции Бесселя широко используются при решении задач акустики, радиофизики, гидродинамики, задач атомной и ядерной физики. Многочисленные приложения функций Бесселя к теории теплопроводности и теории упругости (задачи о колебаниях пластинок, задачи теории оболочек, задачи определения концентрации напряжения вблизи трещин).
Такая популярность функций Бесселя объясняется тем, что решение уравнений математической физики, содержащих оператор Лапласа в цилиндрических координатах, классическим методом разделения переменных приводит к обыкновенному дифференциальному уравнению, служащему для определения этих функций[1].
Функции Бесселя названы по имени немецкого астронома Фридриха Бесселя, который в работе 1824 года, изучая движение планет вокруг солнца, вывел рекуррентные соотношения для функций Бесселя , получил для целых интегральное представление функции , доказал наличие бесчисленного множества нулей функции и составил первые таблицы для функций и .
Однако, впервые одна из функций Бесселя была рассмотрена еще в 1732 году Даниилом Бернулли в работе, посвященной колебанию тяжелых цепей. Д. Бернулли нашел выражение функции в виде степенного ряда и заметил (без доказательства), что уравнение имеет бесчисленное множество действительных корней.
Следующей работой, в которой встречаются функции Бесселя, была работа Леонардо Эйлера 1738 года, посвященная изучению колебаний круглой мембраны. В этой работе Л. Эйлер нашел для целых выражение функции Бесселя в виде ряда по степеням , а в последующих работах распространил это выражение на случай произвольных значений индекса . Кроме того, Л. Эйлер доказал, что для , равного целому числу с половиной, функции выражаются через элементарные функции.
Он заметил (без доказательства), что при действительных функции имеют бесчисленное множество действительных нулей и дал интегральное представление для . Некоторые исследователи считают, что основные результаты, связанные с функциями Бесселя и их приложениями в математической физике, связаны с именем Л. Эйлера.
Изучить свойство функций Бесселя и одновременно освоить методы решения уравнений, сводящихся к функциям Бесселя, позволяет свободно распространяемая программа символьной математики SymPy — библиотеки Python.
В программе символьной математики SymPy графики функций Бесселя первого рода целых порядков можно построить, пользуясь соотношением для суммы ряда:
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()
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()
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()
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]
from sympy import*
from mpmath import*
cplot(lambda z: besselj(1,z), [-8,8], [-8,8], points=50000)
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))
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)])
from mpmath import*
print(quadosc(j0, [0, inf], period=2*pi))
print(quadosc(j1, [0, inf], period=2*pi))
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]))
from mpmath import*
mp.dps = 15
print(besselj(1, 3.5, 0.75))
print(differint(lambda x: besselj(1, x), 3.5, 0.75))
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])
from sympy import*
from mpmath import*
cplot(lambda z: bessely(1,z), [-8,8], [-8,8], points=50000)
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)
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])
from mpmath import*
cplot(lambda z: besseli(1,z), [-8,8], [-8,8], points=50000)
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))
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)
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])
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)))
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)
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))
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))
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)
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))
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))
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])
from mpmath import*
mp.dps = 6; mp.pretty = True
f=lambda x: besselj(-1/3,x)
plot(f, [0, 15])
from mpmath import*
mp.dps = 6; mp.pretty = True
f=lambda x: besselj(-1/3,x)
print("z0=%s"%findroot(f, 1)
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),"м")
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()
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()
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()
К сожалению, не доступен сервер mySQL