Hướng dẫn midpoint integration in python - tích hợp điểm giữa trong python

Phương pháp điểm giữa tổng hợp

Ý tưởng

Thay vì xấp xỉ khu vực dưới một đường cong bằng hình thang, chúng ta có thể sử dụng hình chữ nhật đơn giản. Nghe có vẻ kém chính xác khi sử dụng các đường ngang và không lệch các đường theo hàm được tích hợp, nhưng một phương pháp tích hợp dựa trên hình chữ nhật (phương pháp trung điểm) trên thực tế là chính xác hơn một chút so với phương pháp dựa trên hình thang!

Trong phương thức điểm giữa, chúng tôi xây dựng một hình chữ nhật cho mọi khoảng thời gian phụ trong đó chiều cao bằng \ (f \) ở điểm giữa của khoảng thời gian phụ. Chúng ta hãy làm điều này cho bốn hình chữ nhật, sử dụng các giao dịch phụ giống như chúng ta đã có để tính toán bằng tay với phương pháp hình thang: \ ([0,0.2) \), \ ([0.2,0.6) \), \ ([0.6, 0,8) \) và \ ([0,8,1.0] \). Chúng tôi nhận được $$ \ BẮT ĐẦU {Align} \ int_0^1 f (t) dt & \ xấp xỉ h_1 f \ left (\ frac {0 + 0.2} 0.6} {2} \ right) \ nonumber \\ & + h_3 f \ left (\ frac {0.6 + 0.8} {2} \ right) + h_4 f \ left (\ frac {0.8 + 1.0} ), \ tag {19} \ end {align} $$ trong đó \ (h_1 \), \ (h_2 \), \ (h_3 \) và \ (h_4 \) với phương pháp hình thang và được xác định trong (11)-(14).

Hình 9: Tính toán khoảng tích phân của một hàm là tổng của các khu vực của các hình chữ nhật.

Với \ (f (t) = 3t^{2} e^{t^3} \), xấp xỉ trở thành \ (1.632 \). So với câu trả lời thực sự (\ (1.718 \)), đây là khoảng \ (5 \% \) quá nhỏ, nhưng nó tốt hơn những gì chúng ta có với phương pháp hình thang (\ (10 ​​\% \)) Sub-Intervals. Nhiều hình chữ nhật cho một xấp xỉ tốt hơn.

Công thức chung

