Bai Giang PPT2017 [PDF]

Zitiervorschau

Phương pháp tính Đào Huy Cường Năm 2017

Mục lục 1 Phép tính đúng 1.1 Các dạng biểu diễn của đa thức 1.2 Đa thức Taylor . . . . . . . . . 1.3 Đa thức nội suy Lagrange . . . 1.4 Đa thức nội suy Newton . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

3 3 7 8 9

2 Phép tính gần đúng 12 2.1 Sai số . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2 Tính gần đúng giá trị hàm sơ cấp cơ bản . . . . . . . . . . . . 14 3 Phương pháp số 3.1 Phương pháp 3.2 Phương pháp 3.3 Phương pháp 3.4 Phương pháp 3.5 Phương pháp

giải phương trình chia đôi . . . . . . . lặp điểm bất động . lặp Newton . . . . dây cung . . . . . . Regula Falsi . . . .

một . . . . . . . . . . . . . . .

ẩn . . . . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

4 Tính gần đúng đạo hàm và tích phân 4.1 Tính gần đúng đạo hàm cấp một . . . . . . . . . . . . . 4.1.1 Công thức hai nút . . . . . . . . . . . . . . . . . 4.1.2 Công thức ba nút . . . . . . . . . . . . . . . . . . 4.1.3 Tính không ổn định của sai số khi tính xấp xỉ đạo 4.2 Tính gần đúng đạo hàm cấp hai . . . . . . . . . . . . . . 4.3 Tính gần đúng tích phân xác định . . . . . . . . . . . . . 4.3.1 Công thức Hình thang . . . . . . . . . . . . . . . 4.3.2 Công thức Trung điểm . . . . . . . . . . . . . . . 4.3.3 Công thức Simpson . . . . . . . . . . . . . . . . . 4.3.4 Công thức Hình thang kết hợp . . . . . . . . . . . 4.3.5 Công thức Simpson kết hợp . . . . . . . . . . . . 4.3.6 Công thức Trung điểm kết hợp . . . . . . . . . . 4.3.7 Tính ổn định của sai số khi tính xấp xỉ tích phân 1

. . . . .

17 18 19 21 23 23

. . . . . . . . . hàm . . . . . . . . . . . . . . . . . . . . . . . . . . .

25 25 25 26 27 28 29 29 30 30 31 32 32 32

. . . . .

. . . . .

5 Phương pháp lặp giải gần đúng hệ phương trình tuyến tính bậc nhất 34 5.1 Phương pháp lặp Jacobi . . . . . . . . . . . . . . . . . . . . . 35 5.2 Phương pháp lặp Gauss-Seidel . . . . . . . . . . . . . . . . . . 36 6 Phương pháp số giải phương trình vi phân bậc nhất 38 6.1 Phương pháp Euler . . . . . . . . . . . . . . . . . . . . . . . . 38 6.2 Phương pháp Runge-Kutta bậc hai . . . . . . . . . . . . . . . 41

2

Chương 1 Phép tính đúng 1.1

Các dạng biểu diễn của đa thức

Định nghĩa 1.1. Cho x, h ∈ R và n ∈ N. Lũy thừa suy rộng bậc n bước nhảy h của x, ký hiệu x[n;h] , là một số thực được cho bởi công thức     x(x − h) . . . x − (n − 1)h , nếu n > 0, x[n;h] = (1.1) 1 , nếu n = 0. Khi bước nhảy h = 1, ta viết x[n] thay cho x[n;1] . Ví dụ 1.1. Tính giá trị của các lũy thừa suy rộng: 2[3;0] , 3[3;1] , 1[3;2] , (−2)[4;−1] . Nhận xét 1.1. • Nếu bước nhảy h = 0 thì x[n;0] = xn . • Nếu bậc n = 1 thì x[1;h] = x. • Công thức giảm bậc x[n;h] = x[k;h] × (x − kh)[n−k;h] ,

k = 0, 1, . . . , n.

Đặc biệt x[n;h] = x × (x − h)[n−1;h] ,   x[n;h] = x[n−1;h] × x − (n − 1)h .

3

• Lũy thừa suy rộng của số đối (−x)[n;h] = (−1)n x[n;−h] . Định nghĩa 1.2. Cho x0 , h ∈ R, n ∈ N và P (x) là một đa thức có bậc không vượt quá n. Ta nói P (x) được biểu diễn dưới dạng • chính tắc, nếu P (x) được viết dưới dạng P (x) = an xn + . . . + aj xj + . . . + a0 .

(1.2)

• chuẩn tắc, nếu P (x) được viết dưới dạng P (x) = an (x − x0 )n + . . . + aj (x − x0 )j + . . . + a0 .

(1.3)

• chính tắc suy rộng, nếu P (x) được viết dưới dạng P (x) = an x[n] + . . . + aj x[j] + . . . + a0 .

(1.4)

• chuẩn tắc suy rộng, nếu P (x) được viết dưới dạng P (x) = an (x − x0 )[n;h] + . . . + aj (x − x0 )[j;h] + . . . + a0 .

(1.5)

Nhận xét 1.2. • Nếu h = 0 và x0 = 0 thì dạng chuẩn tắc suy rộng (1.5) trở thành dạng chính tắc (1.2). • Nếu h = 0 thì dạng chuẩn tắc suy rộng (1.5) trở thành dạng chuẩn tắc (1.3). • Nếu h = 1 thì dạng chuẩn tắc suy rộng (1.5) trở thành dạng chính tắc suy rộng (1.4). Ví dụ 1.2. Cho đa thức P (x) = 2x[3;2] + 3x[2;2] − 4x[1;2] + 5. Tính giá trị P (1).

4

Định lý 1.1 (Hoocne suy rộng). Cho α, x0 , h ∈ R và đa thức P (x) ở dạng chuẩn tắc suy rộng P (x) = an (x − x0 )[n;h] + . . . + aj (x − x0 )[j;h] + . . . + a0 . Xét dãy số thực {bn , bn−1 , . . . , b0 } được cho bởi   b n = an ,    bj = α − x0 − jh · bj+1 + aj , j = n − 1, . . . , 0. Khi đó P (x) = (x − α) · Q(x) + b0 ,

(1.6)

trong đó Q(x) = bn (x − x0 )[n−1;h] + . . . + bj (x − x0 )[j−1;h] + . . . + b1 .

Nhận xét 1.3. Ta có thể dùng sơ đồ Hoocne suy rộng sau đây để tính các hệ số {bn , bn−1 , . . . , b0 }: an an−1 an−2 an−3 a1 a0 α − x0 − (n − 1)h bn bn−1 α − x0 − (n − 2)h bn−1 bn−2 bn−2 bn−3 α − x0 − (n − 3)h ... ... ... α − x0 − h b2 b1 α − x0 b1 b0 Hơn nữa, nếu P (x) ở dạng chính tắc (x0 = h = 0), thì ta có sơ đồ Hoocne quen thuộc: an an−1 an−2 . . . a1 a0 α bn bn−1 bn−2 . . . b1 b0 Ví dụ 1.3. Cho đa thức P (x) = 4(x − 1)[3;2] − 3(x − 1)[2;2] + (x − 1)[1;2] − 7. Tính P (2), P (−2). 5

Ví dụ 1.4. Cho đa thức P (x) = 3x4 − 4x3 + 5x2 + 2x − 7. Tìm thương và số dư khi chia đa thức P (x) cho (x − 1). Ví dụ 1.5. Cho đa thức Q(x) = 3x3 − x2 + 4x + 6. Tìm đa thức P (x) = (x − 1)Q(x) − 1. Ví dụ 1.6. Cho đa thức P (x) = 3x3 − 5x2 + 6x − 7. Hãy viết lại đa thức P (x) dưới dạng chuẩn tắc P (x) = b3 (x − 1)3 + b2 (x − 1)2 + b1 (x − 1) + b0 . Nhận xét 1.4. Cho x0 , y0 , h, k ∈ R và đa thức P (x) ở dạng chuẩn tắc suy rộng P (x) = an (x − x0 )[n;h] + . . . + aj (x − x0 )[j;h] + . . . + a0 . Xét dãy số {b0 , . . . , bn−1 } là các số dư được sinh ra do thực hiện liên tiếp các phép chia đa thức P (x) = (x − y0 ) · Qn−1 (x) + b0 , Qn−1 (x) = (x − y0 − k) · Qn−2 (x) + b1 , Qn−2 (x) = (x − y0 − 2k) · Qn−3 (x) + b2 , ...   Q2 (x) = x − y0 − (n − 2)k · Q1 (x) + bn−2 ,   Q1 (x) = x − y0 − (n − 1)k · Q0 (x) + bn−1 . Khi đó, đa thức P (x) trở thành P (x) = bn (x − y0 )[n;k] + . . . + bj (x − y0 )[j;k] + . . . + b0 , trong đó b n = an . 6

Ví dụ 1.7. Cho đa thức P (x) = 4x3 − 2x2 + 4x − 7. Hãy viết lại đa thức P (x) dưới dạng chuẩn tắc suy rộng P (x) = b3 (x − 1)[3;2] + b2 (x − 1)[2;2] + b1 (x − 1)[1;2] + b0 . Ví dụ 1.8. Cho đa thức P (x) = 4(x − 1)[3;−2] − 2(x − 1)[2;−2] + 4(x − 1)[1;−2] − 7. Hãy viết lại đa thức P (x) dưới dạng chính tắc.

1.2

