Python 数值积分与辛普森规则

Question

我已经开始研究这本书 ( Computational Physics Exercise 5.4 ) 及其练习，但我遇到了以下问题:

编写一个 Python 函数 J(m,x)，使用 N = 1000 点的辛普森法则计算 Jm(x) 的值。在程序中使用您的函数在单个图形上绘制贝塞尔函数 J0、J1 和 J2 作为 x 从 x = 0 到 x = 20 的函数。

我创建了以下代码来评估问题的第一部分，但不确定这是否正确:

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)

谁能帮我解决这个问题，我以前从未使用过贝塞尔函数？

Answer 1

你的代码有点乱。您已经定义了一个生成器函数 float_range 但您从未使用过它，并且您的其余代码附加到该函数上，但应该与它分开，并且不在同一级别缩进。

您还弄乱了需要集成以获得贝塞尔函数的核心函数的定义。您遗漏了常量 m，并且混淆了 x 和 t (theta)。当您评估积分时，只有参数 t 应该变化 - m 和 x 是固定的。

无论如何，这里有一些工作代码:

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))

输出

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

---

这种近似的准确性出奇地好。以下是使用 mpmath 生成的一些值模块，提供各种贝塞尔函数。

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

我会让你弄清楚练习的绘图部分。 :)

---

这是我用来生成上述贝塞尔函数值的 mpmath 代码，它精确到 20 位有效数字:

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))