Chúng ta hãy lấy một công thức cho phương thức điểm giữa dựa trên các hình chữ nhật \ (n \) có chiều rộng bằng nhau: $$ \ start . H f \ left (\ frac {x_0 + x_1}} {2} \ right) + h f \ left (\ frac {x_1 + x_2} {2} 1} + x_n} {2} \ right), \ tag {20} \\ & \ xấp xỉ h \ left (f \ left (\ frac {x_0 + x_1} A \ TAG {21} \ end {align} $$ Tổng số này có thể được viết nhỏ gọn hơn là $$ \ BEGIN {Công thức} \ int_a^b f (x) } f (x_i), \ tag {22} \ end {phương trình} $$ trong đó \ (x_i = \ left (a + \ frac {h} {2} \ right) + ih \).

Thực hiện

Chúng tôi làm theo lời khuyên và bài học kinh nghiệm từ việc thực hiện phương pháp hình thang và tạo hàm

    n        midpoint          trapezoidal
      2 0.8842000076332692 0.8770372606158094
      4 0.8827889485397279 0.8806186341245393
      8 0.8822686991994210 0.8817037913321336
     16 0.8821288703366458 0.8819862452657772
     32 0.8820933014203766 0.8820575578012112
     64 0.8820843709743319 0.8820754296107942
    128 0.8820821359746071 0.8820799002925637
    256 0.8820815770754198 0.8820810181335849
    512 0.8820814373412922 0.8820812976045025
   1024 0.8820814024071774 0.8820813674728968
   2048 0.8820813936736116 0.8820813849400392
   4096 0.8820813914902204 0.8820813893068272
   8192 0.8820813909443684 0.8820813903985197
  16384 0.8820813908079066 0.8820813906714446
  32768 0.8820813907737911 0.8820813907396778
  65536 0.8820813907652575 0.8820813907567422
 131072 0.8820813907631487 0.8820813907610036
 262144 0.8820813907625702 0.8820813907620528
 524288 0.8820813907624605 0.8820813907623183
1048576 0.8820813907624268 0.8820813907623890

def midpoint(f, a, b, n):
    h = float(b-a)/n
    result = 0
    for i in range(n):
        result += f((a + h/2.0) + i*h)
    result *= h
    return result

Chúng tôi có thể kiểm tra chức năng như chúng tôi đã giải thích cho phương pháp

>>> a = 1; b = 2; expected = 3
>>> a + b == expected
True

So sánh các phương pháp hình thang và trung điểm

Ví dụ tiếp theo cho thấy chúng ta có thể kết hợp các chức năng

>>> a = 1; b = 2; expected = 3
>>> a + b == expected
True

>>> a = 1; b = 2; expected = 3
>>> a + b == expected
True

from trapezoidal import trapezoidal
from midpoint import midpoint
from math import exp

g = lambda y: exp(-y**2)
a = 0
b = 2
print '    n        midpoint          trapezoidal'
for i in range(1, 21):
    n = 2**i
    m = midpoint(g, a, b, n)
    t = trapezoidal(g, a, b, n)
    print '%7d %.16f %.16f' % (n, m, t)

Lưu ý những nỗ lực được đặt vào định dạng đẹp - đầu ra trở thành

    n        midpoint          trapezoidal
      2 0.8842000076332692 0.8770372606158094
      4 0.8827889485397279 0.8806186341245393
      8 0.8822686991994210 0.8817037913321336
     16 0.8821288703366458 0.8819862452657772
     32 0.8820933014203766 0.8820575578012112
     64 0.8820843709743319 0.8820754296107942
    128 0.8820821359746071 0.8820799002925637
    256 0.8820815770754198 0.8820810181335849
    512 0.8820814373412922 0.8820812976045025
   1024 0.8820814024071774 0.8820813674728968
   2048 0.8820813936736116 0.8820813849400392
   4096 0.8820813914902204 0.8820813893068272
   8192 0.8820813909443684 0.8820813903985197
  16384 0.8820813908079066 0.8820813906714446
  32768 0.8820813907737911 0.8820813907396778
  65536 0.8820813907652575 0.8820813907567422
 131072 0.8820813907631487 0.8820813907610036
 262144 0.8820813907625702 0.8820813907620528
 524288 0.8820813907624605 0.8820813907623183
1048576 0.8820813907624268 0.8820813907623890

Remark.

Các phương pháp hình thang và trung điểm chỉ là hai ví dụ trong một khu rừng của các quy tắc tích hợp số. Các phương pháp nổi tiếng khác là quy tắc của Simpson và Gauss quaratre. Tất cả đều hoạt động theo cùng một cách: $$ \ int_a^b f (x) dx \ xấp xỉ \ sum_ {i = 0}^{n-1} w_if (x_i) Một tổng các đánh giá chức năng, trong đó mỗi đánh giá \ (f (x_i) \) được cho trọng lượng \ (w_i \). Các phương thức khác nhau khác nhau trong cách chúng xây dựng các điểm đánh giá \ (x_i \) và trọng số \ (w_i \). Chúng tôi đã sử dụng các điểm cách đều nhau \ (x_i \), nhưng độ chính xác cao hơn có thể đạt được bằng cách tối ưu hóa vị trí của \ (x_i \).

Kiểm tra

Các vấn đề với quy trình kiểm tra ngắn gọn

Kiểm tra các chương trình để tích hợp số cho đến nay đã sử dụng hai chiến lược. Nếu chúng tôi có câu trả lời chính xác, chúng tôi tính toán lỗi và thấy rằng việc tăng \ (n \) sẽ giảm lỗi. Khi câu trả lời chính xác không có sẵn, chúng ta có thể (như trong ví dụ so sánh trong phần trước), hãy xem các giá trị tích phân và thấy rằng chúng ổn định khi \ (n \) phát triển. Thật không may, đây là những quy trình thử nghiệm rất yếu và hoàn toàn không thỏa đáng khi tuyên bố rằng phần mềm chúng tôi đã sản xuất được thực hiện chính xác.

Để thấy điều này, chúng ta có thể giới thiệu một lỗi trong hàm

>>> a = 1; b = 2; expected = 3
>>> a + b == expected
True

>>> a = 1; b = 2; expected = 3
>>> a + b == expected
True

Chúng ta hãy nhìn vào một lỗi khác, lần này là trong thuật toán toán học: thay vì tính toán \ (\ frac {1} {2} (f (a) + f (b)) \) như chúng ta nên, chúng ta quên mất thứ hai \ ( \ frac {1} {2} \) và viết

>>> a = 1; b = 2; expected = 3
>>> a + b == expected
True

Kiểm tra đơn vị.

Một thói quen tốt là kiểm tra các mảnh nhỏ của một mã lớn hơn, từng lần một. Điều này được gọi là thử nghiệm đơn vị. Một xác định một đơn vị (nhỏ) của mã, và sau đó một đơn vị thực hiện một thử nghiệm riêng cho đơn vị này. Bài kiểm tra đơn vị phải độc lập theo nghĩa là nó có thể được chạy mà không có kết quả của các bài kiểm tra khác. Thông thường, một thuật toán trong các chương trình khoa học được coi là một đơn vị. Thách thức với các bài kiểm tra đơn vị trong điện toán số là đối phó với các lỗi xấp xỉ bằng số. Một tác dụng phụ may mắn của thử nghiệm đơn vị là lập trình viên buộc phải sử dụng các chức năng để mô đun hóa mã thành các phần nhỏ hơn, logic.

Thủ tục kiểm tra thích hợp

Có ba cách nghiêm túc để kiểm tra việc thực hiện các phương pháp số thông qua các bài kiểm tra đơn vị:

1. So sánh với kết quả tính toán bằng tay trong một vấn đề với một vài hoạt động số học, tức là, nhỏ \ (n \).
2. Giải quyết một vấn đề mà không có lỗi số. Chúng tôi biết rằng quy tắc hình thang phải chính xác cho các hàm tuyến tính. Lỗi do chương trình tạo ra sau đó phải bằng 0 (về độ chính xác của máy).
3. Thể hiện tỷ lệ hội tụ chính xác. Một bài kiểm tra mạnh mẽ khi chúng ta có thể tính toán các lỗi chính xác, là để xem lỗi sẽ tăng nhanh như thế nào khi \ (n \) tăng lên. Trong các quy tắc hình thang và trung điểm, người ta biết rằng lỗi phụ thuộc vào \ (n \) là \ (n^{-2} \) là \ (n \ rightarrow \ infty \).

Kết quả tính toán bằng tay

Chúng ta hãy sử dụng hai hình thang và tính toán tích phân \ (\ int_0^1V (t) \), \ (v (t) = 3t^2e^{t^3} \): $$ \ frac {1} {2} H (V (0) + V (0,5)) + \ frac {1} {2} H (v (0,5) + V (1)) = 2.463642041244344, $$ khi \ (h = 0,5 \) là chiều rộng của Hai hình thang. Chạy chương trình cho kết quả chính xác.

Giải quyết vấn đề mà không có lỗi số

Các bài kiểm tra đơn vị tốt nhất cho các thuật toán số liên quan đến các vấn đề toán học nơi chúng ta biết kết quả số trước. Thông thường, kết quả số có chứa các lỗi xấp xỉ không xác định, vì vậy biết kết quả số ngụ ý rằng chúng ta có vấn đề trong đó các lỗi xấp xỉ biến mất. Tính năng này có thể có mặt trong các vấn đề toán học rất đơn giản. Ví dụ, phương pháp hình thang là chính xác để tích hợp các hàm tuyến tính \ (f (x) = ax+b \). Do đó, chúng ta có thể chọn một số hàm tuyến tính và xây dựng hàm thử nghiệm kiểm tra sự bình đẳng giữa biểu thức phân tích chính xác cho tích phân và số được tính toán bằng cách thực hiện phương pháp hình thang.

Một trường hợp thử nghiệm cụ thể có thể là \ (\ int_ {1.2}^{4.4} (6x-4) dx \). Tích phân này liên quan đến khoảng thời gian "tùy ý" \ ([1.2, 4.4] \) và hàm tuyến tính "tùy ý" \ (f (x) = 6x-4 \). Bằng cách "tùy ý", chúng tôi có nghĩa là các biểu thức trong đó chúng tôi tránh các số đặc biệt 0 và 1 vì chúng có các thuộc tính đặc biệt trong các hoạt động số học (ví dụ: quên nhân tương đương với nhân với 1 và quên thêm tương đương với việc thêm 0).

Thể hiện tỷ lệ hội tụ chính xác

Thông thường, các bài kiểm tra đơn vị phải dựa trên các vấn đề trong đó các lỗi xấp xỉ bằng số trong việc thực hiện của chúng tôi vẫn chưa được biết. Tuy nhiên, chúng ta thường biết hoặc có thể giả định một hành vi tiệm cận nhất định của lỗi. Chúng ta có thể thực hiện một số lần chạy thử nghiệm với bài kiểm tra \ (\ int_0^1 3t^2e^{t^3} dt \) trong đó \ (n \) được nhân đôi trong mỗi lần chạy: \ (n = 4,8,16 \ ). Các lỗi tương ứng sau đó lần lượt là 12%, 3%và 0,77%. Những con số này chỉ ra rằng lỗi giảm gần như là hệ số 4 khi nhân đôi \ (n \). Do đó, lỗi hội tụ về 0 là \ (n^{-2} \) và chúng tôi nói rằng tốc độ hội tụ là 2. Trên thực tế, kết quả này cũng có thể được hiển thị về mặt toán học cho phương pháp hình thang và trung điểm. Các phương thức tích hợp số thường có một lỗi hội tụ về 0 là \ (n^{-p} \) cho một số \ (p \) phụ thuộc vào phương thức. Kết quả như vậy, không có vấn đề gì nếu chúng ta không biết lỗi xấp xỉ thực tế là gì: chúng ta biết với tốc độ nào được giảm vị trí để đo lường tỷ lệ dự kiến ​​và xem nó có đạt được không.

Ý tưởng về một thử nghiệm đơn vị tương ứng sau đó là chạy thuật toán cho một số giá trị \ (n \), tính toán lỗi (giá trị tuyệt đối của chênh lệch giữa kết quả phân tích chính xác và phương pháp được tạo bằng phương pháp số) và kiểm tra xem liệu Lỗi có hành vi tiệm cận xấp xỉ chính xác, tức là, lỗi tỷ lệ thuận với \ (n^{-2} \) trong trường hợp phương pháp hình thang và trung điểm.

Chúng ta hãy phát triển một phương pháp chính xác hơn cho các thử nghiệm đơn vị như vậy dựa trên tỷ lệ hội tụ. Chúng tôi giả sử rằng lỗi \ (e \) phụ thuộc vào \ (n \) theo $$ e = cn^r, $$

trong đó \ (c \) là một hằng số chưa biết và \ (r \) là tốc độ hội tụ. Hãy xem xét một tập hợp các thí nghiệm với các \ (n \) khác nhau: \ (n_0, n_1, n_2, \ ldots, n_q \). Chúng tôi tính toán các lỗi tương ứng \ (e_0, ​​\ ldots, e_q \). Đối với hai thử nghiệm liên tiếp, số \ (i \) và \ (i-1 \), chúng tôi có mô hình lỗi $$ \ start \\ e_ {i-1} & = cn_ {i-1}^r \ thinspace. \ TAG {24} \ end {align} $$ Đây là hai phương trình cho hai chưa biết \ (c \) và \ (r \). Chúng ta có thể dễ dàng loại bỏ \ (c \) bằng cách chia các phương trình cho nhau. Sau đó giải quyết cho \ (r \) cho $$ \ start })} \ Thinspace. \ TAG {25} \ end {Công thức} $$ Chúng tôi đã giới thiệu một chỉ số \ (i-1 \) trong \ (r \) vì giá trị ước tính cho \ (r \) thay đổi theo \ (i \). Hy vọng, \ (r_ {i-1} \) tiếp cận tốc độ hội tụ chính xác khi số khoảng thời gian tăng và \ (i \ rightarrow q \).

Độ chính xác hữu hạn của các số điểm nổi

Các quy trình kiểm tra ở trên dẫn đến so sánh các số để kiểm tra xem các tính toán có chính xác không. So sánh như vậy phức tạp hơn những gì một người mới có thể nghĩ. Giả sử chúng ta có tính toán