Đa thức Taylor

Định nghĩa 1.3. Giả sử f ∈ C n [a, b] và x0 ∈ (a, b). Ta gọi Pn (x) = f (x0 ) + f 0 (x0 )(x − x0 ) +

f 00 (x0 ) f (n) (x0 ) (x − x0 )2 + . . . + (x − x0 )n 2! n!

là đa thức Taylor (Taylor polynomial) cấp n của hàm số f tại điểm x0 . Nhận xét 1.5. • Đa thức Taylor của hàm f tại điểm x0 = 0 được gọi là đa thức Maclaurin (Maclaurin polynomial). • Đa thức Taylor cấp n của một đa thức bất kỳ có bậc không quá n tại điểm x0 tùy ý là chính đa thức đó. Ví dụ 1.9. Xác định đa thức Maclaurin cấp 3 của hàm số sau: f (x) = ln(1 + x). Ví dụ 1.10. Xác định đa thức Taylor cấp 4 của hàm số sau tại điểm x0 = 1: f (x) = x4 − 3x3 + 4x2 − 5x + 6. 7

1.3

Đa thức nội suy Lagrange

Định lý 1.2. Giả sử x0 , x1 , . . . , xn là (n+1) số thực khác nhau và y0 , y1 , . . . , yn là các giá trị thực tùy ý. Khi đó, tồn tại duy nhất một đa thức P (x) có bậc không vượt quá n thỏa mãn P (xj ) = yj ,

j = 0, 1, . . . , n.

Định nghĩa 1.4. Cho hàm số f : R → R và x0 , x1 , . . . , xn là (n + 1) số thực khác nhau. Đặt Lk (x) =

n Y

x − xj , x k − xj j=0,j6=k

k = 0, 1, . . . , n.

Đa thức P (x) = f (x0 )L0 (x) + f (x1 )L1 (x) + . . . + f (xn )Ln (x)

(1.7)

