How to write Bessel function using power series method in Python without Sympy?-CodePudding

I am studying Computational Physics with a lecturer who always ask me to write Python and Matlab code without using instant code (a library that gives me final answer without showing mathematical expression). So I try to write Bessel function for first kind using power series because I thought it was easy compare to other method (I am not sure). I dont know why the result is still very different? Far from answer that Sympy.special provided?

Here is my code for x = 5 and n = 3

import math

def bessel_function(n, x, num_terms):
  # Initialize the power series expansion with the first term
  series_sum = (x / 2) ** n

  # Calculate the remaining terms of the power series expansion
  for k in range(0, num_terms):
    term = ((-1) ** k) * ((x / 2) ** (2 * k)) / (math.factorial(k)**2)*(2**2*k)
    series_sum = series_sum   term

  return series_sum

# Test the function with n = 3, x = 5, and num_terms = 10
print(bessel_function(3, 5, 30))
print(bessel_function(3, 5, 15))

And here is the code using sympy library:

from mpmath import *
mp.dps = 15; mp.pretty = True
print(besselj(3, 5))

import sympy

def bessel_function(n, x):
  # Use the besselj function from sympy to calculate the Bessel function
  return sympy.besselj(n, x)

# Calculate the numerical value of the Bessel function using evalf
numerical_value = bessel_function(3, 5).evalf()
print(numerical_value)

CodePudding user response：

It is a big waste to compute the terms like you do, each from scratch with power and factorial. Much more efficient to compute a term from the previous.

For Jn,

Tk / Tk-1 = - (X/2)²/(k(k N))

with T0 = (X/2)^N/N!.

N= 3
X= 5

# First term
X*= 0.5
Term= pow(X, N) / math.factorial(N)
Sum= Term
print(Sum)

# Next terms
X*= -X
for k in range(1, 21):
    Term*= X / (k * (k   N))
    Sum = Term
    print(Sum)

The successive sums are

2.6041666666666665
-1.4648437499999996
1.0782877604166665
0.19525598596643523
0.39236129276336185
0.3615635885763421
0.365128137672062
0.3648098743599441
0.3648324782883616
0.36483117019065225
0.3648312330799652
0.36483123052763916
0.3648312306162616
0.3648312306135987
0.3648312306136686
0.364831230613667
0.36483123061366707
0.36483123061366707
0.36483123061366707
0.36483123061366707
0.36483123061366707

CodePudding user response：

Well i tried to use Matlab with something that similar to gamma function in Python. So I use integration method that I read from my notes.

J_n(x) = (1/π) * ∫_0^π cos(nt - xsin(t)) dt

So here is my matlab code:

n = 3;  % order of the Bessel function
x = 5;  % argument of the Bessel function

Jn = @(t) cos(n*t - x*sin(t));  % Integral form
bessel = (1/pi) * integral(Jn, 0, pi);  % evaluate the integral

Finally i get 0.3648. Is there any something that better than this? Perhaps I still got misunderstanding