>>> a = 1; b = 2; expected = 3
>>> a + b == expected
True

>>> a = 1; b = 2; expected = 3
>>> a + b == expected
True

>>> a = 0.1; b = 0.2; expected = 0.3
>>> a + b == expected
False

Vậy tại sao \ (0,1 + 0,2 \ neq 0,3 \)? Lý do là các số thực nói chung không thể được thể hiện chính xác trên máy tính. Thay vào đó, chúng phải được xấp xỉ bằng một số điểm nổi chỉ có thể lưu trữ một lượng thông tin hữu hạn, thường là khoảng 17 chữ số của một số thực. Hãy để chúng tôi in 0,1, 0,2, 0,1 + 0,2 và 0,3 với 17 số thập phân:

>>> print '%.17f\n%.17f\n%.17f\n%.17f' % (0.1, 0.2, 0.1 + 0.2, 0.3)
0.10000000000000001
0.20000000000000001
0.30000000000000004
0.29999999999999999

Khi chúng tôi tính toán với các số thực, các số này được thể hiện không chính xác trên máy tính và các hoạt động số học với các số không chính xác dẫn đến các lỗi làm tròn nhỏ trong kết quả cuối cùng. Tùy thuộc vào loại thuật toán số, các lỗi làm tròn có thể hoặc không thể tích lũy.

Nếu chúng ta không thể thực hiện các bài kiểm tra như

>>> a = 1; b = 2; expected = 3
>>> a + b == expected
True

>>> a = 0.1; b = 0.2; expected = 0.3
>>> computed = a + b
>>> diff = abs(expected - computed)
>>> tol = 1E-15
>>> diff < tol
True

Ở đây chúng tôi đã thiết lập dung sai để so sánh với \ (10^{-15} \), nhưng tính toán

>>> a = 1; b = 2; expected = 3
>>> a + b == expected
True

>>> a = 1; b = 2; expected = 3
>>> a + b == expected
True

Thực hiện một thử nghiệm với \ (10^k + 0,3 - (10^k + 0,1 + 0,2) \) cho \ (k = 1, \ ldots, 10 \) cho thấy câu trả lời (phải là 0) là xung quanh \ ( 10^{16-K} \). Điều này có nghĩa là dung sai phải lớn hơn nếu chúng ta tính toán với số lượng lớn hơn. Do đó, việc thiết lập một dung sai thích hợp đòi hỏi một số thí nghiệm để xem mức độ chính xác mà người ta có thể mong đợi. Một cách thoát khỏi khó khăn này là làm việc với tương đối thay vì sự khác biệt tuyệt đối. Trong một sự khác biệt tương đối, chúng tôi chia cho một trong các toán hạng, ví dụ: $$ a = 10^k + 0.3, \ Quad B = (10^k + 0.1 + 0.2), \ Quad C = \ frac {A - B} { A} \ Thinspace. $$ Tính toán này \ (c \) cho các \ (k \) khác nhau hiển thị một giá trị xung quanh \ (10^{-16} \). Do đó, một thủ tục an toàn hơn là sử dụng sự khác biệt tương đối.

Xây dựng các bài kiểm tra đơn vị và viết các chức năng kiểm tra

Python có một số khung để tự động chạy và kiểm tra số lượng thử nghiệm rất lớn cho các phần của phần mềm của bạn bằng một lệnh. Đây là một tính năng cực kỳ hữu ích trong quá trình phát triển chương trình: bất cứ khi nào bạn thực hiện một số thay đổi đối với một hoặc nhiều tệp, hãy khởi chạy lệnh kiểm tra và đảm bảo không có gì bị hỏng do chỉnh sửa của bạn.

Các khung thử nghiệm

>>> a = 0.1; b = 0.2; expected = 0.3
>>> a + b == expected
False

>>> a = 0.1; b = 0.2; expected = 0.3
>>> a + b == expected
False

• Tên phải bắt đầu với

  >>> a = 0.1; b = 0.2; expected = 0.3
  >>> a + b == expected
  False
  

  2
• Chức năng kiểm tra không thể có bất kỳ đối số nào
• Các thử nghiệm bên trong các chức năng kiểm tra phải là biểu thức boolean
• Biểu thức boolean

  >>> a = 0.1; b = 0.2; expected = 0.3
  >>> a + b == expected
  False
  

  3 phải được kiểm tra bằng

  >>> a = 0.1; b = 0.2; expected = 0.3
  >>> a + b == expected
  False
  

  4, trong đó

  >>> a = 0.1; b = 0.2; expected = 0.3
  >>> a + b == expected
  False
  

  5 là một đối tượng tùy chọn (chuỗi hoặc số) để được ghi ra khi

  >>> a = 0.1; b = 0.2; expected = 0.3
  >>> a + b == expected
  False
  

  3 là sai

def add(a, b):
    return a + b

def test_add()
    expected = 1 + 1
    computed = add(1, 1)
    assert computed == expected, '1+1=%g' % computed

>>> a = 0.1; b = 0.2; expected = 0.3
>>> a + b == expected
False

>>> a = 0.1; b = 0.2; expected = 0.3
>>> a + b == expected
False

>>> a = 0.1; b = 0.2; expected = 0.3
>>> a + b == expected
False

>>> print '%.17f\n%.17f\n%.17f\n%.17f' % (0.1, 0.2, 0.1 + 0.2, 0.3)
0.10000000000000001
0.20000000000000001
0.30000000000000004
0.29999999999999999

>>> a = 0.1; b = 0.2; expected = 0.3
>>> a + b == expected
False

>>> print '%.17f\n%.17f\n%.17f\n%.17f' % (0.1, 0.2, 0.1 + 0.2, 0.3)
0.10000000000000001
0.20000000000000001
0.30000000000000004
0.29999999999999999

>>> print '%.17f\n%.17f\n%.17f\n%.17f' % (0.1, 0.2, 0.1 + 0.2, 0.3)
0.10000000000000001
0.20000000000000001
0.30000000000000004
0.29999999999999999

Miễn là chúng ta thêm số nguyên, thử nghiệm bình đẳng trong hàm

>>> print '%.17f\n%.17f\n%.17f\n%.17f' % (0.1, 0.2, 0.1 + 0.2, 0.3)
0.10000000000000001
0.20000000000000001
0.30000000000000004
0.29999999999999999

>>> print '%.17f\n%.17f\n%.17f\n%.17f' % (0.1, 0.2, 0.1 + 0.2, 0.3)
0.10000000000000001
0.20000000000000001
0.30000000000000004
0.29999999999999999

def test_add()
    expected = 0.3
    computed = add(0.1, 0.2)
    tol = 1E-14
    diff = abs(expected - computed)
    assert diff < tol, 'diff=%g' % diff

Dưới đây chúng tôi sẽ viết các chức năng kiểm tra cho từng trong ba quy trình kiểm tra mà chúng tôi đề xuất: so sánh với tính toán tay, kiểm tra các vấn đề có thể được giải quyết chính xác và kiểm tra tỷ lệ hội tụ. Chúng tôi gắn bó để kiểm tra mã tích hợp hình thang và thu thập tất cả các chức năng kiểm tra trong một tệp phổ biến bằng tên

>>> print '%.17f\n%.17f\n%.17f\n%.17f' % (0.1, 0.2, 0.1 + 0.2, 0.3)
0.10000000000000001
0.20000000000000001
0.30000000000000004
0.29999999999999999

Kết quả số tính toán bằng tay

Các tính toán tay trước của chúng tôi cho hai hình thang có thể được kiểm tra đối với hàm

>>> a = 1; b = 2; expected = 3
>>> a + b == expected
True

from trapezoidal import trapezoidal
from midpoint import midpoint
from math import exp

g = lambda y: exp(-y**2)
a = 0
b = 2
print '    n        midpoint          trapezoidal'
for i in range(1, 21):
    n = 2**i
    m = midpoint(g, a, b, n)
    t = trapezoidal(g, a, b, n)
    print '%7d %.16f %.16f' % (n, m, t)

>>> print '%.17f\n%.17f\n%.17f\n%.17f' % (0.1, 0.2, 0.1 + 0.2, 0.3)
0.10000000000000001
0.20000000000000001
0.30000000000000004
0.29999999999999999

>>> print '%.17f\n%.17f\n%.17f\n%.17f' % (0.1, 0.2, 0.1 + 0.2, 0.3)
0.10000000000000001
0.20000000000000001
0.30000000000000004
0.29999999999999999