được gọi là đa thức nội suy Lagrange (Lagrange interpolating polynomial) của f (x) tại các điểm x0 , x1 , . . . , xn . • Các đa thức Lk (x) được gọi là các đa thức Lagrange cơ sở. • Các điểm x0 , x1 , . . . , xn được gọi là các mốc nội suy (nodes). Nhận xét 1.6. • Giá trị của đa thức Lagrange cơ sở tại các mốc nội suy là ( 0, nếu k 6= j, Lk (xj ) = với mỗi k, j ∈ {0, 1, . . . , n}. 1, nếu k = j, • Giá trị của đa thức nội suy Lagrange tại các mốc nội suy là P (xj ) = f (xj ),

j = 0, 1, . . . , n.

• Đa thức nội suy Lagrange của một đa thức bất kỳ có bậc không vượt quá n tại các mốc nội suy khác nhau x0 , x1 , . . . , xn là chính đa thức đó. Ví dụ 1.11. Cho hàm số f (x) có giá trị được cho trong bảng sau 1 3 4 −3 x f (x) 5 3 −1 2 Xác định đa thức nội suy Lagrange của f (x). 8

1.4

Đa thức nội suy Newton

Giả sử P (x) là đa thức nội suy Lagrange của hàm f tại (n + 1) mốc nội suy phân biệt x0 , x1 , . . . , xn . Ta biết rằng đa thức P (x) là duy nhất và có dạng biểu diễn được cho trong (1.7). Tuy nhiên, ta muốn biểu diễn lại đa thức P (x) dưới dạng P (x) = a0 +a1 (x−x0 )+a2 (x−x0 )(x−x1 )+. . .+an (x−x0 ) . . . (x−xn−1 ), (1.8) trong đó a0 , a1 , . . . , an là các hằng số. Trước hết, để xác định a0 , ta thay x = x0 vào (1.8) a0 = P (x0 ) = f (x0 ). Tương tự, để xác định a1 , ta thay x = x1 vào (1.8) a1 =

f (x1 ) − f (x0 ) . x1 − x 0

Để thuận tiện cho việc trình bày cách xác định các hệ số ak , ta cần đưa ra ký hiệu tỷ sai phân (divided difference) như sau. Tỷ sai phân cấp 0 của hàm f ứng với xi , ký hiệu f [xi ], đơn giản là giá trị của f tại xi : f [xi ] = f (xi ). Tiếp theo, tỷ sai phân cấp 1 ứng với xi , xi+1 , ký hiệu f [xi , xi+1 ], được định nghĩa bởi f [xi+1 ] − f [xi ] . f [xi , xi+1 ] = xi+1 − xi Tỷ sai phân cấp 2 ứng với xi , xi+1 , xi+2 , ký hiệu f [xi , xi+1 , xi+2 ], được định nghĩa bởi f [xi+1 , xi+2 ] − f [xi , xi+1 ] . f [xi , xi+1 , xi+2 ] = xi+2 − xi Tương tự, sau khi các tỷ sai phân cấp (k − 1) f [xi , . . . , xi+k−1 ],

f [xi+1 , . . . , xi+k ]

được xác định, tỷ sai phân cấp k ứng với xi , . . ., xi+k , ký hiệu f [xi , . . . , xi+k ], được định nghĩa bởi f [xi , . . . , xi+k ] =

f [xi+1 , . . . , xi+k ] − f [xi , . . . , xi+k−1 ] . xi+k − xi

9

Quá trình được kết thúc với tỷ sai phân cấp n f [x0 , x1 , . . . , xn ] =

f [x1 , . . . , xn ] − f [x0 , . . . , xn−1 ] . xn − x0

Sử dụng các ký hiệu tỷ sai phân từ cấp 0 đến cấp n được định nghĩa như trên, ta xác định được các hệ số ak như sau: ak = f [x0 , . . . , xk ],

k = 0, 1, . . . , n.

Định nghĩa 1.5. Cho hàm số f : R → R và x0 , x1 , . . . , xn là (n + 1) số thực khác nhau. Với mỗi k ∈ {1, . . . , n}, đặt Nk (x) = (x − x0 )(x − x1 ) . . . (x − xk−1 ). Đa thức P (x) = f (x0 ) + f [x0 , x1 ]N1 (x) + f [x0 , x1 , x2 ]N2 (x) + . . . + f [x0 , . . . , xn ]Nn (x) (1.9) được gọi là đa thức nội suy Newton của f (x) tại các điểm x0 , x1 , . . . , xn . • Các đa thức Nk (x) được gọi là các đa thức Newton cơ sở. Nhận xét 1.7. • Giá trị của đa thức Newton cơ sở tại một số mốc nội suy là Nk (xj ) = 0,

j = 0, . . . , k − 1.

• Giá trị của đa thức nội suy Newton tại các mốc nội suy là P (xj ) = f (xj ),

j = 0, 1, . . . , n.

• Đa thức nội suy Newton của một đa thức bất kỳ có bậc không vượt quá n tại các mốc nội suy khác nhau x0 , x1 , . . . , xn là chính đa thức đó. Ví dụ 1.12. Cho hàm số f (x) có giá trị được cho trong bảng sau x 1 3 4 −3 f (x) 5 3 −1 2 Xác định đa thức nội suy Newton của f (x). 10

Nhận xét 1.8. • Nếu x0 , x1 , . . . , xn là các mốc nội suy cách đều, tức là xj = x0 + jh,

h 6= 0,

j = 1, 2, . . . , n,

thì các tỷ sai phân f [x0 , . . . , xk ] được biểu diễn thông qua các sai phân (difference) ∆k f như sau: f [x0 , . . . , xk ] =

1 ∆k f (x0 ), k!hk

k = 1, . . . , n,

trong đó ∆f (x) = f (x + h) − f (x), ∆k f = ∆(∆k−1 f ),

k = 2, . . . , n.

• Nếu xj = x0 + jh,

h 6= 0,

j = 1, 2, . . . , n,

thì đa thức nội suy Newton (1.9) có dạng biểu diễn chuẩn tắc suy rộng như sau: P (x) = f (x0 )+

2 ∆n f (x0 ) ∆1 f (x0 ) [1;h] ∆ f (x0 ) [2;h] (x−x ) + (x−x ) +. . .+ (x−x0 )[n;h] . 0 0 1!h1 2!h2 n!hn (1.10)

Ví dụ 1.13. Cho hàm số f (x) có giá trị được cho trong bảng sau 1 3 5 7 x f (x) 5 3 −1 2 Xác định đa thức nội suy Newton của f (x). Định lý 1.3. Giả sử f ∈ C n [a, b] và x0 , x1 , . . . , xn là (n + 1) điểm phân biệt thuộc đoạn [a, b]. Khi đó, tồn tại ξ ∈ (a, b) sao cho f [x0 , . . . , xn ] =

f (n) (ξ) . n!

Chứng minh. Xem [1], Theorem 3.6, trang 128.

11

Chương 2 Phép tính gần đúng 2.1

Sai số

Định nghĩa 2.1. Giả sử a là số gần đúng của số A. Ta gọi • |a − A| là sai số tuyệt đối (absolute error) của số gần đúng a, •

|a − A| là sai số tương đối (relative error) của số gần đúng a. |A|

Nhận xét 2.1. Thông thường, ta không thể biết chính xác sai số tuyệt đối |a − A| mà chỉ có thể ước lượng một chặn trên ∆a của nó, tức là |a − A| ≤ ∆a .

(2.1)

Khi đó, ta viết A = a ± ∆a . Ví dụ 2.1. Nếu xấp xỉ π bởi số gần đúng 3.14, hãy ước lượng sai số tuyệt đối. Định nghĩa 2.2. Cho số thực A   A = ± dm 10m + dm−1 10m−1 + . . . + di 10i + . . . , trong đó m là một số nguyên, dm ∈ {1, 2, . . . , 9} và di ∈ {0, 1, . . . , 9}, 12

i ∈ Z,

i < m.

Cho trước i ∈ Z và i ≤ m. Đặt 

m

i



a = sign(A) · dm 10 + . . . + di 10 , và

( a∗ =

a,

nếu di−1 < 5, i

a + sign(A) · 10 ,

nếu di−1 ≥ 5.

Ta gọi • a là giá trị chặt cụt (chopping) của A đến chữ số hàng i. • a∗ là giá trị làm tròn (rounding) của A đến chữ số hàng i. Nhận xét 2.2. Gọi a, a∗ lần lượt là giá trị chặt cụt và giá trị làm tròn của A đến chữ số hàng i. Khi đó 1 |A − a| ≤ 10i , |A − a∗ | ≤ 10i . 2 Ví dụ 2.2. Tìm giá trị chặt cụt và giá trị làm tròn của π đến chữ số thứ tư sau dấu chấm. Mệnh đề 2.1. Giả sử A = a ± ∆a , Khi đó

B = b ± ∆b .

A + B = (a + b) ± (∆a + ∆b ), A − B = (a − b) ± (∆a + ∆b ), AB = (ab) ± (|b|∆a + |a|∆b + ∆a ∆b ), A a |b|∆a + |a|∆b = ± . B b (|b| − ∆b )|b|

Ví dụ 2.3. Cho A = 1.3456 ± 10−4 và B = 7.4632 ± 10−4 . Tính các giá trị gần đúng A của A + B, A · B, và ước lượng sai số tuyệt đối tương ứng. B

13

2.2

Tính gần đúng giá trị hàm sơ cấp cơ bản

Định lý 2.1 (Taylor’s Theorem). Giả sử f ∈ C n+1 [a, b] và x0 ∈ (a, b). Khi đó, với mỗi x ∈ [a, b], tồn tại ξ(x) nằm giữa x và x0 sao cho f (x) = Pn (x) +

f (n+1) (ξ(x)) (x − x0 )n+1 , (n + 1)!

trong đó Pn (x) là đa thức Taylor cấp n của hàm số f tại điểm x0 . Nhận xét 2.3. • Đa thức Pn (x) được dùng để xấp xỉ f (x) với ước lượng sai số chặt cụt (truncation error) là |f (x) − Pn (x)| ≤

M |(x − x0 )n+1 |, (n + 1)!

trong đó M = max |f (n+1) (x)|. x∈[a,b]

• Ta có khai triển một số hàm số sơ cấp thành chuỗi lũy thừa như sau: x2 xn x + + ... + + . . . , −∞ < x < ∞. 1! 2! n! x3 x5 x2n−1 sin x = x − + + . . . + (−1)n + . . . , −∞ < x < ∞. 3! 5! (2n − 1)! x2 x4 x2n cos x = 1 − + + . . . + (−1)n + . . . , −∞ < x < ∞. 2! 4! (2n)! α(α − 1) 2 α(α − 1) . . . (α − n + 1) n (1 + x)α = 1 + αx + x + ... + x + . . . , −1 < x < 1. 2! n! x x2 x3 x4 xn ln(1 + x) = − + − + . . . + (−1)n−1 + . . . , −1 < x < 1. 1 2 3 4 n 2n−1 x x3 x 5 x7 x arctan x = − + − + . . . + (−1)n−1 + . . . , −1 ≤ x ≤ 1. 1 3 5 7 2n − 1 ex = 1 +

Ví dụ 2.4. Tính giá trị gần đúng của eπ với sai số không vượt quá 10−2 , làm tròn đến chữ số thứ hai sau dấu chấm. Thực hiện theo các bước sau:

14

1. Xấp xỉ ex (với x = π) bởi đa thức Taylor cấp n Pn (x) = 1 +

x2 xn x + + ... + . 1! 2! n!

2. Chọn n để sai số chặt cụt không vượt quá |ex − Pn (x)| =

1 · 10−2 , tức là 6

ξ n+1 1 ≤ · 10−2 , (n + 1)! 6

0 < ξ < x.

3. Làm tròn số đúng x = π thành số gần đúng x∗ đến chữ số thứ l sau dấu chấm. Chọn l để |Pn (x) − Pn (x∗ )| ≤

1 · 10−2 . 6

4. Giả sử Pn (x∗ ) = a0 + a1 + . . . + an . Làm tròn mỗi số hạng ai thành a∗i đến chữ số thứ m sau dấu chấm. Chọn m để n n X 1 X ∗ ai − ai ≤ · 10−2 . 6 i=0 i=0 5. Làm tròn

n P

a∗i đến chữ số thứ hai sau dấu chấm.

i=0

Định lý 2.2. Giả sử x0 , x1 , . . . , xn là (n + 1) số thực khác nhau trong đoạn [a, b] và f ∈ C n+1 [a, b]. Gọi P (x) là đa thức nội suy của f tại các mốc nội suy x0 , x1 , . . . , xn . Khi đó, với mỗi x ∈ [a, b], tồn tại ξ(x) ∈ (a, b) sao cho f (x) = P (x) +

f (n+1) (ξ(x)) (x − x0 )(x − x1 ) . . . (x − xn ). (n + 1)!

Chứng minh. Xem [1], Theorem 3.3, trang 112. Nhận xét 2.4. Do f ∈ C n+1 [a, b] nên hàm |f (n+1) | đạt giá trị lớn nhất trên đoạn [a, b]. Do đó ta có ước lượng của sai số chặt cụt |f (x) − P (x)| ≤

M |(x − x0 )(x − x1 ) . . . (x − xn )|, (n + 1)!

trong đó M = max |f (n+1) (x)|. x∈[a,b]

15

Ví dụ 2.5. Cho hàm f (x) = cos x. Hãy thiết lập đa thức nội suy của f (x) tại các mốc nội suy x0 = 0,

x1 =

π , 2

x2 =

Từ đó tính giá trị gần đúng của cos

π , 3

x3 =

π , 4

x4 =

π . 6

π và ước lượng sai số tuyệt đối. 5

16

Chương 3 Phương pháp số giải phương trình một ẩn Định nghĩa 3.1. Cho hàm số f : D → R, với D ⊂ R. Xét phương trình f (x) = 0.

(3.1)

Số thực r ∈ D được gọi là nghiệm của phương trình (3.1) nếu f (r) = 0. Định nghĩa 3.2. Nếu phương trình (3.1) có nghiệm duy nhất trên khoảng (a, b) thì ta gọi (a, b) là một khoảng phân ly nghiệm của phương trình (3.1). Định nghĩa 3.3. Giả sử lim F (h) = L,

lim G(h) = 0.

h→0

h→0

Nếu tồn tại số dương K sao cho |F (h) − L| ≤ K · |G(h)|, với h đủ nhỏ, thì ta viết F (h) − L = O(G(h)). ∞ Định nghĩa 3.4. Giả sử {xn }∞ n=1 là dãy số hội tụ về r và {αn }n=1 là dãy số hội tụ về 0. Nếu tồn tại số dương K sao cho

|xn − r| ≤ K · |αn |, với n đủ lớn, tức là xn − r = O(αn ), thì ta nói dãy số {xn }∞ n=1 hội tụ về r với tốc độ hội tụ O(αn ). 17

Định nghĩa 3.5. Giả sử {xn }∞ n=1 là dãy số hội tụ về r, với xn 6= r với mọi n. Nếu tồn tại các số dương α, λ sao cho |xn+1 − r| = λ, n→∞ |xn − r|α lim

thì ta nói bậc hội tụ của dãy số {xn }∞ n=1 là α. Ví dụ 3.1. Xác định tốc độ hội tụ và bậc hội tụ của các dãy số sau: xn = 1 +

3.1

n+2 , n2

yn = 1 +

n+5 , n3

zn = 1 +

1 . 22n

Phương pháp chia đôi

Xét phương trình (3.1), với f ∈ C[a, b] thỏa mãn f (a) · f (b) < 0 và (a, b) là khoảng phân ly nghiệm. Ý tưởng của phương pháp chia đôi (Bisection method) là thu nhỏ khoảng phân ly đến một mức độ cần thiết. Đầu tiên, đặt a1 = a, b1 = b và tính trung điểm x1 của đoạn [a1 , b1 ], tức là a1 + b 1 . x1 = 2 Nếu f (x1 ) = 0 thì r = x1 , và ta tìm được nghiệm. Ngược lại, khi đó f (x1 ) hoặc cùng dấu với f (a1 ) hoặc cùng dấu với f (b1 ). • Nếu f (a1 ) · f (x1 ) < 0 thì r ∈ (a1 , x1 ). Đặt a2 = a1 , b2 = x1 . • Nếu f (x1 ) · f (b1 ) < 0 thì r ∈ (x1 , b1 ). Đặt a2 = x1 , b2 = b1 . Khi đó, ta lặp lại quá trình trên với khoảng phân ly nghiệm là (a2 , b2 ). Định lý 3.1. Giả sử f ∈ C[a, b] và (a, b) là khoảng phân ly nghiệm của f . Phương pháp chia đôi sinh ra dãy {xn }∞ n=1 hội tụ về nghiệm duy nhất r của f với b−a (3.2) |xn − r| ≤ n , n ∈ N∗ . 2 Chứng minh. Xem [1], Theorem 2.1, trang 51. Nhận xét 3.1. • Mỗi số hạng của dãy {xn }∞ n=1 đều được xem là nghiệm gần đúng của phương trình (3.1) trên khoảng phân ly nghiệm (a, b) với ước lượng sai số được cho bởi (3.2). 18

• Tốc độ hội tụ của phương pháp chia đôi là O(1/2n ), bậc hội tụ là bậc nhất. • Ưu điểm của phương pháp chia đôi là đơn giản, dễ code. Nhược điểm của phương pháp này là tốc độ hội tụ chậm. Ví dụ 3.2. Sử dụng phương pháp chia đôi để tìm nghiệm gần đúng trên đoạn [0, 1] của phương trình cos x − x = 0. với sai số không vượt quá 10−3 , làm tròn đến chữ số thứ ba sau dấu chấm.

3.2

Phương pháp lặp điểm bất động

Định nghĩa 3.6. Điểm r được gọi là điểm bất động (fixed point) của hàm số g cho trước nếu g(r) = r. Định lý 3.2. Giả sử g ∈ C[a, b] và g(x) ∈ [a, b] với mọi x ∈ [a, b]. Khi đó, g có ít nhất một điểm bất động trong [a, b]. Hơn nữa, nếu g có đạo hàm trên khoảng (a, b) và tồn tại một hằng số k ∈ (0, 1) sao cho |g 0 (x)| ≤ k,

∀x ∈ (a, b),

thì g chỉ có duy nhất một điểm bất động trong [a, b]. Chứng minh. Xem [1], Theorem 2.3, trang 57. Ví dụ 3.3. Chứng minh hàm số sau có một điểm bất động duy nhất trên đoạn [0, 1]: g(x) = cos x. Để tìm giá trị gần đúng cho điểm bất động r của hàm g, ta chọn một giá trị xấp xỉ ban đầu x0 và tạo ra dãy {xn }∞ n=0 như sau: xn+1 = g(xn ),

n ∈ N.

(3.3)

Khi đó, nếu dãy số này hội tụ đến điểm r và g là hàm liên tục, thì   r = lim xn = lim g(xn−1 ) = g lim xn−1 = g(r), n→∞

n→∞

n→∞

và ta đạt được điểm bất động r. Kỹ thuật này được gọi là phương pháp lặp điểm bất động (fixed-point iteration). 19

Định lý 3.3. Cho g ∈ C[a, b] thỏa mãn g(x) ∈ [a, b] với mọi x ∈ [a, b]. Hơn nữa, giả sử g có đạo hàm trên khoảng (a, b) và tồn tại một hằng số k ∈ (0, 1) sao cho |g 0 (x)| ≤ k, ∀x ∈ (a, b). Khi đó, cho trước x0 ∈ [a, b] tùy ý, dãy lặp {xn }∞ n=0 được cho bởi (3.3) luôn hội tụ về điểm bất động duy nhất r của g. Hơn nữa, ta có các đánh giá sai số: |xn − r| ≤ k n · max{x0 − a, b − x0 }, n ∈ N∗ , (3.4) hay |xn − r| ≤

k · |xn − xn−1 |, 1−k

hay |xn − r| ≤

kn · |x1 − x0 |, 1−k

n ∈ N∗ , n ∈ N∗ .

(3.5)

(3.6)

hay Chứng minh. Xem [1], Theorem 2.3, Corollary 2.5, trang 62. Nhận xét 3.2. • Đánh giá (3.4) và (3.6) được dùng để tìm số bước lặp n để sai số tuyệt đối của giá trị gần đúng xn không vượt quá một ước lượng sai số cho trước. • Đánh giá (3.5) được dùng khi cần chỉ ra nghiệm gần đúng xn cùng với ước lượng sai số tương ứng với nó. • Nếu có thêm giả thiết g ∈ C 1 [a, b] và g 0 (r) 6= 0 thì phương pháp lặp điểm bất động có bậc hội tụ là bậc nhất. • Nếu có thêm giả thiết g ∈ C 2 [a, b] và g 0 (r) = 0 thì phương pháp lặp điểm bất động có bậc hội tụ ít nhất là bậc hai. Ví dụ 3.4. Sử dụng phương pháp lặp điểm bất động để tìm nghiệm gần đúng trên đoạn [0, 1] của phương trình cos x − x = 0. với sai số không vượt quá 10−3 , làm tròn đến chữ số thứ ba sau dấu chấm.

20

3.3

Phương pháp lặp Newton

Giả sử phương trình (3.1) có nghiệm là r, với f ∈ C 2 [a, b]. Lấy x0 ∈ [a, b] là một giá trị gần đúng của r sao cho f 0 (x0 ) 6= 0. Ta khai triển Taylor hàm f (x) đến cấp một tại điểm x0 : f (x) = f (x0 ) + f 0 (x0 ) · (x − x0 ) +

f 00 (ξ(x)) (x − x0 )2 , 2

trong đó ξ(x) là nằm giữa x0 và x. Vì f (r) = 0 nên khi thay x = r vào khai triển trên, ta được 0 = f (x0 ) + f 0 (x0 ) · (r − x0 ) +

f 00 (ξ) (r − x0 )2 , 2

trong đó ξ là hằng số nằm giữa x0 và ξ. Giả sử |x0 − r| rất nhỏ, khi đó phần f 00 (ξ) dư (r − x0 )2 cũng rất nhỏ, do đó 2 0 ≈ f (x0 ) + f 0 (x0 ) · (r − x0 ). Giải tìm r, ta được r ≈ x0 −

f (x0 ) . f 0 (x0 )

f (x0 ) là nghiệm xấp xỉ tiếp theo của phương trình (3.1). f 0 (x0 ) Tiếp tục thực hiện quy trình trên bằng cách thay điểm x0 bởi x1 , ta được nghiệm xấp xỉ tiếp theo là Như vậy x1 = x0 −

x2 = x1 −

f (x1 ) . f 0 (x1 )

Cứ lặp lại quy trình như trên, ta thu được dãy các nghiệm xấp xỉ {xn }∞ n=0 như sau: f (xn ) xn+1 = xn − 0 , n ∈ N. (3.7) f (xn ) Phương pháp (3.7) được gọi là phương pháp lặp Newton (Newton’s method). Được xuất phát bởi một nghiệm xấp xỉ ban đầu x0 , phương pháp lặp Newton cho ra nghiệm xấp xỉ x1 tiếp theo là hoành độ giao điểm của tiếp tuyến tại điểm (x0 , f (x0 )) và trục hoành. Nghiệm xấp xỉ x2 tiếp theo là hoành độ giao điểm của tiếp tuyến tại điểm (x1 , f (x1 )) và trục hoành, và cứ tiếp tục như vậy. 21

Định lý 3.4. Cho f ∈ C 2 [a, b]. Giả sử r ∈ (a, b) là nghiệm của phương trình (3.1) và f 0 (r) 6= 0. Khi đó, tồn tại δ > 0 sao cho dãy nghiệm xấp xỉ {xn }∞ n=0 sinh bởi phương pháp Newton (3.7) luôn hội tụ về r với mọi x0 ∈ [r −δ, r +δ]. Chứng minh. Xem [1], Theorem 2.6, trang 70. Nhận xét 3.3. • Nếu mọi số hạng của dãy lặp (3.7) đều nằm trong đoạn [a, b] và min |f 0 (x)| = m > 0, x∈[a,b]

thì ta có ước lượng sai số |xn − r| ≤

f (xn ) , m

n ∈ N∗ .

(3.8)

• Nếu mọi số hạng của dãy lặp (3.7) đều nằm trong đoạn [a, b] và min |f 0 (x)| = m > 0, x∈[a,b]

max |f 00 (x)| = M, x∈[a,b]

thì ta có ước lượng sai số |xn − r| ≤

M · (xn − xn−1 )2 , 2m

n ∈ N∗ .

(3.9)

• Với các giả thiết như trong Định lý trên, dãy lặp {xn }∞ n=0 có bậc hội tụ ít nhất là bậc hai nếu nghiệm xấp xỉ ban đầu x0 được chọn đủ gần nghiệm chính xác r. Ví dụ 3.5. Sử dụng phương pháp lặp Newton để tìm nghiệm gần đúng của phương trình x2 − 2 = 0, với sai số không vượt quá 10−3 , làm tròn đến chữ số thứ ba sau dấu chấm.

22

3.4

Phương pháp dây cung

Trong phương pháp lặp Newton (3.7) ở trên, ta thay các đạo hàm f 0 (xn ) bằng các xấp xỉ f 0 (xn ) ≈

f (xn ) − f (xn−1 ) , xn − xn−1

n ∈ N∗ .

Khi đó, dãy lặp (3.7) trở thành xn+1 = xn −

f (xn )(xn − xn−1 ) , f (xn ) − f (xn−1 )

n ∈ N∗ .

(3.10)

Phương pháp (3.10) được gọi là phương pháp dây cung (Secant method). Được bắt đầu bởi hai nghiệm xấp xỉ x0 và x1 cho trước, phương pháp dây cung cho ra nghiệm xấp xỉ x2 tiếp theo là hoành độ giao điểm của đường thẳng đi qua hai điểm (x0 , f (x0 )), (x1 , f (x1 )) và trục hoành. Nghiệm xấp xỉ x3 tiếp theo là hoành độ giao điểm của đường thẳng đi qua hai điểm (x1 , f (x1 )), (x2 , f (x2 )) và trục hoành, và tiếp tục như vậy.

3.5

Phương pháp Regula Falsi

Quan sát phương pháp lặp Newton (3.7) và phương pháp dây cung (3.10), ta thấy rằng nghiệm xấp xỉ xn+1 được sinh ra bởi hai phương pháp này có thể không nằm giữa hai nghiệm xấp xỉ trước đó là xn và xn−1 . Tính chất này trái ngược với phương pháp chia đôi. Phương pháp Regula Falsi (còn gọi là False Position) tương tự như phương pháp dây cung (3.10) nhưng bao gồm thêm các thủ tục để chắc chắn rằng nghiệm xấp xỉ xn+1 được sinh ra bởi phương pháp này phải nằm giữa hai nghiệm xấp xỉ trước đó là xn và xn−1 . Giả sử x0 , x1 là hai nghiệm xấp xỉ ban đầu cho trước sao cho f (x0 )·f (x1 ) < 0. Khi đó, nghiệm xấp xỉ tiếp theo x2 được tính giống như phương pháp dây cung (3.10), đó là hoành độ giao điểm của dây cung (x0 , f (x0 )), (x1 , f (x1 )) và trục hoành. Tiếp theo, để đi đến quyết định chọn dây cung nào để tính nghiệm xấp xỉ x3 , ta phải xét dấu tích f (x2 ) · f (x1 ). • Nếu f (x2 ) · f (x1 ) < 0 thì r ∈ (x2 , x1 ). Do đó ta chọn nghiệm xấp xỉ tiếp theo x3 là hoành độ giao điểm của dây cung (x2 , f (x2 )), (x1 , f (x1 )) và trục hoành. • Ngược lại, ta chọn nghiệm xấp xỉ tiếp theo x3 là hoành độ giao điểm của dây cung (x0 , f (x0 )), (x2 , f (x2 )) và trục hoành. 23

Tương tự, khi đã có x3 , ta tiếp tục xét dấu tích f (x3 ) · f (x2 ) để xác định xem ta phải tính x4 dựa trên x2 , x3 hay dựa trên x3 , x1 . Và cứ tiếp tục như vậy.

24

Chương 4 Tính gần đúng đạo hàm và tích phân 4.1 4.1.1

Tính gần đúng đạo hàm cấp một Công thức hai nút

Giả sử f ∈ C 2 [a, b] và x0 ∈ (a, b). Lấy h 6= 0 đủ nhỏ sao cho x0 + h thuộc (a, b). Để tính gần đúng f 0 (x0 ), trước hết ta xấp xỉ f (x) bởi đa thức nội suy bậc nhất tại hai mốc nội suy x0 , x0 + h: f (2) (ξ(x)) ∆f (x0 ) (x − x0 ) + (x − x0 )(x − x0 − h), f (x) = f (x0 ) + h 2!

(4.1)

trong đó ξ(x) nằm giữa x0 và x0 + h. Lấy đạo hàm (4.1), ta được ∆f (x0 )  f (2) (ξ(x)) 0 f 0 (x) = + (x − x0 )(x − x0 − h) h 2!  f (2) (ξ(x))  + (x − x0 − h) + (x − x0 ) . 2! Thay x bởi x0 , ta được f (x0 + h) − f (x0 ) h 00 − f (ξ). (4.2) h 2 f (x0 + h) − f (x0 ) Như vậy, ta có thể dùng tỷ số để xấp xỉ cho giá trị f 0 (x0 ) h M |h| với ước lượng sai số chặt cụt (truncation error) là , trong đó M = 2 max |f 00 (x)|. Công thức này được gọi là sai phân tiến (forward-difference f 0 (x0 ) =

x∈[a,b]

formula) nếu h > 0 và sai phân lùi (backward-difference formula) nếu h < 0. 25

Nhận xét 4.1. Ta có thể thu được công thức hai nút (4.2) bằng cách dùng khai triển Taylor tại điểm x0 đến cấp một. Ví dụ 4.1. Tính gần đúng đạo hàm của hàm số f (x) = ln x tại điểm x0 = 1.8 bằng công thức sai phân tiến với h = 0.1 và ước lượng sai số chặt cụt.

4.1.2

Công thức ba nút

Giả sử f ∈ C 3 [a, b] và x0 ∈ (a, b). Lấy h 6= 0 đủ nhỏ sao cho x0 + h và x0 + 2h đều thuộc (a, b). Xấp xỉ f (x) bởi đa thức nội suy bậc hai tại các mốc nội suy cách đều x0 , x0 + h, x0 + 2h: ∆f (x0 ) ∆2 f (x0 ) (x − x0 ) + (x − x0 )(x − x0 − h) h 2!h2 (4.3) f (3) (ξ(x)) + (x − x0 )(x − x0 − h)(x − x0 − 2h), 3! trong đó ξ(x) nằm giữa x0 và x0 + 2h. Lấy đạo hàm (4.3), ta được  ∆f (x0 ) ∆2 f (x0 )  0 f (x) = + (x − x0 − h) + (x − x0 ) h 2!h2  f (3) (ξ(x)) 0 + (x − x0 )(x − x0 − h)(x − x0 − 2h) 3!  f (3) (ξ(x))  + (x − x0 − h)(x − x0 − 2h) + (x − x0 )(x − x0 − 2h) + (x − x0 )(x − x0 − h) . 3! Thay lần lượt x bởi x0 , x0 + h, x0 + 2h, ta được  h2 1 1 3 − f (x0 ) + 2f (x0 + h) − f (x0 + 2h) + f (3) (ξ0 ), f 0 (x0 ) = h 2 2 3  1  h2 1 1 f 0 (x0 + h) = − f (x0 − h) + f (x0 + h) − f (3) (ξ1 ), h 2 2 6   h2 1 1 3 0 f (x0 + 2h) = f (x0 ) − 2f (x0 + h) + f (x0 + 2h) + f (3) (ξ2 ). h 2 2 3 Bằng cách đổi biến thích hợp, ta có ba công thức tính gần đúng f 0 (x0 ) như sau:  h2 1 0 f (x0 ) = − 3f (x0 ) + 4f (x0 + h) − f (x0 + 2h) + f (3) (ξ0 ), 2h 3   h2 1 f 0 (x0 ) = − f (x0 − h) + f (x0 + h) − f (3) (ξ1 ), 2h 6   h2 1 0 f (x0 ) = f (x0 ) − 4f (x0 − h) + 3f (x0 − 2h) + f (3) (ξ2 ). 2h 3 f (x) =f (x0 ) +

26

Để ý rằng công thức thứ ba có thể đạt được từ công thức thứ nhất bằng cách thay h bởi −h. Do đó, ta có hai công thức ba nút tính xấp xỉ f 0 (x0 ), đó là • công thức ba nút điểm cuối (three-point endpoint formula): f 0 (x0 ) =

 h2 1 − 3f (x0 ) + 4f (x0 + h) − f (x0 + 2h) + f (3) (ξ0 ), (4.4) 2h 3

trong đó ξ0 nằm giữa x0 và x0 + 2h, và • công thức ba nút điểm giữa (three-point midpoint formula): f 0 (x0 ) =

 h2 1 f (x0 + h) − f (x0 − h) − f (3) (ξ1 ), 2h 6

(4.5)

trong đó ξ1 nằm giữa x0 − h và x0 + h. Nhận xét 4.2. Ta có thể thu được các công thức ba nút (4.4), (4.5) bằng cách dùng khai triển Taylor tại điểm x0 đến cấp 2. Ví dụ 4.2. Tính gần đúng đạo hàm của hàm số f (x) = xex tại điểm x0 = 2 bằng các công thức ba nút điểm cuối và ba nút điểm giữa với h = 0.1 và ước lượng các sai số chặt cụt tương ứng.

4.1.3

Tính không ổn định của sai số khi tính xấp xỉ đạo hàm

Xét công thức ba nút điểm giữa để xấp xỉ đạo hàm (4.5): f 0 (x0 ) =

 h2 1 f (x0 + h) − f (x0 − h) − f (3) (ξ). 2h 6

Trong quá trình tính toán, ta thay thế f (x0 + h) và f (x0 − h) bằng các giá trị gần đúng f (x0 + h)∗ , f (x0 − h)∗ và mắc phải các sai số (round-off error) e(x0 + h), e(x0 − h), tức là f (x0 + h) = f (x0 + h)∗ + e(x0 + h),

f (x0 − h) = f (x0 − h)∗ + e(x0 − h).

Khi đó, tổng sai số của phép tính gần đúng đạo hàm là   h2 1 1 ∗ ∗ f (x0 +h) −f (x0 −h) = e(x0 +h)−e(x0 −h) − f (3) (ξ). f (x0 )− 2h 2h 6 0

27

Giả sử các sai số e(x0 ± h) đều bị chặn bởi một số dương ε và đạo hàm cấp ba bị chặn bởi số dương M . Khi đó  ε h2 1 0 ∗ ∗ f (x0 + h) − f (x0 − h) ≤ + M. f (x0 ) − 2h h 6 h2 Để làm giảm sai số chặt cụt (truncation error), M , ta chỉ cần giảm h. Tuy 6 ε nhiên, khi h giảm thì sai số tăng lên. Do đó, trong thực hành, ta ít khi lấy h h quá nhỏ. Ví dụ 4.3. Dùng bảng giá trị làm tròn đến chữ số thứ năm sau dấu chấm của hàm f (x) = sin x để tính xấp xỉ f 0 (0.9) ứng với h = 0.1, h = 0.05, h = 0.01, h = 0.005, h = 0.001. Lấy cos 0.9 = 0.62161 và so sánh các sai số.

4.2

Tính gần đúng đạo hàm cấp hai

Giả sử f ∈ C 4 [a, b] và x0 ∈ (a, b). Lấy h 6= 0 đủ nhỏ sao cho x0 + h, x0 − h đều thuộc (a, b). Để tính xấp xỉ f 00 (x0 ) bởi các giá trị của f tại ba nút x0 − h, x0 , x0 + h, ta sử dụng khai triển Taylor của hàm f tại điểm x0 đến cấp ba: 1 1 1 00 f (x0 )h2 + f 000 (x0 )h3 + f (4) (ξ1 )h4 , 2! 3! 4! 1 00 1 000 1 0 2 3 f (x0 − h) = f (x0 ) − f (x0 )h + f (x0 )h − f (x0 )h + f (4) (ξ2 )h4 , 2! 3! 4!

f (x0 + h) = f (x0 ) + f 0 (x0 )h +

trong đó ξ1 nằm giữa x0 và x0 + h, ξ2 nằm giữa x0 và x0 − h. Từ đó, ta thu được công thức tính gần đúng f 00 (x0 ) như sau:  h2 1 f (x0 ) = 2 f (x0 + h) − 2f (x0 ) + f (x0 − h) − f (4) (ξ), h 12 00

(4.6)

trong đó ξ nằm giữa x0 − h và x0 + h. Ví dụ 4.4. Tính gần đúng đạo hàm cấp hai của hàm số f (x) = xex tại điểm x0 = 2 với h = 0.1 và ước lượng sai số chặt cụt.

28

4.3

Tính gần đúng tích phân xác định

Định lý 4.1 (Weighted Mean Value Theorem for Integrals). Giả sử f ∈ C[a, b], g khả tích Riemann trên đoạn [a, b] và g(x) không đổi dấu trên [a, b]. Khi đó, tồn tại hằng số c ∈ (a, b) sao cho Z b Z b g(x)dx. f (x)g(x)dx = f (c) · a

a

Trong trường hợp g(x) ≡ 1, Định lý 4.1 chính là Định lý trung bình tích phân thông thường: Z b 1 f (x)dx. f (c) = b−a a

4.3.1

Công thức Hình thang

R x +h Để tính gần đúng tích phân x00 f (x)dx, với h > 0, trước hết ta xấp xỉ f (x) bởi đa thức nội suy bậc nhất tại hai nút x0 và x0 + h, với giả thiết f thuộc lớp C 2 : ∆f (x0 ) f 00 (ξ(x)) (x − x0 ) + (x − x0 )(x − x0 − h), f (x) = f (x0 ) + h 2 trong đó ξ(x) nằm giữa x0 và x0 + h. Lấy tích phân trên đoạn [x0 , x0 + h], ta được Z x0 +h f (x)dx x0 x0 +h

Z  ∆f (x0 ) 1 x0 +h 00 = f (x0 ) + (x − x0 ) dx + f (ξ(x))(x − x0 )(x − x0 − h)dx h 2 x0 x0  1 Z x0 +h h  = · f (x0 ) + f (x0 + h) + f 00 (ξ(x))(x − x0 )(x − x0 − h)dx. 2 2 x0 (4.7) Vì biểu thức (x − x0 )(x − x0 − h) không đổi dấu trên đoạn [x0 , x0 + h] nên ta áp dụng Định lý giá trị trung bình 4.1 Z x0 +h Z x0 +h 00 00 f (ξ(x))(x − x0 )(x − x0 − h)dx = f (ξ) (x − x0 )(x − x0 − h)dx, Z



x0

x0

trong đó ξ là hằng số nằm giữa x0 và x0 + h. Thay vào (4.7), ta được công thức Hình thang (Trapezoideal rule) như sau: Z x0 +h  h3 h  (4.8) f (x)dx = · f (x0 ) + f (x0 + h) − f 00 (ξ). 2 12 x0 29

Ví dụ 4.5. Tính gần đúng tích phân Z

1

I= 0

x dx 2+x

bằng công thức hình thang (4.8) và ước lượng sai số chặt cụt.

4.3.2

Công thức Trung điểm

R x +h Để tính gần đúng tích phân x00−h f (x)dx, với h > 0, trước hết ta khai triển Taylor f (x) đến cấp một tại điểm x0 , với giả thiết f thuộc lớp C 2 : f (x) = f (x0 ) + f 0 (x0 )(x − x0 ) +

f 00 (ξ(x)) (x − x0 )2 , 2

trong đó ξ(x) nằm giữa x0 −h và x0 +h. Lấy tích phân trên đoạn [x0 −h, x0 +h] và áp dụng Định lý giá trị trung bình 4.1, ta được công thức Trung điểm (Midpoint rule) như sau: Z x0 +h f 00 (ξ) 3 f (x)dx = 2hf (x0 ) + h, (4.9) 3 x0 −h trong đó ξ nằm giữa x0 − h và x0 + h. Ví dụ 4.6. Tính gần đúng tích phân Z

1

I= 0

x dx 2+x

bằng công thức trung điểm (4.9) và ước lượng sai số chặt cụt.

4.3.3

Công thức Simpson

R x +2h Để tính gần đúng tích phân x00 f (x)dx, với h > 0, trước hết ta khai triển Taylor f (x) đến cấp ba tại điểm x0 + h, với giả thiết f thuộc lớp C 4 : f (x) = f (x0 + h) + f 0 (x0 + h)(x − x0 − h) +

f 00 (x0 + h) (x − x0 − h)2 2

f 000 (x0 + h) f (4) (ξ(x)) 3 (x − x0 − h) + (x − x0 − h)4 , + 6 24 30

trong đó ξ(x) nằm giữa x0 và x0 + 2h. Lấy tích phân trên đoạn [x0 , x0 + 2h] và áp dụng Định lý giá trị trung bình 4.1, ta được Z x0 +2h Z x0 +2h f 00 (x0 + h) 1 3 f (x)dx = f (x0 + h) · 2h + · 2h + f (4) (ξ(x))(x − x0 − h)4 dx 6 24 x0 x0 = 2f (x0 + h)h +

f 00 (x0 + h) 3 f (4) (ξ1 ) 5 ·h + h, 3 60

(4.10) trong đó ξ1 là hằng số nằm giữa x0 và x0 + 2h. Mặt khác, theo công thức tính gần đúng đạo hàm cấp hai (4.6), ta có  h2 1 f 00 (x0 + h) = 2 f (x0 + 2h) − 2f (x0 + h) + f (x0 ) − f (4) (ξ2 ), (4.11) h 12 trong đó ξ2 là hằng số nằm giữa x0 và x0 + 2h. Thay (4.11) vào (4.10), ta được Z x0 +2h  h5  5  h 3 f (x)dx = f (x0 )+4f (x0 +h)+f (x0 +2h) − f (4) (ξ2 )− f (4) (ξ1 ) . 3 90 2 2 x0 Ngoài ra, bằng một phương pháp tiếp cận khác khi đánh giá sai số chặt cụt, các hằng số ξ1 , ξ2 trong biểu thức trên có thể được thay thế bởi một hằng số ξ chung trong khoảng (x0 , x0 + 2h). Điều này dẫn đến công thức Simpson (Simpson’s rule) như sau: Z x0 +2h  h5 h f (x0 ) + 4f (x0 + h) + f (x0 + 2h) − f (4) (ξ). (4.12) f (x)dx = 3 90 x0 Ví dụ 4.7. Tính gần đúng tích phân Z I= 0

1

x dx 2+x

bằng công thức Simpson (4.12) và ước lượng sai số chặt cụt.

4.3.4

Công thức Hình thang kết hợp

b−a và xj = a + jh với j = N 0, 1, . . . , N , với N là số nguyên dương. Khi đó, ta có công thức Hình thang kết hợp (Composite Trapezoidal rule) như sau: Z b  b−a h f (x)dx = f (a)+2(f (x1 )+f (x2 )+. . .+f (xN −1 ))+f (b) − h2 f 00 (ξ), 2 12 a (4.13) trong đó ξ là một hằng số thuộc khoảng (a, b). Định lý 4.2. Giả sử f ∈ C 2 [a, b], h =

31

4.3.5

Công thức Simpson kết hợp

b−a và xj = a + jh với j = 2N 0, 1, . . . , 2N , với N là số nguyên dương. Khi đó, ta có công thức Simpson kết hợp (Composite Simpson rule) như sau: Z b  h f (a) + 2(f (x2 ) + . . . + f (x2N −2 )) + 4(f (x1 ) + . . . + f (x2N −1 )) + f (b) f (x)dx = 3 a b − a 4 (4) − h f (ξ), 180 (4.14) trong đó ξ là một hằng số thuộc khoảng (a, b).

Định lý 4.3. Giả sử f ∈ C 4 [a, b], h =

4.3.6

Công thức Trung điểm kết hợp

b−a và xj = a + (j + 1)h với 2N + 2 j = −1, 0, 1, . . . , 2N, 2N + 1, với N là số nguyên dương. Khi đó, ta có công thức Trung điểm kết hợp (Composite Midpoint rule) như sau: Z b   b−a f (x)dx = 2h f (x0 ) + f (x2 ) + . . . + f (x2N ) + h2 f 00 (ξ), (4.15) 6 a Định lý 4.4. Giả sử f ∈ C 2 [a, b], h =

trong đó ξ là một hằng số thuộc khoảng (a, b).

4.3.7

Tính ổn định của sai số khi tính xấp xỉ tích phân

Xét công thức Hình thang kết hợp (4.13) ứng với hàm f và chia đoạn [a, b] thành N đoạn bằng nhau. Giả sử mỗi giá trị f (xi ) được xấp xỉ bởi giá trị gần đúng là f (xi )∗ với sai số là ei , tức là f (xi ) = f (xi )∗ + ei ,

i = 0, 1, . . . , n.

Giả sử các sai số tuyệt đối |ei | bị chặn đều bởi số dương ε. Khi đó, ta có đánh giá sai số tích lũy e(h) trong công thức Hình thang kết hợp như sau: h  |e(h)| = e0 + 2(e1 + . . . + eN −1 ) + eN ≤ N hε = (b − a)ε. 2 Như vậy |e(h)| bị chặn bởi một giá trị độc lập với h và N . Điều này có nghĩa là mặc dù ta tăng số đoạn chia N để tăng độ chính xác nhưng vẫn không làm tăng sai số. Kết quả này cho thấy tính ổn định của các phương pháp tính gần đúng tích phân khi h tiến về không. 32

Ví dụ 4.8. Tính gần đúng tích phân Z I= 0

1

x dx 2+x

bằng công thức Hình thang kết hợp, với sai số không vượt quá 10−2 , làm tròn đến chữ số thứ hai sau dấu chấm.

33

Chương 5 Phương pháp lặp giải gần đúng hệ phương trình tuyến tính bậc nhất Định nghĩa 5.1. Chuẩn vectơ trong không gian Rn là một ánh xạ đi từ Rn vào R, ký hiệu k · k, thỏa mãn các tính chất sau: a) kxk ≥ 0 với mọi x ∈ Rn . b) kxk = 0 khi và chỉ khi x = 0. c) kα · xk = |α|kxk với mọi α ∈ R và x ∈ Rn . d) kx + yk ≤ kxk + kyk với mọi x, y ∈ Rn . Định nghĩa 5.2. Chuẩn l2 và l∞ cho vectơ x = (x1 , . . . , xn ) được định nghĩa bởi q kxk2 = x21 + . . . + x2n , kxk∞ = max |xi |. 1≤i≤n

