與辛普森的規則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（θ）。當您評估積分時，只有參數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))