>>> a = 1; b = 2; expected = 3
>>> a + b == expected
True

Giải quyết vấn đề mà không có lỗi số

Chúng tôi biết rằng quy tắc hình thang là chính xác cho các tích phân tuyến tính. Chọn tích phân \ (\ int_ {1.2}^{4.4} (6x-4) dx \) làm trường hợp thử nghiệm, chức năng kiểm tra tương ứng cho bài kiểm tra đơn vị này có thể trông giống như

from trapezoidal import trapezoidal
from midpoint import midpoint
from math import exp

g = lambda y: exp(-y**2)
a = 0
b = 2
print '    n        midpoint          trapezoidal'
for i in range(1, 21):
    n = 2**i
    m = midpoint(g, a, b, n)
    t = trapezoidal(g, a, b, n)
    print '%7d %.16f %.16f' % (n, m, t)

Thể hiện tỷ lệ hội tụ chính xác

Trong ví dụ hiện tại với tích hợp, người ta biết rằng các lỗi xấp xỉ trong quy tắc hình thang tỷ lệ với \ (n^{-2} \), \ (n \) là số lượng phụ được sử dụng trong quy tắc tổng hợp.

Tỷ lệ hội tụ tính toán đòi hỏi lập trình tẻ nhạt hơn một chút so với các thử nghiệm trước đó, nhưng có thể được áp dụng cho các tích phân chung hơn. Thuật toán thường giống như

• cho \ (i = 0,1,2, \ ldots, q \)
• \ (n_i = 2^{i+1} \)
• Tính toán tích phân với khoảng thời gian \ (n_i \)
• Tính toán lỗi \ (e_i \)
• Ước tính \ (r_i \) từ (25) nếu \ (i> 0 \)

from trapezoidal import trapezoidal
from midpoint import midpoint
from math import exp

g = lambda y: exp(-y**2)
a = 0
b = 2
print '    n        midpoint          trapezoidal'
for i in range(1, 21):
    n = 2**i
    m = midpoint(g, a, b, n)
    t = trapezoidal(g, a, b, n)
    print '%7d %.16f %.16f' % (n, m, t)

Tạo chức năng kiểm tra là vấn đề chọn

>>> a = 0.1; b = 0.2; expected = 0.3
>>> computed = a + b
>>> diff = abs(expected - computed)
>>> tol = 1E-15
>>> diff < tol
True

>>> a = 0.1; b = 0.2; expected = 0.3
>>> computed = a + b
>>> diff = abs(expected - computed)
>>> tol = 1E-15
>>> diff < tol
True

>>> a = 0.1; b = 0.2; expected = 0.3
>>> computed = a + b
>>> diff = abs(expected - computed)
>>> tol = 1E-15
>>> diff < tol
True

>>> a = 0.1; b = 0.2; expected = 0.3
>>> a + b == expected
False

from trapezoidal import trapezoidal
from midpoint import midpoint
from math import exp

g = lambda y: exp(-y**2)
a = 0
b = 2
print '    n        midpoint          trapezoidal'
for i in range(1, 21):
    n = 2**i
    m = midpoint(g, a, b, n)
    t = trapezoidal(g, a, b, n)
    print '%7d %.16f %.16f' % (n, m, t)

Chạy thử nghiệm cho thấy tất cả \ (r_i \), ngoại trừ lần đầu tiên, bằng giới hạn mục tiêu 2 trong vòng hai số thập phân. Quan sát này cho thấy sự dung nạp của \ (10^{-2} \).

Ghi chú về kiểm soát phiên bản của các tập tin.

Có một bộ chức năng kiểm tra để tự động kiểm tra xem phần mềm của bạn có hoạt động được coi là một yêu cầu cơ bản cho điện toán đáng tin cậy. Quan trọng không kém là một hệ thống có thể theo dõi các phiên bản khác nhau của các tệp và các thử nghiệm, được gọi là hệ thống điều khiển phiên bản. Hệ thống kiểm soát phiên bản phổ biến nhất ngày nay là Git, mà các tác giả khuyên người đọc sử dụng cho các báo cáo lập trình và viết. Sự kết hợp của Git và lưu trữ đám mây như GitHub là một cách rất phổ biến để tổ chức công việc khoa học hoặc kỹ thuật. Chúng tôi có một đoạn giới thiệu nhanh chóng để Git và GitHub giúp bạn đứng dậy trong vòng vài phút.

Quy trình làm việc điển hình với Git đi như sau.

1. Trước khi bạn bắt đầu làm việc với các tệp, hãy đảm bảo bạn có phiên bản mới nhất của chúng bằng cách chạy

   >>> a = 0.1; b = 0.2; expected = 0.3
   >>> computed = a + b
   >>> diff = abs(expected - computed)
   >>> tol = 1E-15
   >>> diff < tol
   True
   

   5.
2. Chỉnh sửa các tệp, xóa hoặc tạo tệp (phải đăng ký các tệp mới bằng

   >>> a = 0.1; b = 0.2; expected = 0.3
   >>> computed = a + b
   >>> diff = abs(expected - computed)
   >>> tol = 1E-15
   >>> diff < tol
   True
   

   6).
3. Khi một công việc tự nhiên được thực hiện, hãy thực hiện các thay đổi của bạn bằng lệnh

   >>> a = 0.1; b = 0.2; expected = 0.3
   >>> computed = a + b
   >>> diff = abs(expected - computed)
   >>> tol = 1E-15
   >>> diff < tol
   True
   

   7.
4. Thực hiện các thay đổi của bạn cũng trong đám mây bằng cách thực hiện

   >>> a = 0.1; b = 0.2; expected = 0.3
   >>> computed = a + b
   >>> diff = abs(expected - computed)
   >>> tol = 1E-15
   >>> diff < tol
   True
   

   8.

Một tính năng hay của Git là mọi người có thể chỉnh sửa cùng một tệp cùng một lúc và rất thường xuyên Git sẽ có thể tự động hợp nhất các thay đổi (!). Do đó, kiểm soát phiên bản là rất quan trọng khi bạn làm việc với người khác hoặc khi bạn thực hiện công việc của mình trên các loại máy tính khác nhau. Một tính năng quan trọng khác là bất cứ ai cũng có thể xem lịch sử của một tệp, xem ai đã làm gì khi nào và quay lại toàn bộ bộ sưu tập tệp thành một cam kết trước đó. Tính năng này, tất nhiên, là cơ bản cho công việc đáng tin cậy.

Vector hóa

Các hàm

>>> a = 1; b = 2; expected = 3
>>> a + b == expected
True

def add(a, b):
    return a + b

We have already seen simple examples on vectorization in the section A Python program with vectorization and plotting when we could evaluate a mathematical function \( f(x) \) for a large number of \( x \) values stored in an array. Basically, we can write

from trapezoidal import trapezoidal
from midpoint import midpoint
from math import exp

g = lambda y: exp(-y**2)
a = 0
b = 2
print '    n        midpoint          trapezoidal'
for i in range(1, 21):
    n = 2**i
    m = midpoint(g, a, b, n)
    t = trapezoidal(g, a, b, n)
    print '%7d %.16f %.16f' % (n, m, t)

def add(a, b):
    return a + b

def add(a, b):
    return a + b

def add(a, b):
    return a + b

>>> a = 0.1; b = 0.2; expected = 0.3
>>> computed = a + b
>>> diff = abs(expected - computed)
>>> tol = 1E-15
>>> diff < tol
True

def add(a, b):
    return a + b

>>> a = 0.1; b = 0.2; expected = 0.3
>>> computed = a + b
>>> diff = abs(expected - computed)
>>> tol = 1E-15
>>> diff < tol
True

def add(a, b):
    return a + b

Vectorizing the midpoint rule

The aim of vectorizing the

>>> a = 1; b = 2; expected = 3
>>> a + b == expected
True

>>> a = 1; b = 2; expected = 3
>>> a + b == expected
True

>>> a = 1; b = 2; expected = 3
>>> a + b == expected
True

def add(a, b):
    return a + b

1. compute all the evaluation points in one array

   def add(a, b):
       return a + b
   

   3
2. call

   def test_add()
       expected = 1 + 1
       computed = add(1, 1)
       assert computed == expected, '1+1=%g' % computed
   

   3 to produce an array of corresponding function values