Ví dụ 5.1. Tính chuẩn l2 và l∞ cho vectơ x = (1, −2, 3). n Định nghĩa 5.3. Dãy {x(k) }∞ k=1 các vectơ trong không gian R được gọi là hội tụ về vectơ x đối với chuẩn k · k nếu, cho trước ε > 0, tồn tại số nguyên N (ε) sao cho kx(k) − xk < ε, với mọi k ≥ N (ε).

Định nghĩa 5.4. Chuẩn ma trận trong tập hợp các ma trận cỡ n × n là một ánh xạ đi từ tập hợp các ma trận cỡ n × n vào R, ký hiệu k · k, thỏa mãn các tính chất sau: a) kAk ≥ 0. b) kAk = 0 khi và chỉ khi A là ma trận không. 34

c) kα · Ak = |α|kAk. d) kA + Bk ≤ kAk + kBk. e) kA · Bk ≤ kAkkBk. Định lý 5.1. Nếu k · k là chuẩn vectơ trong không gian Rn , thì kAk = max kA · xk kxk=1

là chuẩn ma trận. Chứng minh. Xem [1], Theorem 7.9, trang 438. Hệ quả 5.1. Với mọi ma trận A cỡ n × n và vectơ z 6= 0, ta có kAzk ≤ kAkkzk. Định lý 5.2. Nếu A = (aij ) là ma trận cỡ n × n, thì kAk∞ = max

