Giải thích phép cầu phương: Các phương pháp số để xấp xỉ tích phân

Có một nhóm tích phân xác định mà hoàn toàn không thể giải được bằng phép tính tay. Không phải vì bạn làm sai. Một số hàm, như e^(−x²) hoặc sin(x)/x, hoàn toàn không có nguyên hàm dạng đóng. Khi điều đó xảy ra, bạn cần cách tiếp cận khác: đánh giá hàm tại một vài điểm được chọn lọc cẩn thận và kết hợp các giá trị đó để ước tính diện tích bên dưới. Đó là phép cầu phương.

Bài viết này sẽ trình bày phép cầu phương là gì, khuôn khổ tổng quát hoạt động ra sao, và các phương pháp chính, từ quy tắc điểm giữa trực quan đến cách tiếp cận Gaussian chính xác đáng ngạc nhiên. Chúng tôi cũng sẽ đề cập đến các phương pháp ghép (composite), nguồn gốc sai số, chiến lược thích nghi, và cách chọn đúng phương pháp cho bài toán trước mắt bạn.

Phép cầu phương là gì?

Từ "quadrature" có từ lâu đời. Nó xuất phát từ bài toán hình học cổ điển về việc dựng một hình vuông có cùng diện tích với một hình đã cho, tức là một cách khác để nói "tìm diện tích". Ngày nay, thuật ngữ này có nghĩa là xấp xỉ số của một tích phân xác định.

Ý tưởng cốt lõi là một sự thay thế: thay vì tìm nguyên hàm chính xác, bạn lấy mẫu hàm f(x) tại một tập điểm gọi là các nút, nhân mỗi giá trị với một trọng số đi kèm, rồi cộng lại:

Các thuật ngữ quadrature và numerical integration về cơ bản có cùng nghĩa. Quadrature thường xuất hiện trong các ngữ cảnh toán học nhiều hơn; numerical integration là thuật ngữ thường ngày trong tính toán khoa học.

Cách hoạt động của tích phân số (quadrature)

Mọi quy tắc cầu phương đều hoạt động trên một khoảng [a, b]. Trong khoảng đó, nó chọn một tập nút xᵢ và gán cho mỗi nút một trọng số wᵢ. Xấp xỉ là tổng có trọng số như trên. Điểm khác nhau giữa các phương pháp là vị trí đặt các nút và cách gán trọng số. Vậy thôi. Quy tắc hình thang, quy tắc Simpson và cầu phương Gaussian đều là các biến thể của cùng một khuôn mẫu này.

Sai số xấp xỉ là khoảng cách giữa tổng có trọng số và tích phân thật. Mỗi phương pháp đều tạo ra một phần sai số, và tốc độ sai số đó giảm khi bạn thêm nhiều nút hơn (hoặc làm nhỏ các khoảng) là điều phân biệt một phương pháp tốt với một phương pháp xuất sắc. Hãy giữ khung tư duy này: chọn một phương pháp cầu phương thực chất là chọn cách "tiêu" các lần đánh giá hàm của bạn.

Các phương pháp cầu phương phổ biến

Các phương pháp dưới đây đi từ trực quan nhất đến tinh vi nhất.

So sánh trực quan các quy tắc tích phân Điểm giữa, Hình thang và Simpson. Ảnh do Tác giả cung cấp.

Quy tắc điểm giữa (Midpoint)

Lấy điểm giữa của [a, b], đánh giá f tại đó, rồi nhân với bề rộng (b − a). Giả định ở đây là giá trị hàm tại trung điểm đại diện cho trung bình trên toàn khoảng.

Nghe có vẻ quá đơn giản, và với một khoảng đơn lẻ thì đúng là vậy. Nhưng đây là bất ngờ đầu tiên: quy tắc điểm giữa là chính xác với các hàm hằng và tuyến tính, và sai số của nó tỉ lệ O(h³) trên mỗi khoảng, với h = b − a. Điều đó thực ra còn tốt hơn quy tắc hình thang trên các hàm trơn, điều thường khiến người ta ngạc nhiên.

Quy tắc hình thang (Trapezoidal)

Dùng cả hai đầu mút thay vì chỉ trung điểm. Nối chúng bằng một đoạn thẳng và tính diện tích hình thang thu được: (b − a) · [f(a) + f(b)] / 2. Chính xác với hàm tuyến tính, có sai số O(h³) trên mỗi khoảng, cùng bậc với quy tắc điểm giữa nhưng với hằng số lớn hơn.

Quy tắc hình thang có lẽ là phương pháp cầu phương được dạy rộng rãi nhất, vì thế danh tiếng của nó có phần vượt trội. Với các hàm "ngoan" thì nó ổn. Nhưng quy tắc Simpson ghép (composite Simpson) gần như luôn làm tốt hơn với cùng số lần đánh giá.