3. use the

   def test_add()
       expected = 1 + 1
       computed = add(1, 1)
       assert computed == expected, '1+1=%g' % computed
   

   4 function to sum the

   def test_add()
       expected = 1 + 1
       computed = add(1, 1)
       assert computed == expected, '1+1=%g' % computed
   

   3 values

def test_add()
    expected = 1 + 1
    computed = add(1, 1)
    assert computed == expected, '1+1=%g' % computed

>>> a = 0.1; b = 0.2; expected = 0.3
>>> computed = a + b
>>> diff = abs(expected - computed)
>>> tol = 1E-15
>>> diff < tol
True

def test_add()
    expected = 1 + 1
    computed = add(1, 1)
    assert computed == expected, '1+1=%g' % computed

def test_add()
    expected = 1 + 1
    computed = add(1, 1)
    assert computed == expected, '1+1=%g' % computed

def test_add()
    expected = 0.3
    computed = add(0.1, 0.2)
    tol = 1E-14
    diff = abs(expected - computed)
    assert diff < tol, 'diff=%g' % diff

from trapezoidal import trapezoidal
from midpoint import midpoint
from math import exp

g = lambda y: exp(-y**2)
a = 0
b = 2
print '    n        midpoint          trapezoidal'
for i in range(1, 21):
    n = 2**i
    m = midpoint(g, a, b, n)
    t = trapezoidal(g, a, b, n)
    print '%7d %.16f %.16f' % (n, m, t)

Let us test the code interactively in a Python shell to compute \( \int_0^1 3t^2dt \). The file with the code above has the name

def test_add()
    expected = 0.3
    computed = add(0.1, 0.2)
    tol = 1E-14
    diff = abs(expected - computed)
    assert diff < tol, 'diff=%g' % diff

from trapezoidal import trapezoidal
from midpoint import midpoint
from math import exp

g = lambda y: exp(-y**2)
a = 0
b = 2
print '    n        midpoint          trapezoidal'
for i in range(1, 21):
    n = 2**i
    m = midpoint(g, a, b, n)
    t = trapezoidal(g, a, b, n)
    print '%7d %.16f %.16f' % (n, m, t)

def test_add()
    expected = 0.3
    computed = add(0.1, 0.2)
    tol = 1E-14
    diff = abs(expected - computed)
    assert diff < tol, 'diff=%g' % diff

def test_add()
    expected = 0.3
    computed = add(0.1, 0.2)
    tol = 1E-14
    diff = abs(expected - computed)
    assert diff < tol, 'diff=%g' % diff

def test_add()
    expected = 0.3
    computed = add(0.1, 0.2)
    tol = 1E-14
    diff = abs(expected - computed)
    assert diff < tol, 'diff=%g' % diff

def add(a, b):
    return a + b

def test_add()
    expected = 0.3
    computed = add(0.1, 0.2)
    tol = 1E-14
    diff = abs(expected - computed)
    assert diff < tol, 'diff=%g' % diff

The vectorized code performs all loops very efficiently in compiled code, resulting in much faster execution. Moreover, many readers of the code will also say that the algorithm looks clearer than in the loop-based implementation.

Vectorizing the trapezoidal rule

We can use the same approach to vectorize the

def add(a, b):
    return a + b

def test_add()
    expected = 0.3
    computed = add(0.1, 0.2)
    tol = 1E-14
    diff = abs(expected - computed)
    assert diff < tol, 'diff=%g' % diff

def test_add()
    expected = 0.3
    computed = add(0.1, 0.2)
    tol = 1E-14
    diff = abs(expected - computed)
    assert diff < tol, 'diff=%g' % diff

from trapezoidal import trapezoidal
from midpoint import midpoint
from math import exp

g = lambda y: exp(-y**2)
a = 0
b = 2
print '    n        midpoint          trapezoidal'
for i in range(1, 21):
    n = 2**i
    m = midpoint(g, a, b, n)
    t = trapezoidal(g, a, b, n)
    print '%7d %.16f %.16f' % (n, m, t)

def test_add()
    expected = 0.3
    computed = add(0.1, 0.2)
    tol = 1E-14
    diff = abs(expected - computed)
    assert diff < tol, 'diff=%g' % diff

from trapezoidal import trapezoidal
from midpoint import midpoint
from math import exp

g = lambda y: exp(-y**2)
a = 0
b = 2
print '    n        midpoint          trapezoidal'
for i in range(1, 21):
    n = 2**i
    m = midpoint(g, a, b, n)
    t = trapezoidal(g, a, b, n)
    print '%7d %.16f %.16f' % (n, m, t)

from trapezoidal import trapezoidal
from midpoint import midpoint
from math import exp

g = lambda y: exp(-y**2)
a = 0
b = 2
print '    n        midpoint          trapezoidal'
for i in range(1, 21):
    n = 2**i
    m = midpoint(g, a, b, n)
    t = trapezoidal(g, a, b, n)
    print '%7d %.16f %.16f' % (n, m, t)

Measuring computational speed

Now that we have created faster, vectorized versions of functions in the previous section, it is interesting to measure how much faster they are. The purpose of the present section is therefore to explain how we can record the CPU time consumed by a function so we can answer this question. There are many techniques for measuring the CPU time in Python, and here we shall just explain the simplest and most convenient one: the

from trapezoidal import trapezoidal
from midpoint import midpoint
from math import exp

g = lambda y: exp(-y**2)
a = 0
b = 2
print '    n        midpoint          trapezoidal'
for i in range(1, 21):
    n = 2**i
    m = midpoint(g, a, b, n)
    t = trapezoidal(g, a, b, n)
    print '%7d %.16f %.16f' % (n, m, t)

from trapezoidal import trapezoidal
from midpoint import midpoint
from math import exp

g = lambda y: exp(-y**2)
a = 0
b = 2
print '    n        midpoint          trapezoidal'
for i in range(1, 21):
    n = 2**i
    m = midpoint(g, a, b, n)
    t = trapezoidal(g, a, b, n)
    print '%7d %.16f %.16f' % (n, m, t)

Double and triple integrals

The midpoint rule for a double integral

Given a double integral over a rectangular domain \( [a,b]\times [c,d] \), $$ \int_a^b \int_c^d f(x,y) dydx,$$ how can we approximate this integral by numerical methods?

Derivation via one-dimensional integrals

Since we know how to deal with integrals in one variable, a fruitful approach is to view the double integral as two integrals, each in one variable, which can be approximated numerically by previous one-dimensional formulas. To this end, we introduce a help function \( g(x) \) and write $$ \int_a^b \int_c^d f(x,y) dydx = \int_a^b g(x)dx,\quad g(x) = \int_c^d f(x,y) dy\thinspace .$$ Each of the integrals $$ \int_a^b g(x)dx,\quad g(x) = \int_c^d f(x,y) dy$$ can be discretized by any numerical integration rule for an integral in one variable. Let us use the midpoint method (22) and start with \( g(x)=\int_c^d f(x,y)dy \). We introduce \( n_y \) intervals on \( [c,d] \) with length \( h_y \). The midpoint rule for this integral then becomes $$ g(x) = \int_c^d f(x,y) dy \approx h_y \sum_{j=0}^{n_y-1} f(x,y_j), \quad y_j = c + \frac{1}{2}{h_y} + jh_y \thinspace . $$ The expression looks somewhat different from (22), but that is because of the notation: since we integrate in the \( y \) direction and will have to work with both \( x \) and \( y \) as coordinates, we must use \( n_y \) for \( n \), \( h_y \) for \( h \), and the counter \( i \) is more naturally called \( j \) when integrating in \( y \). Integrals in the \( x \) direction will use \( h_x \) and \( n_x \) for \( h \) and \( n \), and \( i \) as counter.

Tích phân kép là \ (\ int_a^b g (x) dx \), có thể được xấp xỉ bằng phương pháp điểm giữa: $$ \ int_a^b g (x) dx \ xấp xỉ h_x \ sum_ {i = 0}^ 1} g (x_i), \ Quad x_i = a + \ frac {1} {2} {h_x} + ih_x \ thinspace. \ Bắt đầu {Align} \ int_a^b \ int_c^d f (x, y) dydx & \ xấp xỉ h_x \ sum_ {i = 0}^{n_x-1} h_y \ sum_ {j = 0}^{n_y-1} f (x_i, y_j) \ nonumber \\ & = h_xh_y \ sum_ {i = 0}^{n_x-1} \ sum_ {j = 0}^{n_y-1} f (a + \ frac {h_x} {2 } + ih_x, c + \ frac {h_y} {2} + jh_y) \ thinspace. \ tag {26} \ end {align} $$

