Python numerical integration with Simson's rule
Question: I have started to work through this book ([Computational Physics Exercise
5.4](http://www-personal.umich.edu/~mejn/cp/chapters/int.pdf)) and its
exercises and I got stuck with the following question:

[![enter image description
here](http://i.stack.imgur.com/A2cLq.png)](http://i.stack.imgur.com/A2cLq.png)

Write a Python function J(m,x) that calculates the value of Jm(x) using
Simpson’s rule with N = 1000 points. Use your function in a program to make a
plot, on a single graph, of the Bessel functions J0, J1, and J2 as a function
of x from x = 0 to x = 20.

I have created the following code to evaluate the first part of the question
but not sure if even this is correct:

    
    
    def f(x, t):
        return 1 / pi * (math.cos(x - t * math.sin(x)))
    
    def float_range(a, b, c):
        while a < b:
            yield a
            a += c  
        N = 1000
        a = 0.0
        b = 20.0
        h = (b - a) / N
        c = 0.0
        d = pi
        h2 = (d - c) / N
        s = 0.5 * f(a, 1) + 0.5 * f(b, 1)
        s / 3
        S1 = 0
        S2 = 0
        for k in range(1, N):
            for j in range(0, N):
                if k%2 == 0:
                    S1 += 2 / 3 * f(a + k * h, c + k * h2)
                else:
                    S2 += 4 / 3 * f(a + k * h, c + k * h2)
        s += S1 + S2
        print(h * s)
    

Could anyone please help me to solve this question, I have never used Bessel
functions before?


Answer: Your code is a little messy. You've defined a generator function `float_range`
but you never use it, and the rest of your code is tacked onto that function,
but should be separated from it, and not indented at the same level.

You've also messed up the definition of the core function that needs to be
integrated to get the Bessel function. You left out the constant `m`, and you
mixed up the `x` and the `t` (theta). When you're evaluating the integral only
the parameter `t` should vary - `m` and `x` are fixed.

Anyway, here's some working code:

    
    
    from math import sin, cos, pi
    
    # The core function of the Bessel integral
    def f(m, x, t):
        return cos(m * t - x * sin(t))
    
    #The number of steps used in the Simpson's integral approximation
    N = 1000
    
    def J(m, x):
        ''' Approximate Bessel function Jm(x) for integer m '''
    
        # lower & upper limits of the integral
        a = 0.0
        b = pi
    
        # step size
        h = (b - a) / N
    
        # Sum the values for Simpson's integration
        s = f(m, x, a) + f(m, x, b)
        for i in range(1, N):
            t = a + i * h
            if i % 2 == 1:
                s += 4.0 * f(m, x, t)
            else:
                s += 2.0 * f(m, x, t)
    
        # multiply by h/3 to get the integral
        # and divide by pi to get the Bessel function.
        return s * h / (3.0 * pi) 
    
    for x in range(21):
        print(x, J(0, x), J(1, x), J(2, x))
    

**output**

    
    
    0 1.0 3.59712259979e-17 5.16623780792e-17
    1 0.765197686558 0.440050585745 0.114903484932
    2 0.223890779141 0.576724807757 0.352834028616
    3 -0.260051954902 0.339058958526 0.486091260586
    4 -0.397149809864 -0.0660433280235 0.364128145852
    5 -0.177596771314 -0.327579137591 0.0465651162778
    6 0.150645257251 -0.276683858128 -0.24287320996
    7 0.30007927052 -0.00468282348235 -0.301417220086
    8 0.171650807138 0.234636346854 -0.112991720424
    9 -0.0903336111829 0.245311786573 0.144847341533
    10 -0.245935764451 0.0434727461689 0.254630313685
    11 -0.171190300407 -0.176785298957 0.139047518779
    12 0.0476893107968 -0.223447104491 -0.0849304948786
    13 0.206926102377 -0.0703180521218 -0.217744264242
    14 0.17107347611 0.133375154699 -0.152019882582
    15 -0.0142244728268 0.205104038614 0.0415716779753
    16 -0.174899073984 0.0903971756613 0.186198720941
    17 -0.169854252151 -0.0976684927578 0.158363841239
    18 -0.013355805722 -0.187994885488 -0.0075325148878
    19 0.14662943966 -0.105701431142 -0.157755906096
    20 0.167024664341 0.0668331241758 -0.160341351923
    

* * *

The accuracy of this approximation is surprisingly good. Here are some values
generated using the [mpmath](http://mpmath.org/) module, which supplies
various Bessel functions.

    
    
    0 1.0 0.0 0.0
    1 0.76519768655796655145 0.44005058574493351596 0.11490348493190048047
    2 0.22389077914123566805 0.5767248077568733872 0.35283402861563771915
    3 -0.26005195490193343762 0.33905895852593645893 0.48609126058589107691
    4 -0.39714980986384737229 -0.066043328023549136143 0.36412814585207280421
    5 -0.17759677131433830435 -0.32757913759146522204 0.046565116277752215532
    6 0.15064525725099693166 -0.27668385812756560817 -0.24287320996018546772
    7 0.30007927051955559665 -0.0046828234823458326991 -0.30141722008594012028
    8 0.17165080713755390609 0.23463634685391462438 -0.11299172042407525
    9 -0.090333611182876134336 0.24531178657332527232 0.14484734153250397263
    10 -0.2459357644513483352 0.04347274616886143667 0.25463031368512062253
    11 -0.17119030040719608835 -0.17678529895672150114 0.13904751877870126996
    12 0.047689310796833536624 -0.22344710449062761237 -0.084930494878604805352
    13 0.206926102377067811 -0.070318052121778371157 -0.21774426424195679117
    14 0.17107347611045865906 0.13337515469879325311 -0.15201988258205962291
    15 -0.014224472826780773234 0.20510403861352276115 0.04157167797525047472
    16 -0.17489907398362918483 0.090397175661304186239 0.18619872094129220811
    17 -0.16985425215118354791 -0.097668492757780650236 0.15836384123850347142
    18 -0.013355805721984110885 -0.18799488548806959401 -0.0075325148878013995603
    19 0.14662943965965120426 -0.1057014311424092668 -0.15775590609569428497
    20 0.16702466434058315473 0.066833124175850045579 -0.16034135192299815017
    

  

I'll let you figure out the plotting part of your exercise. :)

* * *

Here's the mpmath code I used to generate those Bessel function values above,
which are accurate to 20 significant figures:

    
    
    from mpmath import mp
    
    # set precision
    mp.dps = 20
    
    for x in range(21):
        print(x, mp.besselj(0, x), mp.besselj(1, x), mp.besselj(2, x))
    