1≤i≤n

n X

|aij |.

j=1

Chứng minh. Xem [1], Theorem 7.11, trang 440. Ví dụ 5.2. Tính chuẩn l∞ của ma trận 

 1 3 −1 A =  2 −1 2  . −3 2 −2

5.1

Phương pháp lặp Jacobi

Xét hệ phương trình Ax = b.

(5.1)

Từ phương trình thứ i của hệ (5.1), ta rút xi theo các ẩn còn lại (với giả thiết aii 6= 0): xi =

n  X j=1,j6=i

−

aij xj  bi + , aii aii

Hệ (5.1) trở thành x = T x + c. 35

i = 1, . . . , n.

Phương pháp lặp Jacobi (Jacobi iterative method) xây dựng dãy {x(k) }∞ k=0 dùng để xấp xỉ nghiệm của hệ (5.1) như sau: x(k+1) = T x(k) + c,

k ∈ N,

(5.2)

trong đó x(0) ∈ Rn được cho trước. Định lý 5.3. Giả sử kT k < 1. Khi đó, cho trước x(0) ∈ Rn tùy ý, dãy lặp {x(k) }∞ k=0 được cho bởi (5.2) luôn hội tụ về một vectơ duy nhất x là nghiệm của hệ (5.1). Hơn nữa, ta có các đánh giá sai số: kx(k) − xk ≤