Đạo hàm trực tiếp

Công thức (26) cũng có thể được lấy trực tiếp trong trường hợp hai chiều bằng cách áp dụng ý tưởng của phương pháp điểm giữa. Chúng tôi chia hình chữ nhật \ ([A, B] \ Times [C, D] \) thành \ (N_X \ Times N_Y \) Các ô có kích thước bằng nhau. Ý tưởng về phương pháp trung điểm là gần đúng \ (f \) bằng một hằng số trên mỗi ô và đánh giá hằng số tại điểm giữa. Ô \ ((i, j) \) chiếm diện tích $$ [a+ih_x, a+(i+1) h_x] \ lần [c+jh_y, c+(j+1) h_y] \ ((x_i, y_j) \) với $$ x_i = a + ih_x + \ frac {1} {2} {h_x}, \ Quad Y_J = C + JH_Y ​​+ \ frac {1} Do đó, không gian.

Lập trình một tổng số

Công thức (26) liên quan đến một tổng số, thường được thực hiện dưới dạng gấp đôi cho vòng lặp. Triển khai chức năng Python (26) có thể trông giống như

from trapezoidal import trapezoidal
from midpoint import midpoint
from math import exp

g = lambda y: exp(-y**2)
a = 0
b = 2
print '    n        midpoint          trapezoidal'
for i in range(1, 21):
    n = 2**i
    m = midpoint(g, a, b, n)
    t = trapezoidal(g, a, b, n)
    print '%7d %.16f %.16f' % (n, m, t)

    n        midpoint          trapezoidal
      2 0.8842000076332692 0.8770372606158094
      4 0.8827889485397279 0.8806186341245393
      8 0.8822686991994210 0.8817037913321336
     16 0.8821288703366458 0.8819862452657772
     32 0.8820933014203766 0.8820575578012112
     64 0.8820843709743319 0.8820754296107942
    128 0.8820821359746071 0.8820799002925637
    256 0.8820815770754198 0.8820810181335849
    512 0.8820814373412922 0.8820812976045025
   1024 0.8820814024071774 0.8820813674728968
   2048 0.8820813936736116 0.8820813849400392
   4096 0.8820813914902204 0.8820813893068272
   8192 0.8820813909443684 0.8820813903985197
  16384 0.8820813908079066 0.8820813906714446
  32768 0.8820813907737911 0.8820813907396778
  65536 0.8820813907652575 0.8820813907567422
 131072 0.8820813907631487 0.8820813907610036
 262144 0.8820813907625702 0.8820813907620528
 524288 0.8820813907624605 0.8820813907623183
1048576 0.8820813907624268 0.8820813907623890

Tái sử dụng mã cho các tích phân một chiều

Việc viết một phương pháp giữa hai chiều là rất tự nhiên như chúng ta đã làm trong chức năng ____1044 khi chúng ta có công thức (26). Tuy nhiên, chúng tôi có thể hỏi, nhiều như chúng tôi đã làm trong toán học, chúng ta có thể sử dụng lại một triển khai được thử nghiệm tốt cho các tích phân một chiều để tính toán các tích phân kép không? Nghĩa là, chúng ta có thể sử dụng chức năng

>>> a = 1; b = 2; expected = 3
>>> a + b == expected
True

def midpoint(f, a, b, n):
    h = float(b-a)/n
    result = 0
    for i in range(n):
        result += f((a + h/2.0) + i*h)
    result *= h
    return result

    n        midpoint          trapezoidal
      2 0.8842000076332692 0.8770372606158094
      4 0.8827889485397279 0.8806186341245393
      8 0.8822686991994210 0.8817037913321336
     16 0.8821288703366458 0.8819862452657772
     32 0.8820933014203766 0.8820575578012112
     64 0.8820843709743319 0.8820754296107942
    128 0.8820821359746071 0.8820799002925637
    256 0.8820815770754198 0.8820810181335849
    512 0.8820814373412922 0.8820812976045025
   1024 0.8820814024071774 0.8820813674728968
   2048 0.8820813936736116 0.8820813849400392
   4096 0.8820813914902204 0.8820813893068272
   8192 0.8820813909443684 0.8820813903985197
  16384 0.8820813908079066 0.8820813906714446
  32768 0.8820813907737911 0.8820813907396778
  65536 0.8820813907652575 0.8820813907567422
 131072 0.8820813907631487 0.8820813907610036
 262144 0.8820813907625702 0.8820813907620528
 524288 0.8820813907624605 0.8820813907623183
1048576 0.8820813907624268 0.8820813907623890

Xác minh thông qua các chức năng kiểm tra

Làm thế nào chúng ta có thể kiểm tra rằng các chức năng của chúng ta cho công việc tích phân kép? Bài kiểm tra đơn vị tốt nhất là tìm ra một vấn đề trong đó lỗi xấp xỉ bằng số biến mất vì sau đó chúng ta biết chính xác câu trả lời số nên là gì. Quy tắc trung điểm chính xác cho các hàm tuyến tính, bất kể chúng ta sử dụng bao nhiêu con. Ngoài ra, bất kỳ hàm hai chiều tuyến tính nào \ (f (x, y) = px+qy+r \) sẽ được tích hợp chính xác bởi quy tắc trung điểm hai chiều. Chúng tôi có thể chọn \ (f (x, y) = 2x+y \) và tạo chức năng kiểm tra thích hợp có thể tự động xác minh hai triển khai thay thế của chúng tôi về quy tắc điểm giữa hai chiều. Để tính toán tích phân của \ (f (x, y) \), chúng tôi tận dụng Sympy để loại bỏ khả năng của các lỗi trong tính toán tay. Chức năng kiểm tra trở thành

    n        midpoint          trapezoidal
      2 0.8842000076332692 0.8770372606158094
      4 0.8827889485397279 0.8806186341245393
      8 0.8822686991994210 0.8817037913321336
     16 0.8821288703366458 0.8819862452657772
     32 0.8820933014203766 0.8820575578012112
     64 0.8820843709743319 0.8820754296107942
    128 0.8820821359746071 0.8820799002925637
    256 0.8820815770754198 0.8820810181335849
    512 0.8820814373412922 0.8820812976045025
   1024 0.8820814024071774 0.8820813674728968
   2048 0.8820813936736116 0.8820813849400392
   4096 0.8820813914902204 0.8820813893068272
   8192 0.8820813909443684 0.8820813903985197
  16384 0.8820813908079066 0.8820813906714446
  32768 0.8820813907737911 0.8820813907396778
  65536 0.8820813907652575 0.8820813907567422
 131072 0.8820813907631487 0.8820813907610036
 262144 0.8820813907625702 0.8820813907620528
 524288 0.8820813907624605 0.8820813907623183
1048576 0.8820813907624268 0.8820813907623890

Hãy để các chức năng kiểm tra lên tiếng?

Nếu chúng ta gọi hàm

from trapezoidal import trapezoidal
from midpoint import midpoint
from math import exp

g = lambda y: exp(-y**2)
a = 0
b = 2
print '    n        midpoint          trapezoidal'
for i in range(1, 21):
    n = 2**i
    m = midpoint(g, a, b, n)
    t = trapezoidal(g, a, b, n)
    print '%7d %.16f %.16f' % (n, m, t)

Phương pháp hình thang có thể được sử dụng thay thế cho phương pháp trung điểm. Việc tạo ra một công thức cho tích phân kép và việc triển khai tuân theo chính xác các ý tưởng giống như chúng tôi đã giải thích với phương thức điểm giữa, nhưng có nhiều thuật ngữ hơn để viết trong các công thức. Bài tập 42: Xuất phát quy tắc hình thang cho một tích phân kép yêu cầu bạn thực hiện các chi tiết. Bài tập đó là một bài kiểm tra rất tốt về sự hiểu biết của bạn về các ý tưởng toán học và lập trình trong phần hiện tại.

Quy tắc điểm giữa cho một tích phân ba

Học thuyết

