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

Question

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)

I think the sum must be multiplied with (x/2)**n. You add that value to the sum in your code — Stefan, Dec 18 '22 at 06:17
You amplify what `stefan` wrote, initialize `series_sum` to 0, and then return `series_sum * (x / 2)**n`, — Frank Yellin, Dec 18 '22 at 06:27
Thank you. I rewrite it as follows: ```# Initialize the power series expansion with the first term series_sum = 0 # 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 * (x / 2)**n``` But I just got 51.18424024866656. I thought it was supposedly around 0.364831230613667. — Becker Hija, Dec 18 '22 at 07:13
For inspiration of a more efficient way to compute such power series see https://stackoverflow.com/a/39409897/3088138, or in python https://stackoverflow.com/a/22791396/3088138 // Your formula for the current term has some errors. Why are there two powers of 2? The second one has the wrong exponent, parentheses are missing. — Lutz Lehmann, Dec 18 '22 at 09:39

score 3 · Accepted Answer · answered Dec 19 '22 at 18:49

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

Thank you mr Yves, this result is great. To be honest I am still beginner and I dont find literature for this. This is important since gas sensor using circular membrane nanofiber and circular electrode will use Bessel Function for Mass Sensitivity calculation. I am sorry if my question is kinda easy, but I dont find any book or article how to code Bessel function without any library from Python or Sympy — Becker Hija, Feb 12 '23 at 13:31

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

1 Answers1