kT k · kx(k) − x(k−1) k, 1 − kT k

k ∈ N∗ ,

(5.3)

hay kx(k) − xk ≤

kT kk · kx(1) − x(0) k, 1 − kT k

k ∈ N∗ .

(5.4)

Chứng minh. Xem [1], Theorem 7.19, Corollary 7.20, trang 457, 458. Ví dụ 5.3. Sử dụng phương pháp lặp Jacobi để giải gần đúng hệ phương trình sau    10x1 − 2x2 + x3 = 6, − x1 + 11x2 − x3 = 25,   2x − x + 10x = −11, 1 2 3 với điều kiện dừng kx(k) − x(k−1) k∞ ≤ 10−3 .

5.2

Phương pháp lặp Gauss-Seidel

Ở phương pháp lặp Jacobi (5.2) như trên, ta thấy rằng để tính tất cả các tọa độ xi (k+1) của nghiệm gần đúng x(k+1) , ta chỉ sử dụng các tọa độ xi (k) của nghiệm gần đúng x(k) ở bước trước đó. Tuy nhiên, với i > 1, ta có thể tính tọa độ xi (k+1) bằng cách sử dụng các tọa độ x1 (k+1) , . . . , xi−1 (k+1) của nghiệm gần đúng x(k+1) đã được tính và các tọa độ xi+1 (k) , . . . , xn (k) của nghiệm gần đúng x(k) . Cụ thể là xi