Khi một phương pháp hoạt động cho một vấn đề một chiều được khái quát thành hai chiều, thường việc mở rộng phương pháp này thành ba chiều là khá đơn giản. Điều này bây giờ sẽ được chứng minh cho các tích phân. Chúng tôi có tích phân ba $$ \ int_ {a}^{b} \ int_c^d \ int_e^f g (x, y, z) dzdydx $$ và muốn xấp xỉ tích phân theo quy tắc điểm giữa. Theo các ý tưởng cho tích phân kép, chúng tôi chia tích phân này thành tích phân một chiều: $$ \ BEGIN {ALIGN*} P (x, y) & = \ int_e^f g (x, y, z) dz \\ q ( x) & = \ int_c^d p (x, y) dy \\ \ int_ {a}^{b} \ int_c^d \ int_e^f g dx \ end {align*} $$ Đối với mỗi tích phân một chiều này, chúng tôi áp dụng quy tắc điểm giữa: $$ \ start & \ xấp xỉ \ sum_ {k = 0}^{n_z-1} g (x, y, z_k), \\ q (x) = \ int_c^d p (x, y) dy & \ xấp xỉ \ sum_ {j = 0}^{n_y-1} p (x, y_j), \\ \ int_ {a}^{b} \ int_c^d \ int_e^f g (x, y, z) dx & \ xấp xỉ \ sum_ {i = 0}^{n_x-1} q (x_i), \ end {align*} $$ trong đó $$ z_k = e + \ frac {1} {2} h_z + kh_z, \ Quad y_j = C + \ frac {1} {2} H_Y + JH_Y ​​\ Quad X_i = A + \ Frac {1} {2} H_X + IH_X \ Thinspace. $$ Bắt đầu với công thức cho \ (\ int_ {a}^{b} \ int_c^d \ int_e^f g (x, y, z) dzdydx \) và chèn hai công thức trước đó cho $$ \ bắt đầu & \ int_ {a}^{b} \ int_c^d \ int_e^f g (x, y, z) \, dzdydx \ xấp xỉ sum_ {j = 0}^{n_y-1} \ sum_ {k = 0}^{n_z-1} g (a + \ frac {1} {2} h_x + ih_x, c + \ frac {1} {2 } h_y + jh_y, e + \ frac {1} {2} h_z + kh_z) \ thinspace. \ TAG {27} \ end {align} $$ Lưu ý rằng chúng tôi có thể áp dụng các ý tưởng theo đạo hàm trực tiếp ở cuối phần Quy tắc điểm giữa cho một tích phân kép để đến (27) trực tiếp: Chia miền thành \ (N_X \ Times N_y \ Times N_Z \) của các tập \ (h_xh_yh_z \); xấp xỉ \ (g \) bởi một hằng số, được đánh giá tại điểm giữa \ ((x_i, y_j, z_k) \), trong mỗi ô; và tổng hợp các tích phân ô \ (h_xh_yh_zg (x_i, y_j, z_k) \).

Thực hiện

Chúng tôi tuân theo các ý tưởng để thực hiện quy tắc trung điểm cho một tích phân kép. Các chức năng tương ứng được hiển thị bên dưới và được tìm thấy trong tệp midpoint_triple.py.

    n        midpoint          trapezoidal
      2 0.8842000076332692 0.8770372606158094
      4 0.8827889485397279 0.8806186341245393
      8 0.8822686991994210 0.8817037913321336
     16 0.8821288703366458 0.8819862452657772
     32 0.8820933014203766 0.8820575578012112
     64 0.8820843709743319 0.8820754296107942
    128 0.8820821359746071 0.8820799002925637
    256 0.8820815770754198 0.8820810181335849
    512 0.8820814373412922 0.8820812976045025
   1024 0.8820814024071774 0.8820813674728968
   2048 0.8820813936736116 0.8820813849400392
   4096 0.8820813914902204 0.8820813893068272
   8192 0.8820813909443684 0.8820813903985197
  16384 0.8820813908079066 0.8820813906714446
  32768 0.8820813907737911 0.8820813907396778
  65536 0.8820813907652575 0.8820813907567422
 131072 0.8820813907631487 0.8820813907610036
 262144 0.8820813907625702 0.8820813907620528
 524288 0.8820813907624605 0.8820813907623183
1048576 0.8820813907624268 0.8820813907623890

Tích hợp Monte Carlo cho các miền hình phức tạp

Việc sử dụng lặp đi lặp lại các quy tắc tích hợp một chiều để xử lý các tích phân gấp đôi và ba tạo thành một chiến lược hoạt động chỉ khi miền tích hợp là một hình chữ nhật hoặc hộp. Đối với bất kỳ hình dạng nào khác của miền, các phương pháp hoàn toàn khác nhau phải được sử dụng. Một cách tiếp cận phổ biến cho các miền hai và ba chiều là chia miền thành nhiều hình tam giác nhỏ hoặc tứ diện và sử dụng các phương pháp tích hợp số cho mỗi tam giác hoặc tứ diện. Thuật toán tổng thể và triển khai quá phức tạp để được giải quyết trong cuốn sách này. Thay vào đó, chúng ta sẽ sử dụng một phương pháp thay thế, rất đơn giản và chung chung, được gọi là tích hợp Monte Carlo. Nó có thể được thực hiện trong một nửa trang mã, nhưng yêu cầu các đơn đặt hàng có cường độ đánh giá chức năng nhiều hơn trong các tích phân kép so với quy tắc điểm giữa.

Tuy nhiên, tích hợp Monte Carlo có hiệu quả tính toán hơn nhiều so với quy tắc trung điểm khi tính toán các tích phân chiều cao hơn trong hơn ba biến trên các miền HyperCube. Ý tưởng của chúng tôi cho các tích phân gấp đôi và ba có thể dễ dàng được khái quát hóa để xử lý một tích phân trong các biến \ (m \). Một công thức trung điểm sau đó liên quan đến tổng số \ (m \). Với các ô \ (n \) theo mỗi hướng tọa độ, công thức yêu cầu đánh giá chức năng \ (n^m \). Đó là, công việc tính toán bùng nổ như một hàm số mũ của số lượng kích thước không gian. Mặt khác, tích hợp Monte Carlo không chịu sự bùng nổ của công việc tính toán này và là phương pháp ưa thích để tính toán các tích phân chiều cao hơn. Vì vậy, nó có ý nghĩa trong một chương về tích hợp số để giải quyết các phương pháp Monte Carlo, cả để xử lý các miền phức tạp và để xử lý tích phân với nhiều biến.

Thuật toán tích hợp Monte Carlo

Ý tưởng về tích hợp Monte Carlo của \ (\ int_a^b f (x) dx \) là sử dụng định lý giá trị trung bình từ tính toán, trong đó nói rằng tích phân \ (\ int_a^b f (x) dx \) của miền tích hợp, ở đây \ (b-a \), nhân thời gian giá trị trung bình của \ (f \), \ (\ bar f \), trong \ ([a, b] \). Giá trị trung bình có thể được tính bằng cách lấy mẫu \ (f \) tại một tập hợp các điểm ngẫu nhiên bên trong miền và lấy giá trị trung bình của các giá trị hàm. Trong các kích thước cao hơn, một tích phân được ước tính là diện tích/khối lượng của miền thời gian giá trị trung bình và một lần nữa người ta có thể đánh giá tích phân tại một tập hợp các điểm ngẫu nhiên trong miền và tính giá trị trung bình của các đánh giá đó.

Let us introduce some quantities to help us make the specification of the integration algorithm more precise. Suppose we have some two-dimensional integral $$ \int_\Omega f(x,y)dxdx, $$ where \( \Omega \) is a two-dimensional domain defined via a help function \( g(x,y) \): $$ \Omega = \{ (x,y)\,|\, g(x,y) \geq 0\} $$ That is, points \( (x,y) \) for which \( g(x,y)\geq 0 \) lie inside \( \Omega \), and points for which \( g(x,y) < \Omega \) are outside \( \Omega \). The boundary of the domain \( \partial\Omega \) is given by the implicit curve \( g(x,y)=0 \). Such formulations of geometries have been very common during the last couple of decades, and one refers to \( g \) as a level-set function and the boundary \( g=0 \) as the zero-level contour of the level-set function. For simple geometries one can easily construct \( g \) by hand, while in more complicated industrial applications one must resort to mathematical models for constructing \( g \).

Let \( A(\Omega) \) be the area of a domain \( \Omega \). We can estimate the integral by this Monte Carlo integration method:

1. embed the geometry \( \Omega \) in a rectangular area \( R \)
2. draw a large number of random points \( (x,y) \) in \( R \)
3. count the fraction \( q \) of points that are inside \( \Omega \)
4. approximate \( A(\Omega)/A(R) \) by \( q \), i.e., set \( A(\Omega) = qA(R) \)
5. evaluate the mean of \( f \), \( \bar f \), at the points inside \( \Omega \)
6. estimate the integral as \( A(\Omega)\bar f \)

To get an idea of the method, consider a circular domain \( \Omega \) embedded in a rectangle as shown below. A collection of random points is illustrated by black dots.

Implementation

A Python function implementing \( \int_\Omega f(x,y)dxdy \) can be written like this:

    n        midpoint          trapezoidal
      2 0.8842000076332692 0.8770372606158094
      4 0.8827889485397279 0.8806186341245393
      8 0.8822686991994210 0.8817037913321336
     16 0.8821288703366458 0.8819862452657772
     32 0.8820933014203766 0.8820575578012112
     64 0.8820843709743319 0.8820754296107942
    128 0.8820821359746071 0.8820799002925637
    256 0.8820815770754198 0.8820810181335849
    512 0.8820814373412922 0.8820812976045025
   1024 0.8820814024071774 0.8820813674728968
   2048 0.8820813936736116 0.8820813849400392
   4096 0.8820813914902204 0.8820813893068272
   8192 0.8820813909443684 0.8820813903985197
  16384 0.8820813908079066 0.8820813906714446
  32768 0.8820813907737911 0.8820813907396778
  65536 0.8820813907652575 0.8820813907567422
 131072 0.8820813907631487 0.8820813907610036
 262144 0.8820813907625702 0.8820813907620528
 524288 0.8820813907624605 0.8820813907623183
1048576 0.8820813907624268 0.8820813907623890

Verification

A simple test case is to check the area of a rectangle \( [0,2]\times[3,4.5] \) embedded in a rectangle \( [0,3]\times [2,5] \). The right answer is 3, but Monte Carlo integration is, unfortunately, never exact so it is impossible to predict the output of the algorithm. All we know is that the estimated integral should approach 3 as the number of random points goes to infinity. Also, for a fixed number of points, we can run the algorithm several times and get different numbers that fluctuate around the exact value, since different sample points are used in different calls to the Monte Carlo integration algorithm.

The area of the rectangle can be computed by the integral \( \int_0^2\int_3^{4.5} dydx \), so in this case we identify \( f(x,y)=1 \), and the \( g \) function can be specified as (e.g.) 1 if \( (x,y) \) is inside \( [0,2]\times[3,4.5] \) and \( -1 \) otherwise. Here is an example on how we can utilize the

from trapezoidal import trapezoidal
from midpoint import midpoint
from math import exp

g = lambda y: exp(-y**2)
a = 0
b = 2
print '    n        midpoint          trapezoidal'
for i in range(1, 21):
    n = 2**i
    m = midpoint(g, a, b, n)
    t = trapezoidal(g, a, b, n)
    print '%7d %.16f %.16f' % (n, m, t)

    n        midpoint          trapezoidal
      2 0.8842000076332692 0.8770372606158094
      4 0.8827889485397279 0.8806186341245393
      8 0.8822686991994210 0.8817037913321336
     16 0.8821288703366458 0.8819862452657772
     32 0.8820933014203766 0.8820575578012112
     64 0.8820843709743319 0.8820754296107942
    128 0.8820821359746071 0.8820799002925637
    256 0.8820815770754198 0.8820810181335849
    512 0.8820814373412922 0.8820812976045025
   1024 0.8820814024071774 0.8820813674728968
   2048 0.8820813936736116 0.8820813849400392
   4096 0.8820813914902204 0.8820813893068272
   8192 0.8820813909443684 0.8820813903985197
  16384 0.8820813908079066 0.8820813906714446
  32768 0.8820813907737911 0.8820813907396778
  65536 0.8820813907652575 0.8820813907567422
 131072 0.8820813907631487 0.8820813907610036
 262144 0.8820813907625702 0.8820813907620528
 524288 0.8820813907624605 0.8820813907623183
1048576 0.8820813907624268 0.8820813907623890

It is mathematically known that the standard deviation of the Monte Carlo estimate of an integral converges as \( n^{-1/2} \), where \( n \) is the number of samples. This kind of convergence rate estimate could be used to verify the implementation, but this topic is beyond the scope of this book.

Test function for function with random numbers

To make a test function, we need a unit test that has identical behavior each time we run the test. This seems difficult when random numbers are involved, because these numbers are different every time we run the algorithm, and each run hence produces a (slightly) different result. A standard way to test algorithms involving random numbers is to fix the seed of the random number generator. Then the sequence of numbers is the same every time we run the algorithm. Assuming that the

from trapezoidal import trapezoidal
from midpoint import midpoint
from math import exp

g = lambda y: exp(-y**2)
a = 0
b = 2
print '    n        midpoint          trapezoidal'
for i in range(1, 21):
    n = 2**i
    m = midpoint(g, a, b, n)
    t = trapezoidal(g, a, b, n)
    print '%7d %.16f %.16f' % (n, m, t)

from trapezoidal import trapezoidal
from midpoint import midpoint
from math import exp

g = lambda y: exp(-y**2)
a = 0
b = 2
print '    n        midpoint          trapezoidal'
for i in range(1, 21):
    n = 2**i
    m = midpoint(g, a, b, n)
    t = trapezoidal(g, a, b, n)
    print '%7d %.16f %.16f' % (n, m, t)

    n        midpoint          trapezoidal
      2 0.8842000076332692 0.8770372606158094
      4 0.8827889485397279 0.8806186341245393
      8 0.8822686991994210 0.8817037913321336
     16 0.8821288703366458 0.8819862452657772
     32 0.8820933014203766 0.8820575578012112
     64 0.8820843709743319 0.8820754296107942
    128 0.8820821359746071 0.8820799002925637
    256 0.8820815770754198 0.8820810181335849
    512 0.8820814373412922 0.8820812976045025
   1024 0.8820814024071774 0.8820813674728968
   2048 0.8820813936736116 0.8820813849400392
   4096 0.8820813914902204 0.8820813893068272
   8192 0.8820813909443684 0.8820813903985197
  16384 0.8820813908079066 0.8820813906714446
  32768 0.8820813907737911 0.8820813907396778
  65536 0.8820813907652575 0.8820813907567422
 131072 0.8820813907631487 0.8820813907610036
 262144 0.8820813907625702 0.8820813907620528
 524288 0.8820813907624605 0.8820813907623183
1048576 0.8820813907624268 0.8820813907623890

Integral over a circle

The test above involves a trivial function \( f(x,y)=1 \). We should also test a non-constant \( f \) function and a more complicated domain. Let \( \Omega \) be a circle at the origin with radius 2, and let \( f=\sqrt{x^2 + y^2} \). This choice makes it possible to compute an exact result: in polar coordinates, \( \int_\Omega f(x,y)dxdy \) simplifies to \( 2\pi\int_0^2 r^2dr = 16\pi/3 \). We must be prepared for quite crude approximations that fluctuate around this exact result. As in the test case above, we experience better results with larger number of points. When we have such evidence for a working implementation, we can turn the test into a proper test function. Here is an example:

    n        midpoint          trapezoidal
      2 0.8842000076332692 0.8770372606158094
      4 0.8827889485397279 0.8806186341245393
      8 0.8822686991994210 0.8817037913321336
     16 0.8821288703366458 0.8819862452657772
     32 0.8820933014203766 0.8820575578012112
     64 0.8820843709743319 0.8820754296107942
    128 0.8820821359746071 0.8820799002925637
    256 0.8820815770754198 0.8820810181335849
    512 0.8820814373412922 0.8820812976045025
   1024 0.8820814024071774 0.8820813674728968
   2048 0.8820813936736116 0.8820813849400392
   4096 0.8820813914902204 0.8820813893068272
   8192 0.8820813909443684 0.8820813903985197
  16384 0.8820813908079066 0.8820813906714446
  32768 0.8820813907737911 0.8820813907396778
  65536 0.8820813907652575 0.8820813907567422
 131072 0.8820813907631487 0.8820813907610036
 262144 0.8820813907625702 0.8820813907620528
 524288 0.8820813907624605 0.8820813907623183
1048576 0.8820813907624268 0.8820813907623890