Quy tắc Simpson

Khớp một đa thức bậc hai qua ba điểm: hai đầu mút và trung điểm. Tích phân chính xác đa thức đó:

Đánh giá thêm tại trung điểm chỉ tốn một lần gọi hàm, nhưng phần thưởng là có thật: quy tắc Simpson chính xác với đa thức đến bậc ba, và sai số trên mỗi khoảng tỉ lệ O(h⁵), tốt hơn hai bậc so với quy tắc hình thang. Với các hàm trơn, đây thường là phương pháp đầu tiên đáng để bạn sử dụng.

Cầu phương Gaussian

Đây là nơi mọi thứ trở nên thực sự thú vị. Các quy tắc điểm giữa, hình thang và Simpson đều dùng các nút cách đều. Cầu phương Gaussian đặt ra câu hỏi khác: nếu ta chọn cả vị trí nút lẫn trọng số để tối đa hóa độ chính xác thì sao?

Với n nút, cầu phương Gaussian (cụ thể là biến thể Gauss-Legendre) có thể tích phân chính xác mọi đa thức đến bậc 2n − 1. Hai nút cho bạn độ chính xác với đa thức bậc ba; năm nút xử lý đa thức bậc 9. Các nút là nghiệm của các đa thức Legendre, và chúng cụm về phía các đầu mút của khoảng thay vì nằm ở các vị trí cách đều. Độ chính xác tăng thêm đến từ việc cho phép các nút di chuyển, không phải từ vị trí tình cờ của chúng. Việc chúng cụm gần các đầu mút là hệ quả của toán học, và ta sẽ thấy vì sao bên dưới.

Dưới đây là cách cả bốn phương pháp sosánh:

Các phương pháp cầu phương ghép (Composite)

Áp dụng một quy tắc trên toàn khoảng [a, b] thường là không đủ. Cách khắc phục tiêu chuẩn: chia khoảng thành n khoảng con bằng nhau có bề rộng h = (b − a)/n và áp dụng quy tắc cơ bản cho từng khoảng con.

Quy tắc điểm giữa và hình thang ghép

Với quy tắc hình thang ghép, đánh giá f tại n+1 điểm cách đều và cộng n hình thang. Tổng sai số giảm xuống O(h²). Giảm một nửa h thì sai số giảm bốn lần. Quy tắc điểm giữa ghép theo cùng logic và có bậc sai số toàn cục tương tự.

Quy tắc Simpson ghép

Phiên bản ghép yêu cầu số khoảng con chẵn. Đánh giá f tại n+1 điểm và đạt sai số toàn cục O(h⁴). Giảm một nửa h khiến sai số giảm mười sáu lần. Đó là lý do quy tắc Simpson ghép thường là mặc định khi bạn cần độ chính xác vừa phải trên một hàm trơn.

Đoạn mã Python dưới đây cho thấy xấp xỉ chặt hơn như thế nào khi bạn thêm các khoảng con cho ∫[0 đến π] sin(x) dx (giá trị chính xác: 2.0):

import numpy as np

def composite_simpson(func, lower, upper, num_intervals):
    """Approximate integral using composite Simpson's rule."""
    if num_intervals % 2 != 0:
        raise ValueError("num_intervals must be even")
    nodes = np.linspace(lower, upper, num_intervals + 1)
    values = func(nodes)
    h = (upper - lower) / num_intervals
    weights = np.ones(num_intervals + 1)
    weights[1:-1:2] = 4   # odd-indexed interior points
    weights[2:-2:2] = 2   # even-indexed interior points
    return (h / 3) * np.dot(weights, values)

exact = 2.0
for n in [4, 8, 16, 32]:
    approx = composite_simpson(np.sin, 0, np.pi, n)
    print(f"n={n:3d}: approx={approx:.10f}, error={abs(approx - exact):.2e}")

Chỉ với 32 khoảng, sai số giảm xuống dưới 10⁻⁶. Đó là hội tụ O(h⁴).

Cầu phương Gaussian

Ý tưởng đằng sau cầu phương Gaussian xứng đáng nhiều hơn một ô trong bảng.

Vì sao các nút cách đều không phải lúc nào cũng tối ưu

Khi các nút nằm tại các điểm cách đều, bạn đang chi tiêu ngân sách đánh giá theo cách có thể không khớp với độ khó của hàm tích phân. Một số vùng cần độ phân giải cao hơn; số khác lại lãng phí đánh giá trên những đoạn phẳng. Cầu phương Gaussian tránh điều này bằng cách chọn các nút được chứng minh là tối ưu cho bậc đa thức mà bạn có thể tích phân chính xác.