(k+1)

i−1 n  X 1 X (k+1) (k) − aij xj − aij xj + bi , = aii j=1 j=i+1

36

i = 2, . . . , n.

(5.5)

Ví dụ 5.4. sau

Sử dụng phương pháp lặp Gauss-Seidel để giải gần đúng hệ phương trình    10x1 − 2x2 + x3 = 6, − x1 + 11x2 − x3 = 25,   2x1 − x2 + 10x3 = −11,

với điều kiện dừng kx(k) − x(k−1) k∞ ≤ 10−3 .

37

Chương 6 Phương pháp số giải phương trình vi phân bậc nhất Xét bài toán giá trị ban đầu cho phương trình vi phân ( y 0 (t) = f (t, y(t)), a ≤ t ≤ b, y(a) = α.

(6.1)

Định lý 6.1. Cho hàm hai biến f (t, y) liên tục và thỏa mãn điều kiện Lipschitz theo biến y trên miền D = [a, b] × R. Khi đó, bài toán giá trị ban đầu (6.1) có duy nhất một nghiệm y(t) với a ≤ t ≤ b. Giả sử y(t) là nghiệm chính xác duy nhất của bài toán giá trị ban đầu (6.1). Thông thường, nghiệm chính xác y(t) sẽ không dễ tìm được; thay vào đó, giá trị xấp xỉ của y(t) tại một số hữu hạn điểm trên đoạn [a, b], mà ta gọi là các điểm lưới (mesh points), có thể tìm được bằng các phương pháp số. Sau khi tính được các giá trị gần đúng của y(t) tại các nút lưới, ta sẽ đạt được nghiệm xấp xỉ của bài toán (6.1) bằng nội suy.

6.1

Phương pháp Euler

Ta chia đoạn [a, b] thành N đoạn bằng nhau bởi các điểm lưới ti = a + ih,

i = 0, 1, . . . , N,

trong đó h=

b−a . N

Ta gọi h là cỡ bước (step size). 38

Giả sử nghiệm chính xác y(t) thuộc lớp C 2 . Khi đó, với mỗi i = 0, 1, . . . , N − 1, ta xấp xỉ hàm y(t) bởi đa thức Taylor cấp một của nó tại điểm ti và tính giá trị của đa thức Taylor đó tại t = ti+1 : y(ti+1 ) = y(ti ) + h · y 0 (ti ) +

h2 00 · y (ξi ), 2

trong đó ξi nằm giữa ti và ti+1 . Vì y 0 (t) = f (t, y(t)) nên ta có y(ti+1 ) = y(ti ) + h · f (ti , y(ti )) +

h2 00 · y (ξi ). 2

Phương pháp Euler tính các giá trị wi ≈ y(ti ) bằng cách cắt bỏ các phần dư h2 00 · y (ξi ) từ đẳng thức trên. Vậy phương pháp Euler như sau: 2 w0 = α, wi+1 = wi + h · f (ti , wi ),

i = 0, 1, . . . , N − 1.

(6.2)

Ta gọi (6.2) là một phương trình sai phân ứng với bài toán (6.1). Ví dụ 6.1. Xét bài toán giá trị ban đầu ( y 0 (t) = y − t2 + 1, y(0) = 0.5.

0 ≤ t ≤ 2,

a) Bằng cách chia lưới đoạn [0, 2] với cỡ bước h = 0.5, hãy dùng phương pháp Euler tính giá trị gần đúng của nghiệm chính xác y(t) tại các điểm lưới. b) Bằng cách chia lưới đoạn [0, 2] với cỡ bước h = 0.2, hãy dùng phương pháp Euler tính giá trị gần đúng của nghiệm chính xác y(t) tại các điểm lưới. Định lý 6.2. Cho hàm f (t, y) liên tục và thỏa mãn điều kiện Lipschitz theo biến y với hằng số L > 0 trên miền D = [a, b] × R, tức là |f (t, y1 ) − f (t, y2 )| ≤ L · |y1 − y2 |,