Với n nút, quy tắc Gauss-Legendre chọn n nghiệm của đa thức Legendre bậc n, Pₙ(x), làm các nút, rồi tính trọng số sao cho mọi đa thức có bậc ≤ 2n − 1 đều được tích phân không sai số. Trong thực tế, numpy.polynomial.legendre.leggauss(n) trả về cả hai. Các nút cụm về phía hai đầu của [−1, 1], khoảng chuẩn tham chiếu, nơi sai số nội suy đa thức có xu hướng lớn nhất và nơi mật độ lấy mẫu bổ sung mang lại hiệu quả cao nhất.

Khi nào nên dùng

Cầu phương Gaussian phát huy tốt nhất khi bạn có một hàm trơn, giải tích và cần độ chính xác cao với ít lần đánh giá. Năm điểm Gauss-Legendre có thể vượt trội so với quy tắc Simpson ghép với hàng chục khoảng trên các hàm "ngoan". Nhưng nó không phù hợp với các hàm có gấp khúc, gián đoạn, hoặc dữ liệu chỉ có tại các điểm cách đều cố định. Trong những trường hợp đó, việc đặt nút được chọn lọc cẩn thận lại phản tác dụng.

Sai số và độ chính xác trong cầu phương

Không có phương pháp số nào miễn sai số. Hai nguồn gốc riêng biệt ảnh hưởng đến mọi kết quả cầu phương, và hiểu cả hai là điều cốt yếu để biết khi nào nên tin kết quả của mình.

Các nguồn sai số gồm:

• Sai số rời rạc hóa (còn gọi là sai số cắt cụt) là khoảng cách giữa tích phân chính xác và giá trị quy tắc cầu phương trả về trong số học hoàn hảo. Đây là điều mà các công thức O(hᵏ) phản ánh. Nó giảm khi bạn thêm nhiều điểm hơn hoặc giảm bề rộng khoảng.
• Sai số dấu chấm động là vấn đề khác. Cộng nhiều lượng nhỏ trong số học độ chính xác hữu hạn sẽ tạo làm tròn không giảm theo h và thậm chí có thể tăng. Với các quy tắc ghép dùng hàng nghìn khoảng con, bạn có thể chạm mức sàn nơi nhiễu dấu chấm động chi phối và việc thêm khoảng không còn giúp nữa. Với hầu hết tác vụ tính toán khoa học, bạn sẽ đạt chuẩn dung sai mong muốn trước khi chạm mức sàn đó. Nhưng nó tồn tại.

Độ trơn của hàm, bề rộng khoảng và lựa chọn phương pháp đều ảnh hưởng đến sai số rời rạc hóa. Một hàm gồ ghề với đặc trưng sắc nét cần nhiều khoảng hơn hoặc một cách tiếp cận thích nghi.

Cầu phương thích nghi (Adaptive)

Lưới đều lãng phí các lần đánh giá khi hàm tích phân trơn ở hầu hết nơi nhưng phức tạp ở một vài chỗ.

Cách hoạt động

Áp dụng một quy tắc cầu phương cho một khoảng, rồi áp dụng một quy tắc bậc cao hơn (hoặc cùng quy tắc trên hai nửa) và so sánh. Chênh lệch ước tính sai số cục bộ. Nếu nó dưới ngưỡng dung sai, chấp nhận khoảng đó. Nếu không, chia đôi và đệ quy trên mỗi nửa. Lặp lại cho đến khi mọi khoảng con đáp ứng dung sai.

Lợi ích hiệu quả so với tinh chỉnh đều có thể rất lớn khi hàm tích phân có các đặc trưng cô lập: một đỉnh hẹp, một gấp khúc, dao động nhanh bao quanh bởi vùng trơn. Các phương pháp thích nghi tập trung các lần đánh giá đúng vào nơi hàm thực sự cần.

scipy.integrate.quad()

Trong Python, scipy.integrate.quad() triển khai cầu phương Gaussian thích nghi (cụ thể là thuật toán DQAGSE của QUADPACK). Truyền vào hàm và cận; nhận về ước lượng và một chặn sai số. Với hầu hết tích phân một chiều, đây là công cụ đúng đắn, không cần tinh chỉnh thủ công.

from scipy import integrate
import numpy as np

# Integrate a function with a peak near x=0.5
def peaked_func(x):
    return np.exp(-100 * (x - 0.5)**2)

result, error_bound = integrate.quad(peaked_func, 0, 1)
print(f"Estimate: {result:.10f}")
print(f"Error bound: {error_bound:.2e}")

quad() xử lý đỉnh hẹp này gọn gàng. Một quy tắc Simpson đều sẽ cần hàng trăm khoảng con để bắt kịp.