∀(t, y1 ), (t, y2 ) ∈ D.

Giả sử y(t) là nghiệm duy nhất của bài toán (6.1) và tồn tại hằng số M > 0 sao cho |y 00 (t)| ≤ M, ∀t ∈ [a, b]. Khi đó, ta có các ước lượng sai số cho phương pháp Euler như sau:  hM  L(ti −a) |y(ti ) − wi | ≤ · e − 1 , i = 0, 1, . . . , N. 2L 39

(6.3)

Chứng minh. Xem [1], Theorem 5.9, trang 271. Định nghĩa 6.1 (local truncation error). Giả sử w0 = α, wi+1 = wi + h · φ(ti , wi ),

i = 0, 1, . . . , N − 1.

(6.4)

là một phương trình sai phân ứng với phương trình vi phân (6.1). Khi đó, ta nói phương trình sai phân (6.4) có sai số chặt cụt địa phương (local truncation error) tại bước thứ i + 1 là   y(ti+1 ) − y(ti ) + h · φ(ti , y(ti )) τi+1 (h) = , i = 0, 1, . . . , N − 1. h Định nghĩa 6.2. Phương trình sai phân (6.4) được gọi là tương thích (consitent) với phương trình vi phân (6.1) nếu lim max |τi (h)| = 0.

h→0 0≤i≤N

Định nghĩa 6.3. Phương trình sai phân (6.4) được gọi là hội tụ (convergent) đến phương trình vi phân (6.1) nếu lim max |wi − y(ti )| = 0.

h→0 0≤i≤N

Nhận xét 6.1. Phương trình sai phân của phương pháp Euler (6.2) có sai số chặt cụt địa phương là   y(ti+1 ) − y(ti ) + h · f (ti , y(ti )) τi+1 (h) = h y(ti+1 ) − y(ti ) = − f (ti , y(ti )) h h2 00 0 hy (ti ) + · y (ξi ) 2 = − f (ti , y(ti )) h h = · y 00 (ξi ), i = 0, 1, . . . , N − 1. 2

40

6.2

Phương pháp Runge-Kutta bậc hai

Giả sử nghiệm chính xác duy nhất y(t) của bài toán (6.1) thuộc lớp C 3 . Khi đó, với mỗi i = 0, 1, . . . , N − 1, ta xấp xỉ hàm y(t) bởi đa thức Taylor cấp hai của nó tại điểm ti và tính giá trị của đa thức Taylor đó tại t = ti+1 : y(ti+1 ) = y(ti ) + h · y 0 (ti ) +

h2 00 · y (ti ) + O(h3 ), 2!

i = 0, 1, . . . , N − 1.

Ta có y 0 (t) = f (t, y(t)),

y 00 (t) = f 0 (t, y(t)),

do đó   h 0 y(ti+1 ) = y(ti ) + h f (ti , y(ti )) + f (ti , y(ti )) + O(h3 ). 2

(6.5)

Mặt khác, ta có h f (ti , y(ti )) + f 0 (ti , y(ti )) 2  h  ∂f ∂f 0 = f (ti , y(ti )) + (ti , y(ti )) + (ti , y(ti ))y (ti ) 2 ∂t ∂y ∂f h ∂f h = f (ti , y(ti )) + (ti , y(ti )) + (ti , y(ti )) f (ti , y(ti )) ∂t 2 ∂y 2   h h = f ti + , y(ti ) + f (ti , y(ti )) + O(h2 ). 2 2 Dấu bằng cuối cùng có được do áp dụng khai triển Taylor bậc nhất cho hàm hai biến f (t, y) tại điểm (ti , y(ti )). Thay vào (6.5), ta được   h h y(ti+1 ) = y(ti ) + h · f ti + , y(ti ) + f (ti , y(ti )) + O(h3 ). 2 2 Phương pháp Runge-Kutta bậc hai tính các giá trị wi ≈ y(ti ) bằng cách cắt bỏ phần dư O(h3 ) từ đẳng thức trên. Vậy phương pháp Runge-Kutta bậc hai có phương trình sai phân như sau: w0 = α, wi+1

 h h = wi + h · f ti + , wi + f (ti , wi ) , 2 2 

Nhận xét 6.2.

41

i = 0, 1, . . . , N − 1.

(6.6)

Phương pháp Runge-Kutta bậc hai (6.6) có sai số chặt cụt địa phương là  h h y(ti+1 ) − y(ti ) + h · f ti + , wi + f (ti , wi ) 2 2 = h   y(ti+1 ) − y(ti ) h h = − f ti + , wi + · f (ti , wi ) h 2 2 2 = O(h ), i = 0, 1, . . . , N − 1. 

τi+1



Ví dụ 6.2. Xét bài toán giá trị ban đầu ( y 0 (t) = y − t2 + 1, y(0) = 0.5.

0 ≤ t ≤ 2,

Bằng cách chia lưới đoạn [0, 2] với cỡ bước h = 0.2, hãy dùng phương pháp Runge-Kutta bậc hai tính giá trị gần đúng của nghiệm chính xác y(t) tại các điểm lưới.

42

Bài tập 1. Cho α là một số thực và P (x) là một hàm đa thức. Trình bày cách tính giá trị P (α) khi P (x) được biểu diễn ở một trong bốn dạng: chính tắc, chuẩn tắc, chính tắc suy rộng, chuẩn tắc suy rộng. 2. Trình bày cách chuyển đổi dạng biễu diễn của đa thức từ dạng này sang dạng khác. 3. Cho đa thức P (x) có bậc không vượt quá n và x0 , x1 , . . . , xn là (n + 1) số thực khác nhau. Giả sử đã biết các giá trị P (x0 ), P (x1 ), . . ., P (xn ). Xác định dạng biểu diễn chính tắc của đa thức P (x). 4. Tính gần đúng giá trị của các hàm số sơ cấp với sai số không vượt quá 10−d cho trước, làm tròn đến d chữ số sau dấu chấm. 5. Trình bày các phương pháp tìm nghiệm gần đúng của phương trình một ẩn: phương pháp chia đôi đoạn chứa nghiệm, phương pháp lặp điểm bất động, phương pháp lặp Newton, phương pháp dây cung, phương pháp Regula Falsi. Đánh giá bậc hội tụ của từng phương pháp. 6. Thiết lập các công thức tính gần đúng đạo hàm và ước lượng sai số chặt cụt ứng với từng công thức. 7. Thiết lập các công thức Hình thang, Trung điểm, Simpson để tính gần đúng giá trị của tích phân xác định và ước lượng sai số chặt cụt ứng với từng công thức. 8. Trình bày các phương pháp lặp tìm nghiệm gần đúng của hệ phương trình bậc nhất: phương pháp lặp Jacobi, phương pháp lặp Gauss-Seidel. 9. Trình bày các phương pháp tìm nghiệm xấp xỉ của bài toán giá trị ban đầu cho phương trình vi phân thường bậc nhất: phương pháp Euler, phương pháp Runge-Kutta bậc hai. Ước lượng sai số chặt cụt địa phương của từng phương pháp.

43

Yêu cầu chung: Mỗi vấn đề khảo sát của mỗi nhóm được viết thành một bài thu hoạch và trình bày theo dàn ý sau: 1. Đặt vấn đề. 2. Nội dung vấn đề: • Cơ sở Toán học của vấn đề khảo sát: định nghĩa, định lý, chứng minh định lý. • Đưa ra phương pháp giải quyết vấn đề. • Đưa ra thuật toán tương ứng với phương pháp. • Ví dụ minh họa cho việc vận hành thuật toán. • Đoạn chương trình Matlab/Octave tương ứng. 3. Nhận xét thêm (nếu có) về vấn đề được khảo sát. 4. Tài liệu tham khảo. Lưu ý: Các kết quả bằng số trong thuật toán và ví dụ phải được viết gần đúng dưới dạng chính tắc. Bài thu hoạch của mỗi nhóm sẽ được trình bày trong giờ lý thuyết hoặc bài tập. Sau khi trình bày thành công trên lớp xong phải hoàn chỉnh thành một bài viết.

44

Tài liệu tham khảo [1] Richard L. Burden and J. Douglas Faires, Numerical Analysis, Ninth Edition, 2010. [2] Alexander Stanoyevitch, Introduction to MATLAB with Numerical Preliminaries, Wiley-Interscience, 2005.

45