Cách chọn phương pháp cầu phương

Bắt đầu với scipy.integrate.quad() cho các tích phân một chiều trong Python trừ khi bạn có lý do khác. Nó thích nghi, chính xác và đã được kiểm chứng tốt. Mọi điều bên dưới hữu ích khi bạn tự triển khai, làm việc trong ngôn ngữ không có thư viện thích nghi tốt, hoặc tối ưu cho trường hợp cụ thể.

• Hàm trơn, giải tích với công thức biết trước: Cầu phương Gaussian cho độ chính xác cao nhất trên mỗi lần đánh giá. Nếu bạn cần độ chính xác vừa phải và đơn giản, quy tắc Simpson ghép là mặc định vững vàng.
• Hàm có gấp khúc hoặc gián đoạn: Các phương pháp cách đều (hình thang, Simpson) thường vượt trội hơn cầu phương Gaussian ở đây, vì vị trí nút cố định không bị xáo trộn bởi bất thường. Chia khoảng tích phân tại điểm gián đoạn đã biết và xử lý từng phần riêng.
• Dữ liệu tại các điểm cách đều cố định: Cầu phương Gaussian không khả dụng. Bạn không thể đánh giá hàm tại các nút tùy ý, vì vậy hãy dùng quy tắc hình thang ghép hoặc Simpson ghép.
• Tích phân nhiều chiều: Cầu phương vấp phải "lời nguyền chiều" rất nhanh. Với n điểm mỗi chiều và d chiều, bạn cần tổng cộng nᵈ lần đánh giá. Mười chiều với mười điểm mỗi chiều là 10¹⁰ lần đánh giá, điều bất khả thi. Monte Carlo mở rộng tốt hơn nhiều trên khoảng hơn năm chiều, nên đó là nơi cần chuyển hướng.

Cầu phương so với các phương pháp số liên quan

Cầu phương không tồn tại trong chân không. Sẽ hữu ích khi thấy nó nằm ở đâu so với các công cụ khác mà bạn có thể cân nhắc.

Cầu phương vs. tích phân giải tích

Tích phân giải tích tìm nguyên hàm chính xác; cầu phương xấp xỉ một con số. Nếu tồn tại nguyên hàm chính xác và bạn có thể tìm ra, hãy dùng nó. Cầu phương có chỗ đứng khi tích phân chính xác là bất khả hoặc không thực tiễn.

Cầu phương vs. tích phân Monte Carlo

Cả hai đều xấp xỉ tích phân xác định, nhưng theo cách khác nhau. Cầu phương là tất định: nó đánh giá f tại các nút cố định, được chọn trước, với tốc độ hội tụ gắn với độ trơn của hàm và số lượng nút. Monte Carlo rút mẫu ngẫu nhiên và lấy trung bình, với sai số tỉ lệ O(1/√N) bất kể số chiều. Đó là hội tụ chậm hơn ở số chiều thấp, nhưng độc lập với chiều, khiến nó là lựa chọn thực tiễn trên khoảng hơn năm chiều.

Cầu phương vs. nội suy

Một số quy tắc cầu phương, bao gồm quy tắc Simpson, hoạt động bằng cách khớp một đa thức qua các điểm đánh giá rồi tích phân chính xác nó. Nhưng mục tiêu khác với nội suy. Nội suy tái dựng hàm; cầu phương chỉ cần một con số. Đa thức là công cụ, không phải kết quả đầu ra.

Kết luận

Điều khiến tôi thấy thỏa mãn về cầu phương là một ý tưởng tổ chức duy nhất lại mang đến nhiều kết quả đến vậy. Mọi phương pháp ở đây, điểm giữa, hình thang, Simpson, Gaussian, thực chất chỉ là những câu trả lời khác nhau cho cùng một câu hỏi: với một ngân sách số lần đánh giá hàm, bạn đặt chúng ở đâu và cân chúng như thế nào?

Khớp các đa thức tốt hơn qua các điểm đánh giá giúp kéo sai số xuống nhanh. Dịch chuyển sang các nút được chọn tối ưu vắt thêm độ chính xác từ cùng ngân sách. Các phương pháp thích nghi sau đó kết nối mọi thứ, điều phối các lần đánh giá đến đúng nơi hàm thực sự cần.

Trong Python, scipy.integrate.quad() là mặc định đúng cho các bài toán 1D. Mọi thứ khác trong bài viết này dành cho khi bạn cần hiểu chuyện gì đang diễn ra bên dưới, hoặc khi bạn tự triển khai. Nếu bạn muốn tiếp tục xây dựng kỹ năng số với SciPy, hãy học khóa Introduction to Optimization in Python của chúng tôi.