Giải tích số [tài liệu giảng dạy ở bậc đại học]

92 trang vanle 3970

Download

Bạn đang xem 20 trang mẫu của tài liệu "Giải tích số [tài liệu giảng dạy ở bậc đại học]", để tải tài liệu gốc về máy bạn click vào nút DOWNLOAD ở trên

Tài liệu đính kèm:

giai_tich_so_tai_lieu_giang_day_o_bac_dai_hoc.pdf

Nội dung text: Giải tích số [tài liệu giảng dạy ở bậc đại học]

TRƯỜNG ĐẠI HỌC THỦY LỢI GIẢI TÍCH SỐ [Tài liệu giảng dạy ở bậc đại học] Nguyễn Thị Vinh HÀ NỘI 2010 0
1 CHƯƠNG 1: MỞ ĐẦU 4 1.1 GIẢI TÍCH SỐ LÀ GÌ 4 1.2 SỰ KHÁC BIỆT GIỮA TOÁN HỌC LÍ THUYẾT VÀ TOÁN HỌC TÍNH TOÁN 4 1.3 CÁC BƯỚC GIẢI MỘT BÀI TOÁN CỦA GIẢI TÍCH SỐ 5 1.4 THUẬT TOÁN VÀ ĐỘ PHỨC TẠP 6 1.4.1 Thuật toán 6 1.4.2 Độ phức tạp thuật toán 7 1.5 SỐ XẤP XỈ VÀ SAI SỐ 10 1.5.1 Số xấp xỉ, sai số tuyệt đối và sai số tuong đối 10 1.5.2 Cách viết số xấp xỉ 11 1.5.3 Qui tròn số và sai số qui tròn 11 1.5.4 Các công thức tính sai số 12 1.6 BÀI TẬP CHƯƠNG 1 13 2 CHƯƠNG 2: GIẢI HỆ PHƯƠNG TRÌNH TUYẾN TÍNH 14 2.1 PHƯƠNG PHÁP KHỬ GAUSS - JORDAN 14 2.1.1 Thuật toán 14 2.1.2 Ưu, nhược điểm của phương pháp 14 2.1.3 Các ví dụ 14 2.1.4 Sơ đồ khối và chương trình 16 2.1.5 Đánh giá độ phức tạp thời gian 17 2.1.6 Ứng dụng phương pháp khử Gauss vào việc tính định thức 17 2.2 GIẢI HỆ PTTT DẠNG BA ĐƯỜNG CHÉO 18 2.2.1 Đặt vấn đề 18 2.2.2 Áp dụng phương pháp khử Gauss–Jordan: 18 2.2.3 Phương pháp truy đuổi (nắn thẳng) giải hệ ba đường chéo 19 2.3 PHƯƠNG PHÁP LẶP SEIDEL 21 2.3.1 Thuật toán 21 2.3.2 Điều kiện hội tụ và đánh giá sai số của phương pháp 21 2.3.3 Ví dụ 22 2.3.4 Sơ đồ khối và chương trình 24 2.3.5 Sử dụng Solver trong EXCEL giải hệ PTTT 26 2.4 TÍNH MA TRẬN NGHỊCH ĐẢO 26 2.4.1 Ứng dụng phương pháp Gauss tính ma trận nghịch đảo 26 2.4.2 Tính ma trận nghịch đảo A–1 bằng phương pháp lặp Newton 27 2.4.3 Sử dụng hàm MINVERSE trong EXCEL tìm A-1 29 2.5 BÀI TẬP CHƯƠNG 2 31 3 CHƯƠNG 3: PHÉP NỘI SUY VÀ ĐƯỜNG CONG PHÙ HỢP 32 3.1 KHÁI QUÁT VỀ BÀI TOÁN NỘI SUY 32 3.1.1 Đặt vấn đề 32 3.1.2 Đa thức nội suy 32 3.1.3 Sơ đồ Horner tính giá trị của đa thức 33 3.2 ĐA THỨC NỘI SUY LAGRANGE 33 3.2.1 Lập công thức 33 3.2.2 Ví dụ: Tìm giá trị gần đúng của f(2,6) từ bảng số liệu 34 1
3.2.3 Sai số: Người ta đã chứng minh rằng nếu hàm f(x) khả vi liên tục đến cấp N+1 trên đoạn [a,b] chứa tất cả các mốc nội suy xk, k = 0, , N thì sai số của nội suy Lagrange là 34 3.2.4 Sơ đồ khối và chương trình 35 3.3 ĐA THỨC NỘI SUY NEWTON VỚI BƯỚC CÁCH ĐỀU 36 3.3.1 Bảng sai phân hữu hạn 36 B ảng sai phân hữu hạn 36 3.3.2 Đa thức nội suy Newton tiến 37 3.3.3 Đa thức nội suy Newton lùi 38 3.3.4 Công thức nội suy Newton với mốc quan sát bất kỳ 41 3.4 NỘI SUY SPLINE 43 3.4.1 Đặt vấn đề 43 3.4.2 Bài toán 43 3.4.3 Xây dựng công thức 43 3.4.4 Các bước giải bài toán nội suy Spline bậc ba 45 3.4.5 Ví dụ 45 3.4.6 Chương trình tính 45 3.5 PHƯƠNG PHÁP BÌNH PHƯƠNG BÉ NHẤT LÀM KHỚP DỮ LIỆU 46 3.5.1 Đặt vấn đề: 46 3.5.2 Lập công thức 47 3.5.3 Các ví dụ: 47 3.5.4 Các bước giải và chương trình 49 3.6 BÀI TẬP CHƯƠNG 3 50 4 CHƯƠNG 4: TÍNH ĐẠO HÀM VÀ TÍCH PHÂN XÁC ĐỊNH 51 4.1 TÍNH GẦN ĐÚNG ĐẠO HÀM 51 4.1.1 Xấp xỉ giá trị đạo hàm dựa vào bảng sai phân 51 4.1.2 Xấp xỉ đạo hàm bằng công thức nội suy 52 4.2 TÍNH GẦN ĐÚNG TÍCH PHÂN XÁC ĐỊNH 55 4.2.1 Lập công thức chung sử dụng đa thức nội suy Newton tiến 55 4.2.2 Quy tắc làm tăng độ chính xác của việc tính tích phân 59 4.3 BÀI TẬP CHƯƠNG 4 61 5 CHƯƠNG 5: GIẢI PHƯƠNG TRÌNH f(x) = 0 62 5.1 ĐẶT VẤN ĐỀ 62 5.1.1 Bài toán 62 5.1.2 Các bước giải 62 5.1.3 Tách nghiệm 62 5.2 CÁC PHƯƠNG PHÁP KIỆN TOÀN NGHIỆM 63 5.2.1 Phương pháp chia đôi 63 5.2.2 Phương pháp lặp đơn 64 5.2.3 Phương pháp dây cung 65 5.2.4 Phương pháp tiếp tuyến (Newton) 67 5.3 GIẢI HỆ PHƯƠNG TRÌNH PHI TUYẾN 69 5.3.1 Lập công thức: 69 Cho hệ phi tuyến 69 5.3.2 Các bước giải hệ phi tuyến bằng phương pháp lặp Newton-Raphson 70 2
5.3.3 Sơ đồ khối và chương trình 73 5.4 PHƯƠNG PHÁP LẶP SEIDEL 74 5.5 Sử dụng Solver trong EXCEL giải hệ phương trình phi tuyến 76 5.6 BÀI TẬP CHƯƠNG 5 77 6 CHƯƠNG 6: CÁC PHƯƠNG PHÁP SỐ GIẢI PHƯƠNG TRINH 78 VI PHÂN 78 6.1 ĐẶT VẤN ĐỀ 78 6.1.1 Bài toán Cauchy (bài toán giá trị đầu) 78 6.1.2 Bài toán biên hai điểm tuyến tính đối với PTVP cấp hai: 79 6.1.3 Các phương pháp số giải bài toán Cauchy 79 6.2 CÁC PHƯƠNG PHÁP SỐ GIẢI BÀI TOÁN CAUCHY 79 6.2.1 Phương pháp Euler 79 6.2.2 Phương pháp Euler cải tiến 81 6.2.3 Phương pháp Runge-Kutta 83 6.2.4 Giải bài toán Cauchy của hệ PTVP cấp một 86 6.3 PHƯƠNG PHÁP SAI PHÂN GIẢI BÀI TOÁN BIÊN TUYẾN TÍNH 87 6.3.1 Xét bài toán biên hai điểm tuyến tính đối với PTVP cấp hai: 87 6.3.2 Ví dụ: Tìm hàm y(x) trên [0; 1] với bước h = 0,1 là nghiệm của 6.3.3 Sơ đồ khối 89 6.4 BÀI TẬP CHƯƠNG 6 90 3
1 CHƯƠNG 1: MỞ ĐẦU 1.1 GIẢI TÍCH SỐ LÀ GÌ Giải tích số (Numerical Analysis) hay còn gọi là Phương pháp số (Numerical Methods) hay Phương pháp tính (Calculating Methods) là một khoa học nghiên cứu các lời giải số của các bài toán của toán học. Ba nhiệm vụ chính của giải tích số là: 1. Xấp xỉ hàm số: Thay một hàm có dạng phức tạp bằng một hàm hoặc nhiềa hàm có dạng đơn giản hơn. Các bài toán thường gặp là nội suy và xấp xỉ hàm. 2. Giải gần đúng các phương trình: Bao gồm các phương trình đại số và siêu việt, các hệ phương trình đại số tuyến tính và phi tuyến, giải các phương trình và hệ phương trình vi phân thường và vi phân đạo hàm riêng, 3. Giải các bài toán tối ưu. Tuy nhiên trong các giáo trình Giải tích số, người ta chỉ đề cập đến hai nhiệm vụ đầu, còn nhiệm vụ thứ ba dành cho các giáo trình về Qui hoạch toán học hay Tối ưu hoá. “An approximate answer to the right problem is worth a great deal more than a precise answer to the wrong problem.” (John W Turkey 1915-2000) 1.2 SỰ KHÁC BIỆT GIỮA TOÁN HỌC LÍ THUYẾT VÀ TOÁN HỌC TÍNH TOÁN TOÁN HỌC LÍ THUYẾT TOÁN HỌC TÍNH TOÁN Chứng minh sự tồn tại nghiệm Tốc độ hội tụ của nghiệm Khảo sát dáng điệu của nghiệm Sự ổn định của thuật toán Một số tính chất định tính của nghiệm Thời gian tính toán trên máy và dung lượng b nhớ cần sử dụng Ví dụ 1: Tính tích phân 1 n x− 1 In = ∫ x e dx (n≥ 1) 0 Tích phân từng phần ta được 1 1 n x− 1 n− 1 x − 1 In− 1 = x e− n x e dx = 1- nI1-n (1.1) 0 ∫ 0 1 1 1 x− 1 x− 1 1-x 1- I1 = xe dx= xe - e dx= e ≈ 0,36787 ∫ 0 ∫ 0 0 Vậy ta có thể tính tích phân trên và để ý rằng In ≥ 0 với mọi n. Trên thực tế không phải như vậy! Công thức trên cho kết quả không chính xác, khi n = 9 I9 ≈ =0,068480 < 0. dù ta tăng độ chính xác của e-1 dến bao nhiêu đi nữa! Nguyên nhân là do sai số ban đầu mắc phải khi tính e-1 bị khuếch đại lên sau mỗi lần tính. 4
Để khắc phuc, ta biến đổi công thức (1.1) -1 In-1 = n (1 – In) và để ý rằng 1 n 1- 0≤ In ≤ x dx = (1 + n) ⇒ lim In = 0 ∫ n→∞ 0 Giả sử ta cần tính I19 và cho I20 ≈ 0 thì sai số mắc phải là ε20 < 1/21. Khi đó I19 ≈ 1/20 với sai số ε19 < 1/21 x 1/20 , , và đến I9 ≈ 0,091623 Ví dụ 2: Giải hệ phương trình đại số tuyến tính AX = b (1.2) với A là ma trận vuông cấp n không suy biến (detA ≠ 0), b là véc tơ cột n thành phần. Do vậy có thể giải theo qui tắc Cramer: xi = ∆i / ∆ , i = 1, , n (1.3) Để tính nghiệm theo (1.3), ta phải tính (n+1) định thức cấp n. Mỗi định thức có n! số hạng, mỗi số hạng có n thừa số, do đó phải thực hiện n-1 phép nhân. Vậy riêng số phép nhân phải thực hiện đã là (n+1) n! (n-1). Giả sử n = 20, và máy tính thực hiện được 5000 phép nhân trong 1 giây thì để thực hiện số phép nhân trên phải mất 300 000 000 năm! Ví dụ 3: Xét hệ (1.2) với ⎛ 0, 1 0 0 ⎞ ⎜ ⎟ ⎜ 0 0, 1 0 ⎟ A = , n= 100 ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ 0 0 0 , 1⎠ Khi đó detA = 10-100 ≈ 0. Theo quan điểm lí thuyết thi A bị suy biến nên không giải được. Tuy nhiên ta có thể nhẩm ra ngay A-1 = 10 E, với E là ma trận đơn vị cấp n và dễ dàng tính được nghiệm của hệ bằng phương pháp ma trận nghịch đảo! Kết luận: Trong quá trình giải các bài toán, có thể nảy sinh nhiều vấn đề mà toán học lí thuyết không quan tâm hoặc không giải quyết được. Vậy cần có một khoa học riêng chuyên nghiên cứu các phương pháp số giải các bài toán. Đó là khoa học tính toán mà giải tích số là một môn học 1.3 CÁC BƯỚC GIẢI MỘT BÀI TOÁN CỦA GIẢI TÍCH SỐ Để giải quyết một bài toán, người ta phải thực hiện quá trình mô phỏng sau đây: 1. Xây dựng mô hình toán học của bài toán thực tế. 2. Phân tích mô hình: tính tương thích của mô hình với bài toán thực tế, sự tồn tại nghiệm của bài toán. 3. Rời rạc hoá mô hình: dùng các phương pháp tính toán để qui bài toán liên tục về bài toán với số ẩn hữu hạn. 4. Xây dựng thuật toán: chú ý đến độ phức tạp thuật toán, tính hội tụ, tính ổn định của thuật toán 5
5. Cài đặt chương trình và chạy thử, kiểm tra kết quả, sửa lỗi và nâng cấp chương trình 1.4 THUẬT TOÁN VÀ ĐỘ PHỨC TẠP 1.4.1 Thuật toán ¾ Khái niệm thuật toán. Thuật toán là một dãy hữu hạn các bước, mỗi bước mô tả chính xác các phép toán hoặc hành động cần thực hiện, để giải quyết một vấn đề. Năm đặc trưng của thuật toán: • Input. Mỗi thuật toán cần có một số (có thể bằng không) dữ liệu vào (input). Đó là các giá trị cần đưa vào khi thuật toán bắt đầu làm việc. Các dữ liệu này cần được lấy từ các tập hợp giá trị cụ thể nào đó. • Output. Mỗi thuật toán cần có một hoặc nhiều dữ liệu ra (output). Đó là các giá trị có quan hệ hoàn toàn xác định với các dữ liệu vào và là kết quả của sự thực hiện thuật toán. • Tính xác định. Mỗi bước của thuật toán cần phải được mô tả một cách chính xác, chỉ có cách hiểu duy nhất • Tính khả thi. Tất cả các phép toán có mặt trong các bước của thuật toán phải đủ đơn giản. • Tính dừng. Mọi bộ dữ liệu vào thoả mãn các điều kiện của dữ liệu vào (tức là được lấy ra từ tập giá trị của các dữ liệu vào), qua xử lí bằng thuật toán phải dừng sau một số hữu hạn bước thực hiện. ¾ Các vấn đề liên quan đến thuật toán. • Tính đúng đắn của thuật toán Khi thuật toán được tạo ra chúng ta cần phải chứng minh rằng thuật toán được thực hiện sẽ cho ta kết quả đúng với mọi dữ liệu vào hợp lệ. Điều này gọi là chứng minh tính đúng đắn của thuật toán. • Các vấn đề nảy sinh Khi giải một bài toán người ta có thể xét nhiều thuật toán khác nhau xem độ phức tạp của chúng ra sao, dùng ngôn ngữ lập trình nào hay cài đặt phần mềm nào để chạy chương trình, cấu trúc dữ liệu nào phù hợp? • Các yêu cầu cơ bản khi giải một bài toán 1 Hiểu đúng bài toán 2 Tìm đúng thuật toán 3 Không nhầm lẫn trong lập trình 4 Dữ liệu quét hết các trường hợp 5 Tốc độ tính toán nhanh 6 Bộ nhớ phù hợp 7 Phần mềm dễ sử dụng, dễ nâng cấp theo yêu cầu 6
1.4.2 Độ phức tạp thuật toán ¾ Định nghĩa Độ phức tạp thuật toán là một công cụ đo, so sánh, lựa chọn các thuật toán khác nhau để tìm ra thuật toán tốt nhất cho lời giải bài toán ¾ Hai tiêu chuẩn để đánh giá độ phức tạp thuật toán - Độ phức tạp về thời gian tính toán - Độ phức tạp về phạm vi bộ nhớ dùng cho thuật toán (dung lượng của không gian nhớ cần thiết để lưu trữ dữ liệu vào, các kết quả tính toán trung gian và các kết quả của thuật toán) Chúng ta sẽ chỉ quan tâm đến thời gian thực hiện thuật toán, vì vậy một thuật toán có hiệu quả được xem là thuật toán có thời gian chạy ít hơn các thuật toán khác. ¾ Cách đánh giá thời gian thực hiện thuật toán - Cỡ của các dữ liệu vào. - Chương trình dịch để chuyển chương trình nguồn thành mã máy. - Thời gian thực hiện các phép toán của máy tính được sử dụng để chạy chương trình. Thời gian chạy chương trình phụ thuộc vào rất nhiều nhân tố, nên ta không thể biết được chính xác thời gian chạy là bao nhiêu đơn vị thời gian chuẩn(bao nhiêu giây) Thời gian thực hiện thuật toán như là hàm số của cỡ dữ liệu vào n. Cỡ của dữ liệu vào là một tham số đặc trưng cho dữ liệu vào, nó ảnh hưởng quyết định đến thời gian thực hiện chương trình. Cỡ dữ liệu phụ thuộc vào thuật toán cụ thể, thường là số nguyên dương n. Ta sử dụng hàm T(n), trong đó n là cỡ dữ liệu vào để biểu diễn thời gian thực hiện một thuật toán. Khi đánh giá thời gian thực hiện bằng phương pháp toán học, chúng ta sẽ bỏ qua nhân tố phụ thuộc vào cách cài đặt chỉ tập trung vào xác định độ lớn của thời gian thực hiện T(n). Kí hiệu toán học ô lớn được sử dụng để mô tả độ lớn của hàmT(n). Giả sử n > 0 (n ∈ Z), T(n),f(n) > 0. Ta viết T(n) = O(f(n)), nếu và chỉ nếu tồn tại các hằng số dương c và n0 sao cho T(n) ≤ c f(n), với mọi n > n0. Ta có thể xem rằng O(f(n)) là cận trên của T(n). Thông thường thời gian chạy một thuật toán tỉ lệ với 1, logn, n, nlogn, nk, hoặc an với a là hằng số ¾ Các quy tắc để đánh giá thời gian thực hiện thuật toán • Chia độ phức tạp thuật toán ra thành nhiều đoạn mà mỗi đoạn có độ phức tạp T1 (n)= O(f1 (n) ⎫ ⎪ ⎬ ⇒ T1(n) + + Tq(n) = O(max(f1(n), ,fq(n))) ⎪ Tq (n)= O(fq (n) ⎭ Thật vậy vì T1(n), ,Tq(n) là ô lớn của f1(n), , fq(n) tương ứng do đó tồn tại hằng số c1, , cq, n1, , nq sao cho T1(n) ≤ c1 f1(n) với mọi n > n1 và T2(n) ≤ c2 f2(n) với mọi n > n2, 7
Đặt n0 = max(n1,n2, , nq), khi đó với mọi n > n0, ta có T1(n) +T2(n) + + Tq(n) ≤ (c1+c2+ + cq) max (f1(n),f2(n), , fq(n)). • Các loai lệnh - Các phép gán, đọc, in, goto là câu lệnh. Các lệnh này được gọi là lệnh đơn và có độ phức tạp thời gian là T(n) = O(C) = O(1) - Nếu S1, S2, .Sn là câu lệnh thì { S1; S2; ; Sn } là câu lệnh hợp thành (câu lệnh ghép khối) - Nếu S, S1, S2, , Sn là các câu lệnh và E1, E2, là biểu thức logic thì if (E1) S; else if (E2) S1; else Sn; gọi là câu lệnh if - Nếu S1, S2, .Sn là câu lệnh, E là biểu thức có kiểu thứ tự đếm được, và v1, v2 , vn là các giá trị cùng kiểu với E thì switch case (constant) { v1: S1; break; v2: S2; break; vn : Sn; break; default: Sn+1; } lệnh này đựoc gọi là câu lệnh case - Nếu S là câu lệnh và E là biểu thức logic thì while (E) S; là câu lệnh while - Nếu S là câu lệnh và E là biểu thức logic thì do S while (E); là câu lệnh do while - Nếu S là câu lệnh, E là biến kiểu thứ tự đếm được thì for (i = E ; i <= n; i++) S; là câu lệnh for • Đánh giá thời gian thực hiện một chương trình 8
Nhờ định nghĩa đệ quy các lệnh, chúng ta có thể phân tích một chương trình xuất phát từ các lệnh đơn, rồi từng bước đánh giá các lệnh phức tạp hơn, cuối cùng đánh giá được thời gian thực hiện chương trình. Giả sử lệnh gán không chứa lời gọi các hàm có thời gian thực hiện. Khi đó để đánh giá thời gian thực hiện một chương trình, ta có thể áp dụng phương pháp đệ quy sau đây: - Thời gian thực hiện các lệnh đơn: gán, đọc, viết, goto là O(1). - Lệnh hợp thành. Thời gian thực hiện lệnh hợp thành được xác định bởi luật tổng. - Lệnh if: Giả sử thời gian thực hiện các lệnh S1, S2 là O(f1(n)) và O(f2(n)) tương ứng. Khi đó thời gian thực hiện lệnh if là O(max(f1(n),f2(n))). - Lệnh case được đánh giá như lệnh if. - Lệnh while: Giả sử thời gian thực hiện lệnh S (thân của lệnh while) là O(f(n). Giả sử g(n) là số tối đa các lần thực hiện lệnh, khi thực hiện lệnh while. Khi đó thời gian thực hiện lệnh while là O(f(n)g(n)). - Lệnh for được đánh giá tương tự lệnh while. Nếu lệnh gán có chứa các lời gọi hàm thì thời gian thực hiện nó không thể xem là O(1) được, vì khi nó thực hiện lệnh gán còn phụ thuộc vào thời gian thực hiện các hàm có trong lệnh gán. việc đánh giá các hàm không đệ quy được tiến hành bằng cách áp dụng các quy tắc trên (việc đánh giá thời gian thực hiện các hàm đệ quy sẽ khó khăn hơn nhiều). Chú ý: - Nếu thuật toán được chia thành nhiều đoạn thì trong trường hợp Đoạn chia có chu trình lồng thắt: độ phức tạp tính là tích Đoạn rời rạc: độ phức tạp tính là max - Người ta chia các bài toán thành 3 lớp, có Độ phức tạp là hàm đa thức (của cỡ dữ liệu vào n) Độ phức tạp là hàm mũ Độ phức tạp là NP - đầy đủ - Hai loại bài toán sau là không khả thi, đối với loại đầu nên giảm bậc của đa thức về càng gần 1 càng tốt. Để hạ bậc của đa thức, người ta thường sử dụng cấu trúc dữ liệu đầu vào (input). Các dữ liệu đầu vào là bình đẳng với nhau. Ví dụ 1: Lệnh lăp for (i = 1; i <= n-1; i++) for (j = 1; j <= n–1; j++) S = S + x[i,j]; có T(n) = O(Cnn) = O(Cn2) Ví dụ 2: Đánh giá thời gian thực hiện của 2 thuật toán tính max của dãy a[1], , a[n] max:=a[1]; for(i=1; i<=n-1; i++) max=a[1]; for(j=i+1; j<=n; j++) for(j=2; j<=n; j++) 9
if(a[i]>a[j]) if (a[j]>max) max=a[j]; { tg=a[i]; a[i]=a[j]; a[j]=tg; } max=a[n]; Thuật toán này làm thay đổi dãy đã cho Thuật toán này không làm thay đổi dãy đã cho Số phép so sánh ≤ (n –1)2 Số phép so sánh = n –1 Số phép gán ≤ 3 (n –1)2 Số phép gán ≤ n T(n) = O(Cn2) T(n) = O(Cn) 1.5 SỐ XẤP XỈ VÀ SAI SỐ 1.5.1 Số xấp xỉ, sai số tuyệt đối và sai số tuong đối ¾ Số xấp xỉ Định nghĩa 1: Số a được gọi là số xấp xỉ của số A nếu a sai khác A không đáng kể và được dùng thay A trong tính toán. Ví dụ: 3,14 và 3,15 là 2 số xấp xỉ của số π ¾ Sai số tuyệt đối và sai số tương đối Định nghĩa 2: |A - a| gọi là sai số tuyệt đối của số xấp xỉ a Nhận xét 1: Do không bết số đúng A nên không biết được |A – a|, người ta đưa thêm vào khái niệm sai số tuyệt đối giới hạn như sau: Định nghĩa 3: Sai số tuyệt đối giới hạn của số a là số Δa: |A-a | ≤ Δa. Vậy. a – Δa ≤ A ≤ a + Δa hay A = a ± Δa (1.4) Ví dụ: Ta biết rằng 3,14 ≤ π ≤ 3,15 nên ta cố thể chọn Δπ ≤ 0,01. Nếu đê ý rằng 3,14 < π < 3,15 thì ta nhận được giá trị tốt hơn. Người ta thường chọn Δπ là số nhỏ nhất thỏa mãn (1.4). Nhận xét 2: Sai số tuyệt đối không thể hiện mức độ chính xác của phép đo hoặc tính toán. Để thể hiện điều đó, người ta đưa vào khái niệm sau: Định nghĩa 4: Sai số tương đối của số xấp xỉ a là số | A - a | δ = | A | Định nghĩa 5: Sai số tương đối giới hạn của số xấp xỉ a là số δa: δ ≤ δa Δa Nghĩa là δa ≥ hay Δa≤ A δ a (1.5) | A | Nhận xét 3: Do không biết số đúng A và vì a ≈ A nên thay vì dùng (1.5), ta dùng ∆a = |a| δa hay A = a (1 ± δa) (1.6) Ví dụ: 10
Mảnh vải A đo được chiều dài a = 6 m với sai số tuyệt đối giới hạn ∆a = 0,2m. Mảnh vải B đo được chiều dài b = 3 m với sai số tuyệt đối giới hạn ∆b = 0,2m. Rõ ràng 2 mảnh vải có cùng sai số tuyệt đối giới hạn là 0,2m và A = 6 ± 0,2; B = 3 ± 0,2 nhưng sai số tương đối giới hạn của phép đo mảnh vải A là Δa 2,0 1 1 δa = = = = δb a 6 30 2 Vậy phép đo mảnh vải A chính xác hơn. 1.5.2 Cách viết số xấp xỉ ¾ Chữ số có nghĩa của một số: Chữ số có nghĩa của một số là những chữ số của số đó kể từ chữ số khác không đầu tiên tính từ trái sang phải Ví dụ: 0,0052 có 2 chữ số có nghĩa và 23,05 có 4 chữ số có nghĩa ¾ Chữ số đáng tin: Chữ số đáng tin của một số thực là những chữ số ở vị trí thứ i (i = n, n-1, ,1, 0, -1, -2, , m) thoả mãn điều kiện: ∆a ≤ 0,5 10 i-n+1 với n ≥ 0, m ≤ 0 là các số nguyên. 2 -1 Ví dụ: Số a = 325,4 = 3 x 10 + 2 x 10 + 5 + 4 x 10 có ∆a = 0,01 ≤ 0,5 x 10 -1-2+1 = 0,5 x 10-2 Vậy chữ số 4 ở hàng thập phân và các chữ số ở bên trái nó là các chữ số đáng tin. Qui ước: nếu viết số gần đúng mà không kèm theo sai số thì mọi chữ số đều là chữ số đáng tin ¾ Cách viết số xấp xỉ: Cách 1: A = a ± ∆a Cách 2: mọi chữ số có nghĩa đều đáng tin, tức là ∆a ≤ 0,5.10-n với n là vị trí của chữ số cuối cùng nằm bên phải, kể từ dấu phẩy. Ví dụ: Số a = 3,14 có chữ số 4 là chữ số đáng tin và do đó các chữ số đứng trước nó đều đáng tin và số này có sai số tuyệt đối ∆a ≤ 0,5.10-2 1.5.3 Qui tròn số và sai số qui tròn Qui tròn một số là vứt bỏ một vài chữ số sao cho số nhận được a1 được gần đúng nhất với số a ban đầu ¾ Qui tắc qui tròn Sai số qui tròn tuyệt đối ≤ 0,5 đơn vị của chữ số ở hàng giữ lại cuối cùng ¾ Sai số qui tròn: Gọi sai số qui tròn tuyệt đối của a là Өa1 = |a1 - a|, ta có |A - a1| = |A – a + a - a1| ≤ ∆a + Өa1 Vậy ta có thể chọn được ∆a1 = ∆a+ Өa1 điều đó có nghĩa là sau qui tròn, sai số tuyệt đối tăng lên Ví dụ: Số a = 3,141 được qui tròn thành 3,14, còn số 3,1415 được qui tròn thành 3,142 11
1.5.4 Các công thức tính sai số ¾ Công thức chung Bài toán: Cho hàm u = f(x1, x2, , xn) khả vi và biết các sai số tuyệt đối giới hạn ∆xi của các đối số xi và giá trị gần đúng u. Hãy tính sai số tuyệt đối ∆u và sai số tương đối δu Giải: Sử dụng các định nghĩa về sai số và công thức toán ta có n ∂u n ∂ Δ = ∑ Δ , δ = ∑ ln u Δ (1.7) u ∂x x u x i= 1 i i i= 1 ∂xi i Áp dụng công thức (1.7) cho các hàm số học ta có các công thức tính sai số sau. ¾ Sai số của tổng đại số: Nếu u = x1 + x2 + + xn thì Δu = Δ + Δ + + Δ x1 x 2 xn Δ + Δ + + Δ Δ x x x δ = u = 1 2 n u u u ¾ Sai số của tích: Nếu u = x1 x2 xn thì δx + δx + + δx , Δu= u δ u 1 2 n ¾ Sai số của thương: Nếu u = x1 / x2 (x2 ≠ 0) thì δu = δ x + δx , Δu= u δ u 1 2 Ví dụ: Tính các sai số tuyệt đối và tương đối của thể tích hình cầu , biết V = πd3/6; d = 3,7 cm ± 0,05 cm và π = 3,14 + 0,0016. Giải: Theo công thức tính sai số của tích và thương ta có 0,0016 0,05 δ= δ +3 δ = +3 =0.0411 ⇒Δ =V. δ =1,0882 ≈ 1,1 (cm3 ) V π d 3,14 3,7 V V 12
1.6 BÀI TẬP CHƯƠNG 1 1> Giải hệ phương trình 0,461x1 + 0,311x2 = 0.150 0,209x1 + 0,141x2 = 0.068 Bằng cách sử dụng phép toán số học làm tròn đến ba chữ số thập phân. Nghiệm chính xác là x1 = 1, x2 = –1; Bạn có so sánh gì giữa hai kết quả? 2> Xét thuật toán tạo nhiễu y cho đại lượng x: A = x × 10n B = A + x y = B – A Tính y khi x = 0,123456 và n = 3 bằng cách sử dụng phép toán số học làm tròn đến sáu chữ số thập phân. Sai số x – y bằng bao nhiêu? 3> Một hình vuông có cạnh a đo được bằng 12,04 (m) với độ chính xác của thước đo là ∆a = 0,005 (m). Đánh giá sai số tương đối, sai số tuyệt đối của diện tích của hình vuông S = a2 và tính ra diện tích này (chỉ giữ lại các chữ số có nghĩa). 4> Giả sử z = 0,180 × 102 là nghiệm gần đúng của phương trình ax = b với a = 0,111 × 100, b = 0,200 × 101. Hãy sử dụng phép tính số học với ba chữ số thập phân để tính sai số ∆ = |b – az|. Sau đó tính sai số trên bằng phép tính số học chính xác và so sánh các kết quả. 13
2 CH ƯƠNG 2: GIẢI HỆ PHƯƠNG TRÌNH TUYẾN TÍNH 2.1 PHƯƠNG PHÁP KHỬ GAUSS - JORDAN 2.1.1 Thuật toán ¾ Bước khử xuôi: Đưa hệ ⎧ a11 x 1+ a 12 x 2 + + a1n x n = b 1 ⎪ ⎪ a21 x 1+ a 22 x 2 + + a2 n x n = b 2 ⎨ (2.1) ⎪ ⎪ ⎩ an1 x 1+ a n 2 x n + + an n x n = b n về hệ dạng tam giác trên bằng các phép biến đổi tương đương qua n lần tác động lên A ()()n n ()n ()n ⎧a11 x1+ a 12 x 2 + +a1n xn = b 1 ⎪ ⎪ a()n x + +a()n x = b()n ⎨ 22 2 2n n 2 (2.2) ⎪ ⎪ ()n ()n ⎩ ann xn= b n Các phép biến đổi trên thực hiện được ở lần tác động thứ i lên ma trận A nếu (i− 1) a ≠0 ∀ i = 1, n ii và là phần tử có trị tuyệt đối lớn nhất trên cột i, cụ thể ta có công thức biến đổi (i) (i− 1) (i − 1) (i-1) (i− 1) ∀i = 1, n; ∀ j = i + 1, n; ∀ k = i, n ajk = ajk − aji aik /aii ¾ Bước quét ngược: Giải lần lượt xn, xn-1, , x1 từ hệ tương đương đã được chéo hóa (2.2) 2.1.2 Ưu, nhược điểm của phương pháp (i− 1) Phương pháp lặp Gauss-Jordan đơn giản, dễ lập trình, tuy nhiên nếu aii ≈ 0 thi kết quả có thể không chính xác. Do vậy ở ma trận biến đổi thứ i ta phải tìm phần tử trội trên cột i và thực hiện phép đổi chỗ 2 hàng nếu cần thiết để đảm bảo hàng thứ i là hàng chứa phần tử trội 2.1.3 Các ví dụ Ví dụ 1: Giải hệ ⎧ 5 x1 − x2 + 2 x3 = 6 ⎪ ⎨ x1 − 4 x2 + x3 = − 2 ⎪ ⎩ −2 x1 − x2 + 4 x3 = 1 14
Bước khử xuôi: Đưa hệ đã cho về hệ tương đương ⎛5− 1 2 6 ⎞ ⎜ ⎟ ⎜0− 3,8 0,6 −3,2 ⎟ ⎜ ⎟ ⎝0 0 4,579 4,579⎠ Bước quét ngược: Giải ra ta được x3 = 1; x2 = 1; x1 = 1 Ví dụ 2: Cho mạch điện như hình vẽ hãy áp dụng định luật Kirhoft để tìm giá trị các dòng điện qua các điện trở trong các mạch nhánh. Các giá trị cho trên hình vẽ: Giải: Áp dụng định luật Kirhoft 1 cho nút C ta có : I1 + I2 - I3 = 0 (2.3) Áp dụng định luật Kirrhoft 2 cho: Vòng ACDB: V1 = R1I1 + R3 I3 (2.4) Vòng AEFB : V1 - V2 = R1I1 - R2 I2 (2.5) Thay các giá trị bằng số đã biết vào các phương trình (2.3), (2.4), (2.5) ta có hệ I1 + I2 - I3 = 0 2I1 +5I3 = 6 2I1 - 4I2 = 4 Áp dụng phương pháp giải Gauss ta tìm được kết quả : I1 = 1,158 ; I2 = -0,421; I3 = 0,737(A) 15
2.1.4 S ơ đồ khối và chương trình Bắt đầu Nhập A Cấu tạo ma trận tam giác (khử xuôi) Tìm nghiệm (khử ngược) In KQ Dừng void Gauss(int n,matran a, mang1 x) { int i, j, k, max; float t; // Buoc khu xuoi for(i = 0; i fabs(a[max][i])) max = j; if (max != i) // Doi cho 2 hang max va i neu can for(k = i; k = i; k––) a[j][k] = a[j][k] – a[i][k] * a[j][i] / a[i][i]; // Buoc quet nguoc if (!(a[n–1][n–1])) { if(a[n–1][n]) printf(" HE VO NGHIEM"); else printf(" HE VO SO NGHIEM"); 16
} else { for(j = n–1; j >= 0; j ) { t = 0.; for (k = j + 1; k < n; k++) t = t + a[j][k] * x[k]; x[j] = (a[j][n] - t) / a[j][j]; } printf(" \n\r Nghiem cua he la: \n"); for (i = 0; i<n; i++) printf("x[%d]= %8.4f\n",i,x[i]); } } 2.1.5 Đánh giá độ phức tạp thời gian Thời gian chạy của bước khử ngược là O(n2) nên hầu hết công việc được làm ở bước khử xuôi. Với mỗi giá trị của vòng lặp theo i, vòng lặp theo k lặp n – i + 2 lần và vòng lặp theo j lặp n – i lần. Vậy số lần lặp tổng cộng là n n 3 ∑ (n− i + 2)(n − i) + O(n2 ) = + O(n2 ) i= 1 3 Giá trị của –a[i][j]/a[i][i]được tính ngoài vòng lặp k nên mỗi lần lặp chỉ bao gồm một phép nhân và một phép cộng. Từ đó ta có kết luận: Một hệ n phương trình, n ẩn giải được bằng phương pháp Gauss với khoảng n3 / 3 phép cộng và nhân. 2.1.6 Ứng dụng phương pháp khử Gauss vào việc tính định thức ¾ Thuật toán: Sử dụng bước khử xuôi, ta biến đổi định thức đã cho về định thức dạng tam giác trên. Từ đó giá trị của định thức chính là tích các phần tử trên đường chéo chính của ma trận vế phải sau n lần tác động, nhân với phần dấu là (–1)p, trong đó p là số lần đổi chỗ 2 hàng. ⎡ 1 0 1 2⎤ ⎢−1 2 3 1⎥ Ví dụ: Tính định thức của ma trận A = ⎢ ⎥ ⎢ 4 0− 2 1⎥ ⎢ ⎥ ⎣ 0 2 1 2⎦ Sau bước khử xuôi ta đưa được A về dạng tam giác trên với p = 2 lần đổi chỗ các hàng ⎡4 0− 2 1 ⎤ ⎢0 2 2,5 1,25⎥ ⎢ ⎥ ⎢0 0 1,5 1,75⎥ ⎢ ⎥ ⎣0 0 0 2,5 ⎦ Lấy tích các phần tử trên đường chéo chính với (–1)2 ta được det(A) = 30 17
¾ Sử dụng hàm MDETERM trong EXCEL tìm det(A) 2.2 GIẢI HỆ PTTT DẠNG BA ĐƯỜNG CHÉO 2.2.1 Đặt vấn đề Khi ma trận vế phải của hệ A X = b có dạng ba đường chéo, một số thuật toán đơn giản giúp giảm được đáng kể thời gian tính toán và dung lượng bộ nhớ ⎡ c1 d 1 0 0 0 0 0 e1 ⎤ ⎢ b c d 0 0 0 0 e ⎥ ⎢ 2 2 2 2 ⎥ ⎢ 0 b3 c 3 d 3 0 0 0 e3 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ 0 0 0 0 cn− 2 d n − 2 0 en− 2 ⎥ ⎢ 0 0 0 0 b c d e ⎥ ⎢ n− 1 n − 1 n −1 n −1⎥ ⎣ 0 0 0 0 0 bn c n en ⎦ 2.2.2 Áp dụng phương pháp khử Gauss–Jordan: Nhận xét: Đối với bước khử xuôi, chỉ có trường hợp j = i+1 và k = i+1 cần phải có mặt vì a[j] [k] = 0 với mọi k > i + 1. Đưa hệ đã cho về dạng 2 đường chéo ở tam giác trên. Đối với ma trận có dạng như vậy, bước khử xuôi và bước quét ngược được rút gọn về 1 vòng lặp: for (i = 0; i = 0; i ) x[i]=(a[i][n]–a[i][i+1]*x[i+1])/a[i][i]; 18
Vậy nếu giải hệ ba đường chéo bằng phương pháp Gauss-Jordan thì độ phức tạp thời gian là tuyến tính. 2.2.3 Phương pháp truy đuổi (nắn thẳng) giải hệ ba đường chéo Dung lượng đòi hỏi để chứa mảng hai chiều có thể được thay thế bằng bốn mảng một chiều, mỗi mảng cho một đường chéo và một mảng cho vec tơ cột vế phải. Sau đó dùng phương pháp truy đuổi giải hệ, độ phức tạp thuật toán là thời gian tuyến tính. ¾ Thuật toán Viết lại phương trình thứ i của hệ bi xi–1 + ci xi + di xi+1 = ei , i = 1, , n , (b1 = 0 và dn = 0) (2.6) và đưa vào quan hệ phụ: xi+1= fi xi + hi, i=1, , n (2.7) Thay (2.7) vào (2.6) ta tính được fi và hi bằng công thức truy hồi: b e − d h f= − i ; h = i i i ∀i = 1, n -1 (2.8) i− 1 i− 1 c + d f ci + d i f i i i i Để tính fn-1 và hn-1 ta xét phương trình cuối cùng bn xn-1 + cn xn = en hằng đẳng hệ số với (2.6) khi i = n, ta tính được bn en cn x n= b n x n− 1 + e n ⇔ x n = −xn− 1 + cn cn bn en hay f n− 1 = − ; h n− 1 = cn cn rồi thay vào (2.8) tính ngược lên. Tính x1: xét phương trình đầu tiên của (2.6) c1 x1 + d1 x2 = e2 mà x2 = f1 x1 + h1. Giải ra ta được e1− d 1 h 1 x1 = c1+ d 1 f 1 Thay vào công thức (2.7) ta tìm được các xi còn lại 19
¾ Sơ đồ khối và chương trình Bắt đầu Nhập b, c, d, Tính các hệ số fi và hi (i = n–1, , 1) Tìm nghiệm xi (i = 1, , n) In KQ Dừng void Tridiagonal(int n,mang1 b,mang1 c,mang1 d,mang1 e,mang1 x) { mang1 f,h; int i; // Tinh cac vec to f va h printf("\n\rTinh cac vecto he so f va h cua nghiem x: \n"); f[n-2]=-b[n-1]/c[n-1]; h[n-2]=e[n-1]/c[n-1]; for (i=n-2;i>0;i ) { f[i-1]=-b[i]/(c[i]+d[i]*f[i]); h[i-1]=(e[i]-d[i]*h[i])/(c[i]+d[i]*f[i]); } for (i=0;i<=n-2;i++) printf(" f[%d] = %10.4f h[%d] = %10.4f",i,f[i],i,h[i]); //Tinh nghiem X printf("\n\r Nghiem cua he la: \n"); x[0]=(e[0]-d[0]*h[0])/(c[0]+d[0]*f[0]); printf(" x[0] = %10.4f \n",x[0]); for (i=1;i<n;i++) { x[i]=f[i-1]*x[i-1]+h[i-1]; printf(" x[%d] = %10.4f \n",i,x[i]); } } 20
2.3 PHƯƠNG PHÁP LẶP SEIDEL 2.3.1 Thuật toán Bước 1: Đưa hệ ⎧ a11 x 1+ a 12 x 2 + +a12 x 2 = b 1 ⎪ ⎪ a21 x 1+ a 22 x 2 + + a2 n x n = b 2 ⎨ (2.9) ⎪ ⎪ ⎩ an1 x 1+ a n 2 x n + + an n x n = b n về hệ tương đương ' ⎧ x1= c 12 x 2 + +c1n x n + b 1 ⎪ ' ⎪ x2= c 21 x 1 + c 23 x 3 + + c2 n x n + b 2 ⎨ (2.10) ⎪ ⎪ ' ⎩ xn= c n1 x 1 + + cn n− 1 x n − 1+ b n bằng cách ở phương trình thứ i chuyển các xk ≠ xi sang phải và chia 2 vế cho aii, i=1, , n (0) (0) (0) Bước 2: Thay xấp xỉ ban đầu X = (x1 , ,xn ) của nghiệm vào (2.10) tìm xấp xỉ thứ (1) (1) (1) nhất X = (x1 , ,xn ) của nghiệm ()1 (0 ) (0) ' ⎧x1 = c12 x 2 + +c1n x n + b1 ⎪ ()1 ()1 (0 ) ()0 ' ⎪x2 = c21 x 1 + c23 x 3 + + c2n x n + b1 ⎨ (2.11) ⎪ ⎪ ()1 ()1 ()1 ' ⎩xn = cn1 x 1 + +cnn− 1 x n − 1 + b1 (1) (1) (1) (1) Ta sử dụng các x1 , , xi được tính ở các hàng trên để tính các xi+1 , ,xn , với i = 1, , n –1 ở các hàng dưới . (2) (2) (2) Bước 3: Tìm xấp xỉ thứ 2 X = (x1 , ,xn ) của nghiệm như bước 2 bằng cách thay (0) (1) (1) (2) tương ứng xi bởi xi ,xi bởi xi , i=1, ,n. Tiếp tục như vậy cho đến bước lặp thứ i = M cho trước hoặc đến khi số bước lặp i thoả mãn X(i) − X * ≤ ε , ε > 0 cho trước (2.12) ∞ 2.3.2 Điều kiện hội tụ và đánh giá sai số của phương pháp • Điều kiện đủ để phương pháp lặp Gauss-Seidel hội tụ: Nếu ma trận A có tính chất đường chéo trội, tức là n aii > ∑ aij , i= 1 , , n (2.13) j= 1, j ≠ i 21
thì nghệm X(n) hội tụ đến nghiệm duy nhất X* mà không phụ thuộc vào việc chọn véc tơ X(0) ban đầu. • Đánh giá sai số của phương pháp: Người ta chứng minh được rằng μ XX()k −* ≤ XX()()k− k− 1 (2.14) ∞ 1− μ ∞ μk XX()k −* ≤ XX()()1− 0 ∞ 1− μ ∞ trong đó X = max xi , ∞ i qi μ = max , i= 1, 2, , n, i 1− pi i− 1 n pi= ∑ c ij ; qi= ∑ c ij (2.15) j= 1 j= i 2.3.3 Ví dụ Giải hệ sau với sai số ε = 0,05 và M = 10. ⎧ 5x1 − x2 + 2x3 = 6 ⎪ ⎨ x1− 4x 2 + x3 = − 2 ⎪ ⎩−2x1 − x2 + 4x3 = 1 Do 5> − 1 + 2 −4 > 1 + 1 4> − 2 + − 1 ⎧x =0,2x −0,4x + 1,2 ⎪ 1 2 3 ⎪ ⇒ x =0,25x +0,25x + 0,5 ⎨ 2 1 3 ⎪ x =0,5x +0,25x +0,25 ⎩⎪ 3 1 2 nên A thoả mãn tính chất đường chéo trội, vì vậy phép lặp Seidel hội tụ Bước 1: Đưa hệ đã cho về hệ tương đương Có p1=0, q1=0,6; q1/(1 − p1)=0,6 p2=0,25; q2=0,25; q2/(1 − p2)=1/3 p3=0,75; q3=0; q3/(1 − p3)=0 22
μ Vậy μ=0,6 ⇒ =1,5 1− μ Bước 2: Thay xấp xỉ ban đầu, chẳng hạn X(0) = (0 0 0) của nghiệm (1) (1) (1) (1) vào vế phải tìm xấp xỉ thứ nhất X = (x1 ,x2 ,x3 ) của nghiệm theo (2.11) với n=3 (1) x1 = 0*0,2 - 0*0,4 + 1,2 = 1,2 (1) x2 = 0,25*1,2 + 0*0,25 + 0,5 = 0,8 (1) x3 = 0,5*1,2 + 0,25*0,8 + 0,25 = 1,05 (2) (2) (2) (2) Bước 3: Tìm xấp xỉ thứ hai X = (x1 , x2 , x3 ) của nghiệm (1) (0) max x −x = max ( 1,2; 0,8; 1,05 )= 1,2 i i i (1) Vậy X −X* ≤ 1,5*1,2 = 1,8 (2) x1 = 0,2*0,8 - 0,4*1,05 + 1,2 = 0,94 (2) x2 = 0,25*0,94 + 0,25*1,05 + 0,5 = 0,998 (2) x3 = 0,5*0,94 + 0,25*0,998 + 0,25 = 0,97 (2) (1) max x −x = max ( 0,26; 0,198; 0,08 )= 0,26 i i i (2) Vậy X −X* ≤ 1,5*0,26= 0,39 (3) (3) (3) (3) Bước 4: Tìm xấp xỉ thứ ba X = (x1 , x2 , x3 ) của nghiệm (3) x1 = 0,2*0,998 - 0,4*0,97 +1,2 = 0,95 (3) (2) max x −x = max ( 0,01; 0,012; 0 )= 0,012 i i i (3) Vậy X −X * ≤ 1,5* 0,012= 0,018 < ε (3) x2 = 0,25*0,95 + 0,25*0,97 + 0,5 = 1,01 (3) x3 = 0,5*0,95 + 0,25*1,01 + 0,25 = 0,97 Ta dừng ở xấp xỉ thứ ba của nghiệm vì sai số đã thỏa mãn độ chính xác ε cho trước 23
2.3.4 S ơ đồ khối và chương trình Bắt đầu Nhập A Kiểm tra tính chéo trội của A s đ Cấu tạo ma trận lặp Giải lặp Seidel In KQ Dừng enum BOOLEAN Hoitu (int m,matran a) { int i, j; float s; enum BOOLEAN TL, KT, KTH[nmax]; TL=true; i = 0; do { KT = true; KTH[i] = true; if (a[i][i]) { s = 0.; for (j=0;j = 1.) KT = false; KTH[i] = KTH[i] && KT; } else KTH[i] = false; TL = TL && KTH[i]; i++; } while ((i <= m –1) && (TL==true)); return TL; } void Matranlap (int m, matran a, matran b) { 24
int i, j; for (i = 0; i <= m –1; i++) { b[i][m] = a[i][m] / a[i][i]; for (j = 0; j <= m –1; j++) if (i == j) b[i][j] = 0.; else b[i][j] = –a[i][j] / a[i][i]; } } void LapSeidel(int m, int lanlap, float epsilon, matran b, mang1 xs) { mang1 xt; int i, j, k; double s; enum BOOLEAN t; // Lay vec to ve phai lam vec to nghiem ban dau for (i = 0;i <= m –1; i++) xt[i] = b[i][m]; k = 1; t = false ; do { for (i = 0; i <= m –1; i++) { xs[i] = 0; for (j = 0; j <= m-1; j++) if (j < i) xs[i] = xs[i] + b[i][j] * xs[j]; else xs[i] = xs[i] + b[i][j] * xt[j]; xs[i] = xs[i] + b[i][m]; } printf(" Nghiem lap thu %2d la:\n",k); for (i = 0; i <= m –1; i++) { printf(" x[%d] =%12.6f",i,xs[i]); printf("\n"); } // Kiem tra dieu kien dung s = 0.0 ; for (i = 0; i < m –1; i++) s += fabs(xs[i] – xt[i]) ; if (s <= epsilon) { t = true ; printf("\n"); printf("PHEP LAP HOI TU SAU %d BUOC\n" , k); } else printf("Sai so o lan lap thu %d la %12.2E" ,k,s); for (i = 0; i < m –1;i++) xt[i] = xs[i] ; k++ ; } while((!t) && (k<=lanlap)); } 25
2.3.5 S ử dụng Solver trong EXCEL giải hệ PTTT • Nhập các giá trị đoán nhận ban đầu cho x1, x2, , xn • Nhập các vế trái của các phương trình f1, f2, , fn vào các ô riêng được liên kết với giá trị của các ẩn • Chọn công cụ Solver (Data Tab Æ Solver) a. Đặt “Target Cell” cho ô có công thức = 0 b. Chọn “Equal to Value of 0” c. Đối với “By Changing Cells,” trỏ tới các ô chứa giá trị ban đầu của x1, x2, , xn d. Thêm các ràng buộc (các phương trình) e. Nhấn nút “Options” , chọn “Assume linear model” f. Nhấn nút “Solve” 2.4 TÍNH MA TRẬN NGHỊCH ĐẢO 2.4.1 Ứng dụng phương pháp Gauss tính ma trận nghịch đảo Để tính ma trận nghịch đảo A-1 của ma trận A, ta dựa vào định nghĩa A A–1 = E hay ⎡a11 a 12 a1n ⎤ ⎡x11 x 12 x1n ⎤ ⎡1 0 0 ⎤ ⎢a a a ⎥ ⎢x x x ⎥ ⎢ 0 1 0 ⎥ ⎢ 21 22 2n ⎥ ⎢ 21 22 2n ⎥ = ⎢ ⎥ (2.16) ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣an1 a n2 a nn ⎦ ⎣xn1 x n2 x nn ⎦ ⎣ 0 0 1 ⎦ Để tính các cột của ma trận nghịch đảo ta giải Gauss – Jordan hệ n phương trình có cùng ma trận vế trái, còn các véc tơ cột vế phải lần lượt là các cột của ma trận đơn vị E. 26
⎡1 1 1⎤ ⎢ ⎥ Ví dụ: Tìm ma trận nghịch đảo của ma trận A = ⎢2 3 1⎥ ⎣⎢1 2 1⎦⎥ Các cột của A–1 lần lượt là nghiệm của các hệ phương trình ĐSTT sau đây với cùng ma trận vế phải A ⎡1 1 1 1 0 0⎤ ⎢ ⎥ ⎢2 3 1 0 1 0⎥ ⎣⎢1 2 1 0 0 1⎦⎥ Giải ra ta được ⎡ 1− 1 1 ⎤ −1 ⎢ ⎥ A = ⎢− 1 0 1 ⎥ ⎣⎢ 1 1− 2⎦⎥ 2.4.2 Tính ma trận nghịch đảo A–1 bằng phương pháp lặp Newton ¾ Đặt vấn đề Cho số a ≠ 0, cần tìm số x là số nghịch đảo của a. Đặt f(x) = a – 1/x, áp dụng phương pháp lặp Newton (tiếp tuyến) tìm nghiệm gần đúng của phương trình f(x) = 0 1 a − f(xn ) x n xn+ 1= x n − = xn − = xn (2 − a xn ) (2.17) f'(xn ) 1 2 x n ¾ Thuật toán Áp dụng công thức lặp tương tự (2.17) để tìm ma trận nghịch đảo A-1 của ma trận -1 không suy biến A. Gọi Xk là ma trận gần đúng của A ở bước lặp thứ k, ta có Xn+1 = Xn (2E – AXn) (2.18) với E là ma trận đơn vị cùng cấp. ¾ Sự hội tụ Đặt Gk = E – AXk thay (2.18) vào = E – AXk – 1(2E – AXk - 1) = E – 2AXk – 1 + AXk – 1AXk - 1 2 2 = (E – AXk – 1) = G k-1 2 4 2k ⇒Gk = G k− 1 = Gk − 2 = = G 0 (2.19) Mặt khác do A-1A = E nên -1 -1 -1 -1 2k A – Xk = A (E – AXk) = A Gk = A G0 −1 −1 2k ⇒AXAG −k ≤ 0 (2.20) -1 Vậy nếu xấp xỉ đầu X0 được chọn gần với A , cụ thể là nếu 27
E− AX0 = G 0 < 1 -1 thì từ (2.20) dãy {Xk} hội tụ về A rất nhanh. ¾ Sơ đồ khối và chương trình Bắt đầu Nhập A Cấu tạo ma trận lặp Kiểm tra đ/k hội tụ s đ Tính lặp ma trận nghịch đảo In KQ Dừng void Daomatran(int n,int lanlap,matrix a,matrix xt,matrix xs) { int i, j, k,dem; matrix g,y; double t; FILE *fp; // Tao ma tran lap for (i = 0; i <= n − 1; i++) for (j = 0; j <= n − 1; j++) for (k = 0; k <= n − 1; k++) { if (i==j) { g[i][j] = 1. - a[i][k] * xt[k][j]; y[i][j]=g[i][j]+1; } else { g[i][j] = - a[i][k] * xt[k][j]; y[i][j]=g[i][j]; } } // Kiem tra dieu kien hoi tu cua ma tran lap t=0.; for (i = 0; i <= n − 1; i++) 28
for (j = 0 ; j = 1) { printf("\tPhep lap khong hoi tu! \n"); return; } else //Tinh lap ma tran dao { dem = 1; do { for (i = 0; i <= n − 1; i++) for (j = 0 ; j <= n − 1 ; j++) for (k = 0; k <= n − 1; k++) xs[i][j] = xt[i][k] * y[k][j]; for (i = 0; i <= n − 1; i++) for (j = 0 ; j <= n − 1 ; j++) xt[i][j] = xs[i][j]; dem++; } while (dem <= lanlap); } } 2.4.3 Sử dụng hàm MINVERSE trong EXCEL tìm A-1 1. Chọn mảng các ô chứa A-1. 2. Nhập hàm MINVERSE tính A-1 với các tham số là các ô chứa ma trận A 3. Nhấn phím F2, rồi nhấn tiếp CTRL+SHIFT+ENTER 29
2.5 BÀI TẬP CHƯƠNG 2 1> Giải các hệ phương trình tuyến tính sau bằng phương pháp Gauss-Jordan 2x1 + 5x2 – 2x3 = 3 3x1 + 6x2 + 9x3 = 39 x1 + 3x2 – x3 = 2 xl + 2x2 – x3 = 2 2x1 + 4x2 + x3 = 7 3x1 + 6x2 – 2x3 = 7 2> Sử dụng chương trình con Gauss để giải hệ phương trình tuyến tính sau 0,21 x1 + 0,32 x2 + 0,12 x3 + 0,31x4 = 0,96 0,10 x1 + 0,15 x2 + 0,24 x3 + 0,22x4 = 0,71 0,20 x1 + 0,24 x2 + 0,46 x3+ 0,36x4 = 1,26 0,61 x1 + 0,40 x2 + 0,32 x3 + 0,20x4 = 1,53 Kiểm tra đáp án của bạn bằng cách thay lại vào hệ phương trình ban đầu, và đánh giá độ chính xác. Nghiệm đúng là [l l l l]T 3> Sử dụng chương trình con Tridiagonal để giải hệ phương trình tuyến tính sau 2 –2(1 + h )x1 + x2 = 1 2 –xi–1 – 2(1 + h )xi + xi+1 = 0 i = 2, , n–1 2 xn–1 – 2(1 + h )xn = 1 khi n = 30 và h = 0,1 4> Tìm ma trận nghịch đảo của ma trận sau bằng phương pháp Gauss-Jordan 1 2 3 4 5 6 7 8 9,01 5> Viết chương trình con tính định thức của ma trận vuông cấp n, bằng cách sử dụng bước khử xuôi của phương pháp Gauss-Jordan đưa ma trận ban đầu về dạng tam giác trên. Sử dụng chương trình này, tính định thức của ma trận sau 0,217 0,732 0,414 0,508 0,809 0,376 0,795 0,886 0,338 31
3 CH ƯƠNG 3: PHÉP NỘI SUY VÀ ĐƯỜNG CONG PHÙ HỢP 3.1 KHÁI QUÁT VỀ BÀI TOÁN NỘI SUY 3.1.1 Đặt vấn đề Nhiều bài toán đòi hỏi tìm hàm f(x) từ tập các điểm (xi, yi) cho trước, chẳng hạn như việc thiết kế hệ thống tính toán điều khiển năng lượng cho một toà nhà cao tầng, đòi hỏi số liệu về sự biến đổi nhiệt độ ở đó trong mỗi ngày. Mẫu nhiệt độ thường được đo tại các thời điểm rời rạc. Tuy nhiên khi tính toán lại cần đến các giá trị nhiệt độ không có trong mẫu đo. Một trong các cách khắc phục vấn đề trên là tìm hàm đơn giản, thường là đa thức, đo đúng mẫu và nội suy những giá trị nằm giữa các điểm mẫu. Cũng có trường hợp biểu thức giải tích của f(x) đã biết nhưng quá phức tạp. Khi đó dùng phép nội suy ta có thể tính được giá trị của f(x) tại điểm x bất kì trên [a,b] mà độ chính xác không kém bao nhiêu. Phép nội suy thường dùng trong các ứng dụng của khoa học, công nghệ và kinh tế. Cơ sở của mọi phép nội suy là sự phù hợp của kiểu đường cong hay hàm nội suy đối với tập các điểm mẫu cho dưới dạng bảng. 3.1.2 Đa thức nội suy ¾ Bài toán: Cho tập các điểm thực nghiệm (x0,y0), (x1,y1), ,(xn,yn), hãy xây dựng một đa thức n–1 n Pn(x) = a0 + a1x + + an-1x + anx (3.1) thoả mãn Pn(xi) = yi i = 0, 1, 2, , n (3.2) Các điểm xi, yi i = 0, 1, 2, , n gọi là các mốc nội suy hay các điểm nút nội suy. Về mặt hình học, bài toán nội suy đa thức là tìm một đường cong y = Pn(x) đi qua các điểm M0(x0,y0), M1(x1,y1), , Mn(xn,yn) cho trước Sau đó dùng đa thức này để tính các giá trị Pn(x) tại các giá trị x nằm giữa x0 và xn ¾ Giải: n+1 hệ số ai được xác định từ hệ phương trình ĐSTT n j ∑ aj xi = yi i = 0, n (3.3) j= o Định thức chính của hệ là 2 n 1 x0 x 0 x0 2 n 1 x1 x 1 x1 D = = ∏ (xi− x j ) ≠ 0 0 ≤ i < j ≤ n 2 n 1 xn x n x n nên hệ tồn tại duy nhất nghiệm. ¾ Nhận xét: Phương pháp này có 2 nhược điểm - Số phép tính cần thực hiện cỡ n (logn)2, chủ yếu là phép nhân. 32
- Nếu n lớn thì bậc của đa thức nội suy là quá cỡ! 3.1.3 Sơ đồ Horner tính giá trị của đa thức ¾ Bài toán: Cho đa thức (1), cần tính giá trị của đa thức tại x = c n-1 n Pn(c) = a0 + a1c + + an-1c + anc (3.4) ¾ Giải: Cách tính Pn(c) tiết kiệm nhất về số phép tính là viết lại (3.4) dưới dạng Pn(c) = ( (((anc + an-1)c + an-2)c + an-3)c + + a1)c + a0 (3.5) Vậy để tính Pn(c) ta chỉ cần tính lần lượt các số sau theo công thức truy hồi bn = an bn-1 = an-1 +bnc b0 = a0 + b1c = Pn(c) hay lập thành sơ đồ Horner sau đây an an-1 an-2 . a0 + bnc bn-1c b1c *c bn bn-1 bn-2 b0 = Pn(c) ¾ Ví dụ: 3 2 Tính giá trị của P3(x) = 3x + 2x – 5x + 7 tại x = 3 Giải: Lập sơ đồ Horner 3 2 -5 7 + 9 33 84 *3 3 11 28 91 = P3(3) 3.2 ĐA THỨC NỘI SUY LAGRANGE 3.2.1 hLập công t ức Với các mốc bất kỳ x0 < x1 < < xk < xk+1 < < xN ta có các giá trị quan sát tương ứng yk=f(xk)., k=0, , N Tìm đa thức nội suy LN(x) bậc ≤ N sao cho LN(xk) = yk, k=0,1, , N Đa thức nội suy Lagrange được cho như sau N ∏k (x) f(x)≈ Ln (x) = ∑ yk (3.6) k= 0 ∏k ()x k trong đó N ∏∏k ()x=() x − x j (3.7) j= 0, j ≠ k Kiểm tra điều kiện đường cong đa thức đi qua các điểm cho trước N ∏ k ()x i Ln (x i ) = ∑ yk ≡ yi ∀ i = 0, N k= 0 ∏ k ()x k Chẳng hạn với N=3 ta có khai triển (3.6) như sau 33
()()x− x1 x − x2 ( x− x3 ) f (x) ≈ y0 ()x0− x 1 x 0 −x ()() 2 x 0 − x 3 ()()x− x0 x− x2 ( x− x3 ) + y1 ()x1− x 0 x 1 − x ()() 2 x 1 − x 3 ()()x− x0 x− x1 ( x− x3 ) + y2 ()x2− x 0 x 2 − x ()() 1 x 2 − x 3 ()()x− x0 x− x 1 ( x− x 2 ) + y3 . ()x3− x 0 x 3 − x ()() 1 x 3 − x 2 Ta có thể rút ra cách đơn giản như sau để triển khai (3.6): • Vế phải có N+1 số hạng • Trong số hạng thứ k, (k = 0, ,N ) có 2 thừa số: Thừa số thứ hai là yk Thừa số thứ nhất là phân số mà tử số có N thừa số không có xk, mẫu số cũng có N thừa số như tử số khi thay x bởi xk. 3.2.2 Ví dụ: Tìm giá trị gần đúng của f(2,6) từ bảng số liệu xk –4 –1 2 5 8 11 yk –1 3 5 8 16 25 Áp dụng công thức (3.6) với N = 5, x0 = –4, , x5 = 11, x = 2,6: ()2,,,,, 6+ 1 2 6 − ()() 2 2 6 − 5( 2 6 − 8)( 2 6 − 11) f(,) 2 6 ≈ ()− 1 ( )−4 + () 1 − 4 − 2 − 4 − 5() − 4 − 8 − 4 − ()() 11 ()()2,,,,, 6+ 4 2 6 − 2( 2 6 − 5)( 2 6 − 8)( 2 6 − 11) + 3 ( )−1 + () 4 − 1 − 2 − 1 − 5 () − 1 − 8 − 1 − 11 ()() ()2,,,,, 6+ 4 2 6 + ()() 1 2 6 − 5( 2 6 − 8)( 2 6 − 11) + 5 ( )2+ () 4 2 + 1 2 − 5 () 2 − 8 2 − 11 ()() ()2,,,,, 6+ 4 2 6 + ()() 1 2 6 − 2( 2 6 − 8)( 2 6 − 11) + 8 ( )5+ ()4 5 + 1 5 − 2 () 5 − 8 5 − 11()() ()2,,,,, 6+ 4 2 6 + ()() 1 2 6 − 2( 2 6 − 5)( 2 6 − 11) + 16 ( )8+ ()4 8 + 1 8 − 2 () 8 − 5 8 − 11()() ()2,6+ 4 2,6 + ()() 1 2,6 − 2( 2,6 − 5)( 2,6 − 8) + 25 () 11+ () 4 11 + 1 11 − 2 () 11 − 5 11 − 8 ()() = 5,3494 3.2.3 Sai số: Người ta đã chứng minh rằng nếu hàm f(x) khả vi liên tục đến cấp N+1 trên đoạn [a,b] chứa tất cả các mốc nội suy xk, k = 0, , N thì sai số của nội suy Lagrange là f (N+ 1)(ξ) N R( x )= | f ( x ) − L ( x ) | = ()x− x, x <ξ < x (3.8) n N ()N+ 1 ! ∏ k 0 N k= 0 34
3.2.4 Sơ đồ khối và chương trình Bắt đầu Nhập n, x, xk, fk, k=0, 1, , N f = 0 k=0,1, ,N L = 1 i=0,1, , N i≠k L*=(x-xi) / (xk-xi) f += L * yk In f Dừng printf("DOC DIEM CAN TINH \n " ); scanf("%f" , &x) printf("TINH y = f(%6.2f) \n",x); y = 0 ; for (k = 0 ; k <= n ; k++) { L = 1; for (i = 0 ; i <= n ; i++) if (i != k) L *= (x- x0[i]) / (x0[k] – x0[i]) ; y += L*y0[k]; } Nhận xét: - Thuật toán có 2N(N+1) phép nhân và chia 35
- Khi các mốc nội suy cách đều, tức là xk+1 = xk + h = = x0 + kh, h là hằng số, đặt u = (x – x0) / h với u = 0, 1, , N thì công thức nội suy Lagrange (3.6) có dạng N (− 1)N− k u(u − 1) (u − k + 1)(u − k − 1) (u − N) f (x)≈ LN (x) = LN (x 0 + uh) = ∑ fk (3.9) k= 0 k!(N− k)! - Nhược điểm: Khi tăng hoặc giảm số mốc nội suy, ta phải tính lại tất cả các số hạng của đa thức mới chứ không kế thừa được các kết quả cũ. 3.3 ĐA THỨC NỘI SUY NEWTON VỚI BƯỚC CÁCH ĐỀU 3.3.1 B ảng sai phân hữu hạn Với các mốc quan sát xk = x0 + k.h, h là hằng số, ta có các gía trị quan sát tương ứng fk = f(xk), k = –m, , n . Sau đây ta sẽ trình bày khái niệm về sai phân hữu hạn. ¾ Sai phânhữu hạn Sai phân hữu hạn bậc một ∆fk = fk+1 – fk. 2 Sai phân hữu hạn bậc hai ∆ fk =∆fk+1 – ∆fk. 3 2 2 Sai phân hữu hạn bậc ba ∆ fk =∆ fk+1 – ∆ fk Bảng sai phân hữu hạn 2 3 4 xk = x0+k.h fk ∆fk ∆ fk ∆ fk ∆ fk x-1 f -1 ∆f -1 2 x 0 f 0 ∆ f -1 3 ∆f 0 ∆ f -1 2 4 x 1 f 1 ∆ f 0 ∆ f - 1 3 ∆f 1 ∆ f 0 2 4 x 2 f 2 ∆ f 1 ∆ f 0 3 ∆f 2 ∆ f 1 2 x3 f3 ∆ f 2 ∆f 3 x 4 f 4 Chú ý: 36
* Để thuận tiện cho cho việc tra cứu các số liệu, các giá trị chuẩn đặt trong khung * Giá trị cột thứ 3 trở đi được tính bằng cách lấy giá trị dưới trừ giá trị trên tương ứng của cột trước. ¾ Ví dụ Bảng sai phân hữu hạn cho hàm f(x) = x3 với các mốc quan sát x = 0 + 2k, k = –2, , 4 2 3 4 5 x k = 0 +2k fk ∆fk ∆ fk ∆ fk ∆ fk ∆ fk x - 2 = –4 –64 56 x - 1 = –2 –8 –48 8 48 x 0 = 0 0 0 0 8 48 0 x 1 = 2 8 48 0 56 48 x 2 = 4 64 96 0 152 48 x3 = 6 216 144 296 x 4 = 8 512 2 Từ bảng này ta có thể tìm f–1= –8, ∆f–2 = 56, ∆ f2 = 144, 3.3.2 Đa thức nội suy Newton tiến ¾ Bài toán Với các mốc quan sát cách đều xk = x0 + k.h, k = 0, ,n, h là hằng số và giá trị quan sát tương ứng fk = f(xk) ta cần tìm đa thức nội suy dưới dạng Pn(x) = a0 + a1(x – x0) + a2(x – x0) (x – x1) + + an(x – x0) (x – xn-1) (3.10) ¾ Công thức nội suy Newton tiến i i Cho x = x0, a0 = f0; x = x1, a1 = ∆f0 / h; , x = xi, ai = ∆ f0 / (i! h ) Đổi biến u = (x – x0) / h => x = x0 + uh Δf Δ2f Δ3f Δ n f f(x+ uh) ≈ f +0 u + 0 u(u− 1) + 0 u(u− 1)(u − 2) + + 0 u(u− 1) (u − n + 1) (3.11) 0 0 1! 2! 3! n! Nhận xét: Công thức này cho xấp xỉ tốt nhất của f(x) khi x ở đầu bảng, vì vậy người ta thường áp dụng với x0 ở đầu bảng. ¾ Sai số Người ta đã chứng minh rằng 37
hn+ 1 f () n + 1 ()ξ R() x = u( u− 1 ) ( u − n ), x xn-1 > > x0, xk = x0 + kh, h là hằng số, các giá trị quan sát tương ứng fk = f(xk), k = 0, , n ¾ Lập công thức Đặt x = xn + uh, tìm đa thức nội suy dưới dạng Pn(x) = a0 + a1(x – xn) + a2(x – xn)(x – xn-1) + + an(x – xn) (x – x1) (3.13) Cho x = xn ta được a0 = fn; x = xn-1 ta được a1 = ∆fn-1 / h; 38
i i x = xi ta được ai = ∆ fn-I / (i! h ) Đổi biến u = (x – xn) / h ta được 2 3 Δfn− 1 Δ f n− 2 Δ fn− 3 f(xn + uh) ≈ fn +u + u(u+ 1) + u( u+ 1 )( u + 2 ) + 1! 2! 3! Δnf + o u(u+ 1) (u + n − 1) (3.14) n! Nhận xét: Công thức này cho xấp xỉ tốt nhất của f(x) khi x ở cuối bảng dữ liệu, vì vậy người ta thường áp dụng với xn ở cuối bảng. ¾ Sai số: Người ta đã chứng minh rằng hn+ 1 f ( n + 1) ()ξ Rn () x = u(u+ 1) (u + n), x 1-n <ξ < x n (3.15) ()n+ 1 ! ¾ Ví dụ: Từ bảng số liệu xk –2 –1 0 1 2 3 fk 16 1 0 1 16 81 hãy tính gần đúng giá trị f( 2,2 ). Ở đây ta thấy rằng x k là các mốc cách đều với h = 1 và x = 2,2 cuối bảng giá tri x k nên ta chọn xn = 3 u = ( 2,2 – 3 ) /1 = –0,8. BẢNG SAI PHÂN HỮU HẠN LÙI TƯƠNG ỨNG 2 3 4 5 xk = x0 + kh fk Δfk Δ fk Δ fk Δ fk Δ fk x 0 = –2 16 –15 x 1 = 1 1 14 –1 –12 x 2 = 0 0 2 24 1 12 0 x 3 = 1 1 14 24 15 36 x 4 = 2 16 50 65 x 5 = 3 81 Theo công thức (3.14) ta tính được f(2,2) ≈ 81 + 65.(–0,8/1! + 50.( –0,8).( –0,8 + 1) / 2! +36.(–0,8).( –0,8 + 1)(–0,8 + 2) / 3! + 24.(–0,8).(–0,8 – 1)( –0,8 + 2)( –0,8 + 3) / 4! + 0.( –0,8).( –0,8 + 1)( –0,8 + 2)( –0,5 + 3)( –0,8 + 4) / 5! = 23,4256 39
Chú ý: - Công thức nội suy Newton tiến (lùi) chỉ là cách viết khác của công thức Lagrange với mốc cách đều (3.9) - Công thức nội suy Newton dựa trên các hệ số thứ k là các sai phân bậc k của k f(x), các sai phân này có tính chất đều (đến bậc thứ k ≤ n nào đó thì Δ fi = 0 với mọi i), do đó bậc của đa thức nội suy sẽ thấp đi - Dùng công thức nội suy Newton tiến (lùi) không phải tính lại từ đầu nếu thêm mốc nội suy mới. ¾ Sơ đồ khối và chương trình cho công thức nội suy Newton tiến Bắt đầu Nhập n, x, x0, h, yk k=0, 1, , n u= (x-x0)/h i = 0, 1. , n-1 Di = (yi+1 – yi) f = y0 +u*D0 t=u k=2,3, ,n i =0, 1, , n-k Di=(Di+1 -Di) t *= (u – k +1) f += t * D0/k! In f Dừng u = (x1 – x) / h; // x la moc noi suy dau tien for (i = 0; i <= n – 1; i++) D[i] = y[i+1] – y[i]; f = y[0] + u * D[0]; t = u; for (k = 2; k <= n; k++) { 40
for (i = 0 ; i <= n – k ; i++) D[i] = D[i+1] – D[i]; printf("\nSai phan cap %d = %8.6f",k,D[0]); t *= (u – k + 1); f += t * D[0] / GiaiThua(k); } printf("\n Noi suy Newton tien f(%8.6f) = %8.6f",x1,f); 3.3.4 Công thức nội suy Newton với mốc quan sát bất kỳ ¾ Công thức: Với các mốc quan sát bất kỳ x0 < x1 < <xk < xk+1 < < xn và các giá trị quan sát fk=f(xk), k=0, , n tương ứng ta có công thức sau đây f(x) ≈ D00 + D11(x–x0) + D22(x–x0)(x–x1) + + Dnn(x–x0)(x–x1) (x–xn-1) (3.16) trong đó các tỉ sai phân Dkj được định nghĩa như sau Dk0 = fk , k = 0, ,n (3.17) DD− D = k j− 1 k − 1j − 1 , j= 1,n; k = j,n (3.18) kj x− x k k− j ¾ Sai số f ()n+ 1 ()ξ n Rn (x) = ∏()x− xk , x0 < ξ < x n (3.19) ()n+ 1 ! k= 0 BẢNG TỶ SAI PHÂN xk Dk0 = fk Dk1 Dk2 Dk3 Dk4 x 0 D00 = f0 x 1 D10 = f1 D11 x 2 D20 = f2 D21 D22 x 3 D30 = f3 D31 D32 D33 x 4 D40 = f4 D41 D42 D43 D44 ¾ Ví dụ: Tính gần đúng f(2,6) từ bảng số liệu với 4 chữ số sau dấu phẩy xk –4 –1 2 5 8 11 fk –1 3 5 8 16 25 Theo công thức (3.17):ta có D00 = -1 , D10 = 3 , D20 = 5 , D30 = 8 , D40 = 16 , D50 = 25 Theo công thức (3.18) : 41
D11= (D10 – D00) / (x1 – x0) = (3 + 1) / (–1 + 4) = 1,3333 D21= (D20 – D10) / (x2 – x1) = (5 – 3) / (2 + 1) = 0,6667 D31= (D30 – D20) / (x3 – x2) = (8 – 5) / (5 – 2) = 1 D41= (D40 – D30) / (x4 – x3) = (16 – 8) / (8 – 5) = 2,6667 D51= (D50 – D40) / (x5 – x4) = (25 – 16) / (11 – 8) = 3 D22= (D21 – D11) / (x2 – x0) = (0,6667 – 1,3333) / (2 + 4) = –0,1111 D32= (D31 – D21) / (x3 – x1) = (1 – 0,6667) / (5 + 1) = 0,0556 D42= (D41 – D31) / (x4 – x2) = (2,6667 – 1) / (8 – 2) = 0,2778 D52= (D51 – D41) / (x5 – x3) = (3 – 2,6667) / (11–5) = 0,0556 D33= (D32 – D22) / (x3 – x0) = (0,0556 + 0,1111) / (5 + 4) = 0,0185 D43= (D42 – D32) / (x4 – x1) = (0,2778 – 0,0556) / (8 + 1) = 0,0247 D53= (D52 – D42) / (x5 – x2) = (0,0556 – 0,2778) / (11 – 2) = –0,0247 D44= (D43 – D33) / (x4 – x0) = (0,0247 – 0,0185) / (8 + 4) = 0,0005 D54= (D53 – D43) / (x5 – x1) = (–0,0247 – 0,0247) / (11 + 1) = –0,0041 D55= (D54 – D44) / (x5 – x0) = (–0,0041 – 0,0005) / (11 + 4) = –0,0003 BẢNG TỶ SAI PHÂN xk Dk0 = fk Dk1 Dk2 Dk3 Dk4 Dk5 –4 –1 –1 3 1,3333 2 5 0,6667 –0,1111 5 8 1 0,0556 0,0185 8 16 2,6667 0,2778 0,0247 0,0005 11 25 3 0,0556 –0,0247 –0,0041 –0,0003 Theo công thức (3.16): f(x) ≈ –1 + 1,3333(x + 4) – 0,1111(x + 4)(x + 1) + 0,0185(x + 4)(x + 1)(x – 2) + 0,0005(x + 4)(x + 1)(x – 2)(x – 5) – 0,0003(x + 4)(x + 1)(x – 2)(x – 5)(x – 8) Thay x = 2,6 ta được f(2,6) ≈ 5,3859 42
3.4 NỘI SUY SPLINE 3.4.1 Đặt vấn đề Các đa thức bậc thấp là những đường cong đơn giản mà người ta muốn dùng cho việc nội suy. Ý tưởng của phương pháp nội suy Spline (hàm ghép trơn) là thay vì tìm một dường cong đa thức đi qua mọi điểm nội suy, người ta tìm những đường cong bậc ba khác nhau nối các điểm kề nhau sao cho chúng trơn ở các điểm nối Mi(xi,yi) đã cho, i = 0, 1, , n. y 5 4 6, 4 3 1, 3 4, 3 8, 3 12, 3 2 1 2, 1 10, 1 0 x 0 5 10 15 3.4.2 Bài toán Giả sử có hàm f(x) xác định trên đoạn [x0, xn] và biết được giá tri của nó tại các điểm x0 < x1 < < xn là y0, y1, , yn tương ứng. 3 2 Hãy tìm n đường ghép trơn bậc ba Si(x) = aix + bix + cix + di (3.20) xác định trên [xi-1, xi], i = 1, 2 , n sao cho ⎧ S (x )= y , i = 1, n (3.21) ⎪ i i-1 i− 1 ⎪Si (x i )= y i , i = 1, n (3.22) ⎪ ' ' ⎪Si (xi )= S i+ 1 (xi ), i= 1, n -1 (3.23) ⎨ '' '' ⎪Si (xi )= S i+ 1 (xi ), i= 1, n -1 (3.24) ⎪ '' S1 (x0 )= 0 ⎪ '' ⎩⎪Sn (xn )= 0 3.4.3 Xây dựng công thức Ta thấy hệ phương trình điều kiện có đủ 4n phương trình để tìm 4n ẩn ai, bi, ci, di (i = 1, , n) Đặt '' '' '' Si (xi )= p i , i = 1, n -1, p0= S 1 (x0 ) = 0, pn = S n (xn ) = 0 (3.25) '' Do Si (x) là đa thức bậc nhất trên [xi-1, xi] nên có thể tìm nó dưới dạng '' Si (x) = α i(x i − x) +βi (x − xi-1 ) . Cho x = xi-1 và x = xi ta có 43
pi-1 pi-1 pi pi αi = = , βi = = , hi = x i − x i− 1 xi− x i-1 hi xi− x i-1 hi Khi đó lấy tích phân xác định 2 lần hàm Si’’(x) trên [xi-1, xi] ta được pi-1 3 pi 3 ~ ~ Si (x) =(xi − x) +(x − xi-1 ) +αi(x i − x) +βi(x − x i-1 ) (3.26) 6hi 6h i Lần lượt cho x = xi và x = xi-1 ta tính được 2 2 ~ 1 pi h i ~ 1 pi-1 h i βi =(yi −), αi =(yi-1 − ) hi 6 hi 6 Đổi biến x - xi-1 t = hi thay vào (3.26) ta có h2 {(t 3 − t)p + [(1 − t)3 − (1 − t)]p } S (t)= ty + (1 − t)y + i i i-1 (3.27) i i i-1 6 Khi x [xi, xi+1] thì t [0, 1], i = 1, , n-1 và các điều kiện từ (3.21) đến (3.24) trở thành Si(1) = Si+1(0) = yi, i = 1, , n-1 S'i (1)= S'i+ 1 (0), i= 1, , n -1 '' '' Si (1)= Si+ 1 (0) = pi , i = 1, , n -1 Các pi được xác định từ điều kiện (3.25) bằng cách lấy đạo hàm cấp một theo x của (3.20) viết theo t: 2 2 1 h {(3t −1)p −[3(1 − t) −1]p } y− y [S' (t)] = [S' (t)] t' = [S' (t)] = z + i i i -1 , z = i i-1 i x i t x i t h i 6 i i h i Do S'i (1)= S'i+ 1 (0) ∀ i = 1,n -1 ta có n-1 phương trình với n-1 ẩn p1, p2, , pn-1 hi pi -1 + 2 (hi + hi+1) pi + hi+1 pi+1 = 6 (zi+1 – zi), i = 1, , n –1 (3.28) với ma trận hệ số là ⎛ 2(h+ h ) h 0 0 0 0 0 w ⎞ ⎜ 1 2 2 1 ⎟ ⎜ h2 2(h2+ h 3 ) h3 0 0 0 0 w 2 ⎟ ⎜ 0 h 2(h+ h ) h 0 0 0 w ⎟ ⎜ 3 3 4 4 3 ⎟, ⎜ ⎟ ⎜ ⎟ ⎜ 0 0 0 0 0 hn− 2 2(h n− 2+ h n − 1 ) h n-1 w n− 2 ⎟ ⎜ ⎟ ⎝ 0 0 0 0 0 0 h n-1 2(hn-1+ h n ) w n− 1 ⎠ wi = 6(z i+ 1 − z i ) 44
Đây là hệ n – 1 hệ phương trình tuyến tính n – 1 ẩn dạng 3 đường chéo, đối xứng qua đường chéo chính và đường chéo chính là trội, nên luôn tồn tại duy nhất nghiệm. 3.4.4 Các bước giải bài toán nội suy Spline bậc ba Nhập ma trận hệ số của hệ phương trình tuyến tính Giải hệ phương trình tuyến tính tìm các giá trị pi , i = 1, , n – 1. Kiểm tra giá trị x cần nội suy nằm trên khoảng [xi–1, xi] nào rồi tính x - x t = 1-i h i Sử dụng công thức (3.26) tính được Si(t) 3.4.5 Ví dụ i xi hi yi yi – yi-1 zi = (yi – yi-1)/hi wi = 6(zi+1–zi) Cho bảng các giá trị 0 1 2 của hàm f(x). Hãy tính 1 2 1 1,5 –0,5 –0,5 2,25 2 4 2 1,25 –0,25 –0,125 0,45 các hệ số pi của đa thức ghép trơn Spline bậc ba 3 5 1 1,2 –0,05 –0,05 0,15 trên các đoạn đã cho 4 8 3 1,125 –0,075 –0,025 0,075 rồi tính gần đúng 5 10 2 1,1 –0,025 –0,0125 f(6,46). Giải: Các hệ số p1, p2, p3, p4 là nghiệm của hệ 3 đường chéo sau ⎛6 2 0 0 2,250 ⎞ ⎜ ⎟ ⎜2 6 1 0 0,450 ⎟ ⎜0 1 8 3 0,150 ⎟ ⎜ ⎟ ⎜ ⎟ ⎝0 0 3 10 0,075 ⎠ Giải hệ trên ta tìm được p1= 0,39541; p2 = –0,06123; p3 = 0,02658; p4 = –0,00048 (p0= p5 = 0) Do x = 6,46 thuộc vào đoạn thứ i = 4, cụ thể là x ∈ [5; 8] nên t =(6,46 – 5) / (8 – 5) ≈ 0,48667 Thay giá trị của t và các giá trị pi-1 = p3 = 0,02658 và pi = p4 = –0,00048 vào (3.27) ta tính được f(6,46) ≈ 1,14869. So với giá trị đúng của f(x) = 1 + 1/x tại x = 6,46 là 1,15480 chỉ sai khác 0,6% 3.4.6 Chương trình tính float Spline (int n, mang x, mang y, float v) { mang h, h1, d, w, p; int i; float t, s; for (i = 1; i <= n ; i++) { h[i] = x[i] - x[i-1]; 45
printf(" h[%d] = %10.6f \n",i,h[i]); } for (i = 1; i x[i]) i++; t=(v–x[i–1])/h[i]; s=t*y[i]+(1–t)*y[i–1]+h[i]*h[i]*(Tich(t)*p[i]+Tich(1–t)*p[i–1])/6; return s; } Hàm Tich(t) được định nghĩa như sau: float Tich(float t) { return t*t*t - t; } Nhận xét: thuật toán có độ phức tạp thời gian tuyến tính 3.5 PHƯƠNG PHÁP BÌNH PHƯƠNG BÉ NHẤT LÀM KHỚP DỮ LIỆU 3.5.1 Đặt vấn đề: Phương pháp bình phương bé nhất thường được dùng để lập công thức thực nghiệm. Giả sử cần tìm quan hệ hàm giữa 2 đại lượng x và y dựa vào bảng số liệu thực nghiệm x x1 x2 xn y y1 y2 yn Bên cạnh đó, ta còn biết dạng của hàm phải tìm. Hàm này phụ thuộc vào một số tham số, chẳng hạn f(x) = f(c1, c2, , cm; x) Việc chọn hàm xấp xỉ tốt nhất trở thành việc chọn các tham số phù hợp nhất với các điểm thực nghiệm. Tiêu chuẩn để chọn là tổng các bình phương của sai số ở các điểm * * * thực nghiệm là bé nhất, tức là bộ tham số phù hợp nhất phải là bộ số (c1 ,c2 , ,cm ) làm cực tiểu hàm m biến n 2 E(c1 ,c 2 , ,c m )= ∑ [f(c1 ,c 2 , ,c m ; x i ) - y i ] (3.28) i= 1 Trong phạm vi giáo trình này, ta chỉ xét trường hợp dạng tham số tuyến tính tức là hàm thực nghiệm là tổ hợp tuyến tính của các hàm đơn giản hơn với các hệ số của tổ hợp chính là các tham số phải tìm f(x) = c1f1(x) + c2f2(x) + + cmfm(x) (3.29) chẳng hạn 46
f(x) = a+ bx; f(x) = a + bx + cx2; f(x) = a + b/x; f(x) = a + bcosx + csinx; f(x) = a + blnx hoặc có thể biến đổi về dạng tham số tuyến tính như f(x) = aebx; (a > 0) f(x) = axb; (a > 0, b > 0) 3.5.2 Lập công thức * * * Điểm (c1 ,c2 , , cm ) làm cực tiểu hàm E(c1,c2, ,cm) phải là điểm dừng của hàm này, tức là phải thoả mãn hệ m phương trình n ∂E ∂f(c1 ,c 2 , ,c m ;x i ) = 2∑ [f(c1 ,c 2 , cm ; x i )− y i ] =0, ∀ j = 1,m ∂c j i= 1 ∂c j Khi f(c1,c2, ,cm;xi) = c1f1(xi) + c2f2(xi) + + cmfm(xi), thay vào hệ trên ta được hệ đối xứng ⎧ c1 (f 1 ,f 1 )+ c 2 (f 1 ,f 2 )+ + cm (f 1 ,f m )= (f 1 , y) ⎪ ⎪ c1 (f 2 ,f 1 )+ c 2 (f 2 ,f 2 )+ + cm (f 2 ,f m )= (f 2 , y) ⎨ (3.30) ⎪ ⎪ ⎩c1 (f m ,f 1 )+ c 2 (f m ,f 2 )+ + cm (f m ,f m )= (fm , y) trong đó fi = (fi(x1), , fi(xn)), i = 1, , m và y = (y1, , yn), còn (fi, fj) là tích vô hướng giữa 2 véc tơ fi và fj., cụ thể n n ()fi ,f j= ∑ f i (xk )f j (x k ) , (fi , y)= ∑ fi (xk )y k k= 1 k= 1 Hệ phương trình tuyến tính(3.30) là đối xứng, giải được bằng phương pháp Gauss- Joordan. 3.5.3 Các ví dụ: Ví dụ 1: Giả sử có 11 điểm thực nghiệm x 1 2 3 4 5 6 7 8 9 10 11 y 0 0.6 1.77 1.92 3.31 3.52 4.59 5.31 5.79 7.06 7.17 Tìm đường thẳng y = f(x) = c1 + c2x xấp xỉ tốt nhất tập điểm trên theo tiêu chuẩn bình phương bé nhất. Giải: Trong trường hợp này f1(x) = 1, f2(x) = x 47
Lập bảng i x = f2 y f1 (f1,f1)(f1,f2) (f1,y) (f2,f2)(f2,y) 1 1 0,00 1 1 1 0 1 0 2 2 0,60 1 1 2 0,6 4 1,2 3 3 1,77 1 1 3 1,77 9 5,31 4 4 1,92 1 1 4 1,92 16 7,68 5 5 3,31 1 1 5 3,31 25 16,55 6 6 3,52 1 1 6 3,52 36 21,12 7 7 4,59 1 1 7 4,59 49 32,13 8 8 5,31 1 1 8 5,31 64 42,48 9 9 5,79 1 1 9 5,79 81 52,11 10 10 7,06 1 1 10 7.06 100 70,60 11 11 7,17 1 1 11 7,17 121 78,87 Tổng 11 66 41,04 506 328,05 Thay vào (3.30) ta có hệ 2 phương trình đối xứng của 2 ẩn c1 và c2 là ⎛11 66 41.04 ⎞ ⎜ ⎟ ⎜ ⎟ ⎝66 506 328.05⎠ Giải ra ta được c1 = –0,7314 và c2 = 0,7437. Vậy f(x) = –0,7314 + 0,7437x Ví dụ 2: Giả sử có 6 điểm thực nghiệm x 1 2 4 5 8 10 y 2,0 1,5 1,25 1,2 1,125 1,1 có thể làm khớp bởi đường cong f(x) = c1 + c2 / x. Hãy xác định c1, c2 bằng phương pháp bình phương bé nhất Giải: Trong trường hợp này f1(x) = 1, f2(x)=1/x. Ta tính được Lập bảng i x y f1 = 1 f2 = 1/x (f1,f1) (f1,f2) (f1,y) (f2,f2) (f2,y) 1 1 2,000 1 1,000 1 1,000 2,000 1,000 2,000 2 2 1,500 1 0,500 1 0,500 1,500 0,250 0,750 3 4 1,250 1 0,250 1 0,250 1,250 0,063 0,313 4 5 1,200 1 0,200 1 0,200 1,200 0,040 0,240 5 8 1,125 1 0,125 1 0,125 1,125 0,016 0,141 6 10 1.100 1 0.100 1 0,100 1,100 0,010 0,110 Tổng 6 2,175 8,175 1,378 3,553 Thay vào (3.30) ta có hệ 2 phương trình đối xứng của 2 ẩn c1 và c2 là 48
⎛ 6 2,175 8,175⎞ ⎜ ⎟ ⎜ ⎟ ⎝ 2,175 1,378 3,553 ⎠ Giải ra ta được c1 = 1 và c2 = 1 hay f(x) = 1 + 1/x 3.5.4 Các bước giải và chương trình - Nhập các véc tơ x = (x1, , xn) và y = (y1, , yn) - Lập ma trận Fmxn các vectơ fi = (fi(x1), ,fi(xn)), i = 1, , m - Lập ma trận các hệ số của phương trình AX = b với aij = (fi, fj); i, j = 1, , m và b = ai, m+1 = (fi, y), i = 1, , m - Giải hệ AX = b ta được bộ tham số (c1*, c2*, , cm*) tốt nhất của hàm f(c1, c2, ,cm; x) đã cho for (i = 1; i <= m; i++) for (j = 1; j <= m+1; j++) { t = 0. ; for (k = 1; k <=n; k++) t += f[i][k] * f[j][k]; a[i][j] = t; } for (i = 1; i <= m; i++) { t = 0. ; for (k = 1; k <= n; k++) t += f[i][k]*y[k]; a[i][m+1] = t; } Gauss(m, a, c); 49
3.6 BÀI TẬP CHƯƠNG 3 1> Từ bảng các logarit chúng ta nhận được các giá trị sau đây của log x tại các điểm ấn định của bảng x log x 1,0 0.0 1,5 0,17069 2,0 0,30103 3,0 0,47712 3,5 0,54407 4,0 0,60206 nội suy các giá trị sau đây: log2,5 và log1,25 sử dụng đa thức nội suy Newton. 2> Đối với bảng dữ liệu x –2 –1 0 1 2 y –4 1 1 2 10 Tính các giá trị của đa thức nội suy Newton P5(–0,5) và P5(1,5) 3> Tạo bảng sai phân cho dữ liệu dưới đây, và đánh giá bậc của đa thức nội suy Newton cần thiết để sinh ra các giá trị nội suy đúng cho số chữ số có nghĩa đã biết. X f(x) 1 1,5709 2 1,5713 3 1,5719 4 1,5727 5 1,5738 6 1,5751 7 1,5767 8 1,5785 9 1,5805 4> Đối với bảng dữ liệu dưới đây, sử dụng đa thức ghép trơn spline bậc ba S(x) để tính nội suy tại các điểm x=10 và x=35 x 5 15 30 50 70 y 1,201 1,824 2,624 2,787 2,210 5> Tìm các hệ số a*, b*, c* sao cho hàm f(x) = ax2+bx+c khớp tốt nhất đối với bảng dữ liệu cho trong bài tập 4 theo tiêu chuẩn bình phương bé nhất. 50
4 CH ƯƠNG 4: TÍNH ĐẠO HÀM VÀ TÍCH PHÂN XÁC ĐỊNH 4.1 TÍNH GẦN ĐÚNG ĐẠO HÀM 4.1.1 Xấp xỉ giá trị đạo hàm dựa vào bảng sai phân Với các mốc quan sát xk = x0 + k.h, h là hằng số, ta có các gía trị quan sát tương ứng fk = f(xk), k = –m, , n. Các giá trị của cột ∆fk trong bảng số gia hữu hạn.được i dùng để ước lượng đạo hàm bậc nhất, các giá trị của các cột ∆ fk được dùng để ước lượng các đạo hàm bậc i tương ứng của hàm f(x) tại xk. BẢNG SAI PHÂN HỮU HẠN 2 3 4 fk fk fk fk fk x-1 f –1 f –1 2 x 0 f 0 f –1 3 f 0 f –1 2 4 x 1 f 1 f 0 f –1 3 f 1 f 0 2 4 x 2 f 2 f 1 f 0 3 f 2 f 1 2 x3 f3 f 2 f 3 x 4 f 4 Dựa vào định nghĩa Δf f '(x)= lim Δx→ 0 Δx ta có 2 giá trị gần đúng sau đây với ∆x là hằng số dương: Δf Δf f '(x ) ≈ k , f '(x ) ≈ k− 1 k Δx k Δx lấy trung bình cộng của 2 công thức trên ta có công thức xấp xỉ: Δfk + Δ f k− 1 f '(xk ) ≈ (4.1) 2Δ x với sai số là: Δx 2 ε = f()3 ()ξ , x - Δ x < ξ < x + Δ x (4.2) 1 6 Các xấp xỉ đạo hàm cấp cao được tính tương tự như công thức (4.1) Ví dụ: Tính gần đúng f’(1,06) từ bảng số liệu sau đây với ∆x = 0,01 xk 1,03 1,04 1,05 1,06 1,07 1,08 1,09 1,10 51
fk 0,012837 0,017033 0,021189 0,025306 0,029384 0,033424 0,037426 0,041393 Theo công thức (4.1): Δf ()1,05+ Δ f ()1,06 0,004117 + 0,004078 f'(,) 1 06 ≈ = = 0. 40975 2.0,01 0,02 4.1.2 Xấp xỉ đạo hàm bằng công thức nội suy Ta có thể sử dụng một trong các công thức nội suy của chương III, thay f(x) bằng đa thức nội suy của nó: f(x) = P(x) + Rn(x). Khi đó ' f ’(x) = P’(x) + Rn (x) Nhận xét: Phép tính gần đúng đạo hàm như trên không chính xác bằng phép nội suy trực tiếp vì từ hệ thức f(xi) = P(xi) nói chung không suy ra được f '(xi )= P'(xi ) ¾ Trường hợp các mốc nội suy cách đều • Công thức nội suy Newton tiến: Sử dụng công thức 2 N Δf0 Δ f0 Δ f 0 f() x≈ f0 +u + u(u− 1) + + u( u− 1 ) ( u − N + 1) 1! 2! N! với u =(x-x0)/h, f0=f(x0) để tính gần đúng giá trị của đạo hàm f’(x) khi x = x0: 2 N ' ⎡Δf0 Δ f0 Δ f0 ⎤ ' f'( x 0 ) ≈⎢ u + u( u− 1 ) + + u( u− 1 ) ( u − N + 1)⎥ ux () x 0 (4.3) ⎣ 1! 2! N! ⎦ u N k− 1 k 1 (−1) Δ f0 Sau một số biến đổi ta được f '(x0 ) ≈ ∑ (4.4) h k= 1 k với sai số là N() N+ 1 N N+ 1 N+ 1 ' h f (ξ) hΔ f 0 Δ f 0 rn()() x 0= R n x 0 = ≈ = , x<ξ < x 0 (4.5) N+ 1 ()N+ 1 h N+ 1 ()N+ 1 h Ví dụ: Tính gần đúng giá trị của f '(1,06) theo công thức nội suy Newton tiến từ bảng số liệu trong ví dụ trên Ở đây ta cần tính giá trị đạo hàm tại x = 1,06 vì vậy ta có thể chọn x0 = 1,06 và lập bảng số gia hữu hạn tiến tương ứng 2 3 4 xk=1,06 + 0,01k fk Δfk Δ fk Δ fk Δ fk x0 =1,06 0,025306 0,004078 x1 =1,07 0,029384 –0,000038 0,004040 0 x2 =1,08 0,033424 –0,000038 0,000003 0,004002 0,000003 x3 =1,09 0,037426 –0,000035 0,003967 x4 =1,10 0,041393 52
Áp dụng công thức (4.4) với N = 4, h = 0,01 ta tìm được 1 ⎡Δf Δ2f Δ3f Δ4f ⎤ f '(1,06) ≈ ⎢ 0 − 0 + 0 − 0 ⎥ h ⎣⎢ 1 2 3 4 ⎦⎥ = (1/0,01).(0,004078 + 0,000038/2 + 0 - 0,000003/4) = 0,4096 • Công thức nội suy Newton lùi: cách tính hoàn toàn tương tự • Sơ đồ khối và chương trình tính đạo hàm bằng công thức nội suy Newton tiến Bắt đầu Nhập n, x0, h, yk, k=0, 1, , n dấu = (-1)0 i = 0, 1. , n-1 Di = (yi+1 – yi) dhf = dấu * D0 k=2,3, ,n i =0, 1, , n-k Di=(Di+1 - Di) dấu *= ( -1) dhf + = dau * D0 / k In dhf / h Dừng dau = (–1)*( –1); 53
for (i = 0 ; i <= n–1 ; i++) D[i] = y[i+1] – y[i]; dhf = dau * D[0]; for (k = 2 ; k <= n; k++) { for (i = 0 ; i <= n –k ; i++) D[i] = D[i+1] – D[i]; printf("\nSai phan cap %d = %8.6f",k,D[0]); dau *= (–1); dhf += dau * D[0] / k; } printf("\nDao ham tai x = %8.6f la f = %8.6f",x[0],dhf/h); Nhận xét: Thuật toán có độ phức tạp thời gian cỡ O(n2) ¾ Trường hợp các mốc nội suy không cách đều – Phương pháp Spline bậc ba Dùng Spline bậc ba là là một trong những phương pháp tiện lợi nhất hiện nay. Áp dụng công thức '' pi-1 pi Si (x) = (xi − x) +(x − xi− 1 ), h i = x i − x i− 1, ∀ i = 1,n (4.6) h i h i '' Lấy tích phân xác định hàm Si (x) trên [xi-1, xi] ta được 2 2 ' pi− 1 2 pi 2 1 pi-1 h i 1 pi h i Si (x) = −(xi − x) +(x − xi− 1 ) − (yi− 1 − ) + (yi − ) (4.7) 2h i 2h i h i 6 h i 6 ta tính được f ’(x) ≈ Si’(x), f ’’(x) ≈ Si’’(x) khi xi-1 ≤ x ≤ xi Ví dụ: Cho bảng các giá trị của hàm f(x): x 1.0 2.0 4.0 5.0 8.0 10.0 Hãy tính các hệ số p của Spline bậc ba i y 2.0 1.5 1.25 1.2 1.125 1.1 trên các đoạn đã cho rồi tính gần đúng f ’(6,46) Giải: Các hệ số pi là nghiệm của hệ 3 đường chéo sau ⎛6 2 0 0 2,250 ⎞ ⎜ ⎟ ⎜2 6 1 0 0,450 ⎟ ⎜0 1 8 3 0,150 ⎟ ⎜ ⎟ ⎜ ⎟ ⎝0 0 3 10 0,075 ⎠ Giải hệ trên ta tìm được p1=0,39541; p2=–0,06123; p3=0,02658; p4 –0,00047(p0=p5 =0) Do x = 6,46 thuộc vào đoạn thứ i = 4, cụ thể là x ∈[5, 8] nên pi–1 = p3 = 0,02658 và pi = p4 = –0,00047. Thay vào (4.7) ta có: 54
' ' f (6,46) ≈ S4 (6,46) = 2 2 p3 2 p4 2 1 p3 h 4 1 p4 h 4 = −(x4 − x) +(x − x3 ) − (y3 - ) + (y4 - ) 2h4 2h4 h4 6 h4 6 0,02658 - 0,00047 1 0,02658. 9 1 - 0,00047. 9 = − (8− 6,46)2 + (6,46− 5)2 − (1,2 - ) + (1,125 - ) 2.3 2.3 3 6 3 6 − 0,02658 0,00047 1 0,23922 1 0,00423 = 2,3716− 2,1316− (1,2 − ) + (1,125 + ) 6 6 3 6 3 6 = − 0,01051 − 0,00017 − 0,38671 + 0,37524 = − 0,02215 Nhận xét: - So với kết quả lấy đạo hàm của f(x) = 1 + 1/x, f ’( 6,46) = – 0,02396 chỉ sai khác 0,2%. - Chương trình tính Spline bậc ba chỉ cần bổ sung thêm công thức (4.7) và (4.6) để tính gần đúng đạo hàm cấp một và cấp hai của f(x) tại x bất kì trên đoạn [x0, xn] 4.2 TÍNH GẦN ĐÚNG TÍCH PHÂN XÁC ĐỊNH b Trong nhiều trường hợp ta cần tính giá trị của tích phân xác định I = ∫ f(x)dx a nhưng không thể biểu diễn nguyên hàm dưới dạng hàm sơ cấp hoặc dù có thể tìm nguyên hàm dưới dạng hàm sơ cấp nhưng đó là biểu thức phức tạp. Khi đó ta phải áp dụng các cách tính gần đúng giá trị của tích phân xác định mà các công thức ban đầu thường được thiêt lập từ ý nghĩa hình học của tích phân xác định. 4.2.1 Lập công thức chung sử dụng đa thức nội suy Newton tiến Cho f(x) khả tích trên [a, b], hãy tính gần đúng tích phân I ¾ Lập công thức: Chia [a, b] thành n đoạn nhỏ đều nhau bởi các điểm chia a = x0 x = x0 + uh: 55
f 2 3 Δ 0 Δ f0 Δ f0 f (x)= f (x0 + uh) ≈ f0 + u + u() u− 1 + u( u− 1 )( u − 2 ) 1! 2! 3! n Δ f0 +K +u( u − 1 ) ( u − n + 1 ) n! b xn n 2 3 ⎡ Δf0 Δ f0 Δ f0 f()() x dx= f x dx ≈h ⎢f0 +u + u() u− 1 + u(u− 1)(u − 2) ∫ ∫ ∫ ⎣ 1! 2! 3! a x0 0 Δnf ⎤ + + 0 u(u− 1) (u − n + 1)⎥ du n! ⎦ 2 2 n 3 n ⎡ Δf n Δ f0 Δ f0 =h⎢ nf +0 + u() u− 1 du + u( u− 1 )( u − 2 ) du 0 2.! 1 2! ∫ 3! ∫ ⎣⎢ 0 0 Δnf n ⎤ + + 0 u( u− 1 ) ( u − n + 1) du⎥ (4.8) n! ∫ 0 ⎦⎥ ¾ Công thức hình thang - Khi n=1: Thay vào (4.8) ta có công thức nội suy Newton tuyến tính b ⎡ Δf0 ⎤ f0+ f 1 ∫ f(x)dx ≈ h⎢ f0 + ⎥ = h (4.9) a ⎣ 2.! 1 ⎦ 2 Về mặt hình học, vế phải của (4.9) là diện tích hình thang vuông có đường cao h và hai đáy f0, f1 - Khi n tùy ý: ta tách tích phân ban đầu thành tổng của n tích phân trên các đoạn con. Mỗi tích phân này được tính gần đúng bằng công thức hình thang (4.9) b n xi n f+ f h f(x)dx = f(x)dx≈ h i-1 i =[f + 2f + 2f + f ] ∫ ∑ ∫ ∑ 0 1 n− 1 n (4.10) i== 1 i 1 2 2 a xi-1 Công thức (4.10) là công thức hình thang tổng quát i xi yi Sai số: Nếu hàm f(x) có đạo hàm cấp hai liên tục trên đoạn [a, b] 0 0 1 thì sai số của công thức (4.10) là 1 0,1 0.90909 h2 f''()ξ 2 0,2 0,83333 ε = ()b− a , a <ξ < b (4.11) 3 0,3 0,76923 12 Ví dụ: 4 0,4 0,71429 Dùng công thức hình thang tổng quát, tính gần đúng 5 0,5 0,66667 1 dx 6 0,6 0,62500 I = ∫ 7 0,7 0,58824 01+ x 8 0,8 0,55556 với các điểm chia xk = 0,1xk, k = 0, , 10 và tìm sai số. 9 0,9 0,52632 Giải: n = 10, a = 0, b = 1, h = 1/10 = 0,1 các giá trị được tính 10 1 0,5 theo bảng bên: Theo công thức (4.10) ta có: 56
1 dx I = ∫ 0 1+ x ⎡1+ 0, 5 ⎤ ≈ ⎢ +0,,,,,,,,,, 90909 + 0 83333 + 0 76923 + 0 71429 + 0 66667 + 0 625 + 0 58824 + 0 55556 + 0 52632⎥ 01 ⎣ 2 ⎦ = 0,69377 Sai số ε được đánh giá từ 2 1,0.2 2 )x(''f = ⇒ max | f'' (x) |= 2 ⇒ ε ≤ =0,00167 ≈ 0,002 (1+ x)3 0 #include #include double f(float x) { double y; y = 1/(1+x) ; return y ; } void main() { int i , n; double h, a, b, s; clrscr(); printf("TICH PHAN BANG PHUONG PHAP HINH THANG CUA HAM y=1/(1+x)\n"); printf("Nhap cac can a va b: "); scanf("%lf %lf",&a,&b); printf("Nhap so cac khoang chia n: "); scanf("%d",&n); s = (f(a) + f(b)) / 2 ; h = (b-a) / n; for (i = 1 ; i <= n–1 ;i++) s += f(a+i*h) ; printf("TICH PHAN I= %10.3lf tren [%lf,%lf] voi n=%d\n",s*h,a,b,n); getch(); } Nhận xét: Thuật toán có độ phức tạp thời gian tuyến tính O(n) ¾ Công thức Simpson (công thức parabol) - Khi n = 2: Thay vào (4.8) ta có công thức nội suy Newton với ba mốc cách đều: b 2 2 2 Δf0 n Δ f0 h I = f() x dx ≈ h[ nf0 + + ∫ u(u− 1)du] =(f + 4f + f ) ∫ 2 1!. 0 1 2 a !2 0 3 (4.12) 57
Về mặt hình học, vế phải của (4.12) là diện tích của hình thang cong mà cạnh cong là cung parabol đi qua ba điểm A(x0, y0), B(x1, y1) và C(x2, y2) của đường cong y = f(x) nên (4.12) vừa được gọi là công thức Simpson vừa được gọi là công thức parabol. - Khi n = 2k với k tùy ý: ta tách tích phân ban đầu thành tổng của n tích phân trên các đoạn con chứa ba mốc nội suy liên tiếp( x0, x1, x2), ,(x2i-2, x2i-1, x2i), , (x2n-2, x2n-1, x2n). Mỗi tích phân này được tính gần đúng bằng công thức Simpson (4.12) x b k 2i h k ∫ f(x)dx = ∑ ∫ f(x)dx ≈ ∑ (f2i− 2+ 4f 2i− 1 + f 2i ) a i== 1 3 i 1 x 2i-2 (4.13) h =[(f + f ) + 4(f + f + + f ) + 2(f + f + + f )] 3 0 2k 1 3 2k− 1 2 4 2k− 2 Công thức (4.13) là công thức Simpson tổng quát i xi y2i-1 y2i Sai số: Nếu hàm f(x) có đạo hàm cấp bốn liên tục trên 0 0 đoạn [a, b] thì sai số của công thức (4.13) là 1 0,1 0.90909 4 (4) h f ()ξ b− a 2 0,2 0,83333 ε = ()b− a , a < ξ < b với h = 180 2k 3 0,3 0,76923 (4.14) 4 0,4 0,71429 5 0,5 0,66667 Ví dụ: Dùng công thức Simpson tổng quát, tính 6 0,6 0,62500 1 dx 7 0,7 0,58824 I = ∫ 01+ x 8 0,8 0,55556 với các điểm chia xi = 0,1x i, i = 0, , 10 và tìm sai số. 9 0,9 0,52632 10 1 Giải: với n = 2k = 10, a = 0, b = 1, h = 1/10 = 0,1 ∑ 3,45995 2,72818 các giá trị được tính theo bảng bên. Theo công thức (4.13) ta có: 1 dx 1,0 I = ∫ ≈[1 + 0,0467 + 4.3,45955+ 2.2,272818]= 0,69315 0 1+ x 3 Sai số ε được đánh giá từ (4.14) 24 24.0,14 f(4) (x) = ⇒ max | f(4) (x) |= 24 ⇒ε ≤ = 0,000013 (1+ x)5 0< x < 1 180 Chương trình tính: lập trình tương tự như đối với công thức hình thang Nhận xét: 58
- Từ (4.11) và (4.14) ta thấy công thức Simpson có độ chính xác cao hơn công thức hình thang với cùng số điểm chia n - Việc đánh giá sai số của (4.10) và (4.13) phụ thuộc vào f(x) và gặp khó khăn khi chỉ biết f(x) dưới dạng bảng 4.2.2 Quy tắc làm tăng độ chính xác của việc tính tích phân Đặt vấn đề: Cần tìm số điểm chia [a,b] để đảm bảo độ chính xác cho trước hoặc với số lần lặp đủ lớn. Muốn vậy ta dựa vào quy tắc như sau Tiến hành tính gần đúng các tích phân b Sj ≈ ∫ f (x)dx a j với các bước chia lần lượt là hj = (b-a) / 2 , j = 1,2, bằng thuật toán lặp với điều kiện dừng ở bước lặp thứ j là • với ε > 0 cho trước |S -S |<ε (4.15) j j-1 • hoặc j = M (số lần lặp cho trước) Như vậy công thức hình thang được tính lặp h j STj =[f0 + 2f 1 + + 2fj + fj ] (4.16) 2 2 −1 2 và công thức Simpson được tính lặp h j SSj =[f0 + fj + 4(f 1 + f 3 + + fj )+ 2(f2 + f 4 + + fj )] (4.17) 3 2 2 − 1 2 − 2 Ví dụ : Với ε = 0,001 và số lần lặp tối đa M = 10 hãy tính gần đúng giá trị của tích phân sau 1 3 4 I=∫ 2x + 1 dx −1 theo quy tắc hình thang . Chương trình sau đây cho giá trị gần đúng I = 0,693 đạt độ chính xác đã cho sau 4 lần lặp void main() { int i , j, n, lanlap; double h, a, b, epsilon, st, ss; clrscr(); printf("TICH PHAN BANG PP HINH THANG HAM y=1/(1+x)\n"); printf("Nhap can duoi a va can tren b: "); scanf("%lf %lf",&a,&b); printf("Nhap so epsilon: "); scanf("%lf",&epsilon); printf("Nhap so lan lap: "); scanf("%d",&lanlap); 59
ss = (f(a) + f(b)) * (b - a) / 2 ; j = 1; n = 1; do { st = ss; n *= 2; h = (b-a) / n; ss = (f(a) + f(b)) / 2 ; for (i = 1 ; i epsilon) && (j <= lanlap)); printf("I = %10.6lf tren [%lf,%lf] voi n=%d \n",ss,a,b,n); getch(); } 60
4.3 BÀI TẬP CHƯƠNG 4 1> Hãy sử dụng công thức Simpson tổng hợp tính các tích phân sau đây với h=0,1 1 2 1 1 x1/7 ex dx, dx, dx ∫ ∫2 ∫ 2 0 1 (3x - 2) 0 1+ x 2> Hãy sử dụng công thức hình thang tổng hợp để tính tích phân sau với bước h=0,1 1 1 ∫ dx 0 4+ sin( 20x) Bước h này đã đủ nhỏ để kết quả tính toán đạt được độ chính xác đến 3 chữ số thập phân? 3> Phải chia khoảng [0; 1] ít nhất thành bao nhiêu phần bằng nhau để việc tính tích phân trong bài tập 4.1 bằng công thức Simpson tổng hợp đạt được độ chính xác đến 4 chữ số thập phân? Tính tích phân đó với n tìm được. 4> Viết chương trình tính gần đúng tích phân xác định bằng công thức Simpson tổng hợp có tự động làm tăng độ chính xác theo yêu cầu. Kiểm tra chương trình bằng cách tính các tích phân sau với độ chính xác 10-6 1 2 1 1 x1/7 ex dx, dx, dx ∫ ∫2 ∫ 2 0 1 (3x - 2) 0 1+ x 61
5 CH ƯƠNG 5: GIẢI PHƯƠNG TRÌNH f(x) = 0 5.1 ĐẶT VẤN ĐỀ 5.1.1 Bài toán Cho hàm f(x) liên tục trên đoạn [a, b] hoặc trên khoảng vô hạn và phương trình f(x) = 0 (5.1) chỉ có nghiệm thực cô lập, tức là với mỗi nghiệm thực tồn tại một lân cận không chứa các nghiệm thực khác của phương trình 5.1.2 Các bước giải Việc tính gần đúng nghiệm thực của (5.1) được tiến hành theo hai bước: 1. Tách nghiệm hay tìm khoảng cách li nghiệm (a, b) chỉ chứa một nghiệm của (5.1) 2. Kiện toàn nghiệm: tính gần đúng nghiệm với độ chính xác cho trước 5.1.3 Tách nghiệm ¾ Cơ sở của phương pháp tách nghiệm Hệ quả của định lí Lagrange : Nếu hàm f(x) xác định và liên tục trên [a, b], f(a)f(b) < 0 và có f ‘(x) giữ dấu không đổi trên (a, b) thì tồn tại duy nhất một nghiệm thực x* (a, b) của phương trình (5.1). ¾ Các phương pháp tách nghiệm • Phương pháp giải tích: lập bảng xét dấu của đạo hàm cấp một f ‘(x) rồi tìm các khoảng (a, b) thỏa mãn hệ quả trên Ví dụ: Tìm các khoảng chứa các nghiệm cô lập của phương trình f(x) = x3 – 6x + 2= 0 Giải: f ‘(x) = 3(x2 – 2) có bảng xét dấu sau x -∞ -3 -√2 -1 0 1 √2 3 +∞ f ‘(x) + + - - - - + + f(x) - -7 7,66 7 2 -3 -3,66 11 + Vậy phương trình trên có ba nghiệm x1(-3, -√2), x2(0, 1) và x3(√2, 3) • Phương pháp hình học 12 Vẽ đồ thị hàm y = f(x), từ đó tìm ra 10 8 các khoảng cách li (a, b) chứa các 6 giao điểm của đồ thị với trục hoành 4 2 y Ví dụ: 0 Tìm các khoảng chứa các nghiệm cô -4 -2-2 0 2 4 -4 lập của phương trình -6 f(x) = x3 – 6x + 2= 0 -8 Vẽ đồ thị hàm y = x3 – 6x +2 62
Nhìn vào đồ thị ta thấy đường cong cắt trục hoành tai ba điểm phân biệt trên 3 khoảng rời nhau: (–4, –1), (0, 1) và (1, 3) 5.2 CÁC PHƯƠNG PHÁP KIỆN TOÀN NGHIỆM 5.2.1 Phương pháp chia đôi ¾ Thuật toán: Giả sử tách được khoảng (a, b) chứa nghiệm x*. Chia đôi (a, b) rồi xét dấu của f((a+b)/2) - Nếu f((a+b)/2) = 0 thì x* = (a+b)/2 là nghiệm đúng của phương trình f(x) = 0 (5.2) - Nếu f((a+b)/2) ≠ 0 thì x1 = (a+b)/2 là nghiệm gần đúng ở lần lặp thứ nhất. Chọn khoảng (a1, b1) chứa nghiệm có một trong hai đầu mút là (a+b)/2 còn đầu mút kia là a hay b Lặp lại quá trình trên đến khi tìm được nghiệm đúng x* hoặc nghiệm gân đúng xi đạt được độ chính xác mong muốn. ¾ Sự hội tụ và sai số: Sử dụng phương pháp chia đôi liên tiếp ta nhận được dãy khoảng lồng nhau {(ai, bi)} hữu hạn nếu x* là điểm giữa của khoảng thứ i hay vô hạn nhỏ dần: i ai < x* < bi, f(ai).f(bi) < 0, bi - ai = (b – a)/2 (5.3) Khi i → ∞ , do sự liên tục của f(x) nên lim bi = lim ai = x* a+ b b− a b− a ε =x * − x < i i −x ≤ i i = (5.4) i 2 i 2 2i+ 1 ¾ Sơ đồ khối Bắt đầu Nhập a, b, ε x=(a+b)/2 f(a)f(x) <0 b=x a=x b-a<ε Dừng 63
5.2.2 Phương pháp lặp đơn ¾ Thuật toán Giả sử tách được khoảng (a, b) chứa nghiệm x* Biến đổi tương đương phương trình f(x) = 0 Ù x = φ(x) (5.5) Chọn x0 (a, b) bất kì rồi tính các xâp xỉ liên tiếp nhờ phép lặp xi+1 = φ(xi) (5.6) Nếu dãy nghiệm gần đúng {xi} hội tụ về điểm ξ thì đó chính là nghiệm x* của (5.5) ¾ Điều kiện hội tụ và sai số Định lí: Nếu φ(x) khả vi liên tục trên [a, b], φ(x) ∈ [a, b] và |ω’(x)| ≤ q < 1 thì: Phương trình (5.5) có nghiệm duy nhất x* (a, b) Phép lặp (5.6) hội tụ về x* q q i x− x * ≤ x− x hay x− x * ≤ x− x (5.7) i 1− q i i− 1 1− q 1 0 ¾ Sơ đồ khối Bắt đầu Nhập x0, ε, q x = φ(x0), ss = |x-x0|, x0 = x ss = ss(1-q)/q ss ≤ ε Dừng Ví dụ: Giải phương trình f(x) = x3 – 6x + 2 = 0 trên (0; 1) bằng phương pháp chia đôi và lặp đơn với độ chính xác ε = 0,01 Chia đôi Lặp đơn 3 2 xi 1 x 1 xi+ 1=ϕ() x i = + ,| ϕ '( xi ) |= ≤ < 1, x0 = 0, 5 i xi 6 3 2 2 0 0,5 0,354167 1 0,25 0,340737 2 0,375 0,339927 3 0,3125 0,339820 4 0,34375 0,339877 5 0,328125 0,339877 6 0,335938 Nhận xét: Phương pháp chia đôi đơn giản nhưng hội tụ chậm hơn phương pháp lặp đơn 64
Dưới đây ta luôn luôn giả thiết rằng f(x) khả vi liên tục cấp hai trên [a, b] và f′(x), f″x) giữ dấu không đổi trên đó (không giảm tổng quát, ta coi f″x) > 0 nếu không ta xét phương trình g(x) = 0, với g(x) = –f(x) 5.2.3 Phương pháp dây cung ¾ Thuật toán: Giả sử f”(x) > 0 trên [a, b] Thay cung đường cong trên [a, b] bởi dây trương cung ấy và xem hoành độ giao điểm của dây cung với trục hoành là nghiệm gần đúng của phương trình (1). - Trường hợp 1: f′(x) > 0 Khi đó f(a) 0, đường cong y = f(x) là lõm trên [a, b]. Phương trình dây trương cung qua hai điểm A(a, f(a)) và B(b, f(b)) là y− f (a) x− a = f (b)− f (a) b− a Tại giao điểm với trục hoành (x, 0) ta có 0− f() a x− a f(a).(b− a) = ⇒ x = a − ∈ (a, b) f()() b− f a b− a f(b)− f(a) Thay a mới bởi x vừa tìm được. Quá trình lặp lại cho đến khi x đạt độ chính xác cho trước f (xi )(b− xi ) xi+ 1= x i − (5.8) f (b)− f (xi ) - Trường hợp 2: f′(x) 0, đường cong y = f(x) là lồi trên [a, b]. Phương trình dây trương cung qua hai điểm A(a, f(a)) và B(b, f(b)) là y− f (b) x− b = f (b)− f (a) b− a Tại giao điểm với trục hoành (x, 0) ta có 0− f() b x− b f(b).(b− a) = ⇒ x = b − ∈ (a, b) f()() b− f a b− a f(b)− f(a) Thay b mới bởi x vừa tìm được. Quá trình lặp lại cho đến khi x đạt độ chính xác cho trước f (xi )(x i − a) xi+ 1= x i − (5.9) f (xi )− f (a) Nhận xét: Công thức (5.8) và (5.9) còn đúng khi f′′(x) < 0 65
Nghiệm xấp xỉ ban đầu x0 là một trong hai điểm đầu mút a hoặc b được chọn theo qui tắc: f(x0).f′′(x0) 0 f(a) > 0 ¾ Sự hội tụ và sai số: - Dễ thấy dãy {xi} trong (5.8) hay (5.9) là đơn điệu và bị chặn bởi x* nên luôn có giới hạn.ξ Chuyển qua giới hạn hai công thức này ta thấy giá trị giới hạn ξ thoả mãn phương trình (5.1), mà trên [a, b] nó chỉ có duy nhất nghiệm, vậy ξ ≡ x*. - Sử dụng các giả thiết đã cho về hàm f(x) và công thức số gia hữu hạn ta nhận được đánh giá độ chính xác ở bước lặp thứ i: M− m =x − x * ≤ x− x (5.10) l i i+ 1 m i+ 1 i trong đó 0 0 trên (1,1; 1,5) nên ta chọn xấp xỉ ban đầu x0 =1,1 ( = a) rồi lặp theo (5.8) Với m = f′(1,1) =2,99 và M = f′(1,4) = 5,12 ta ước lượng được sai số sau mỗi bước lặp theo (5.10) 3 2 i xi f(xi) = xi – 0,2xi – 0,2xi - 1,2 |xi –xi-1| εi 0 1,1 –0,08254 1 1,18254 –0,06252 0,08254 0,058802 2 1,19709 –0,01056 0,01455 0,010362 3 1,19952 –0,00175 0,00243 0,001730 Vậy x2=1,19952 là nghiệm gần đúng phải tìm (nghiệm đúng là 1,2). 66
5.2.4 Phương pháp tiếp tuyến (Newton) ¾ Thuật toán: Với các giả thết về f(x) như ở phương pháp dây cung, ta chọn nghiệm xấp xỉ ban đầu x0 là một trong hai điểm đầu mút a hoặc b theo qui tắc f(x0).f′′(x0) > 0 rồi viết phương trình tiếp tuyến của đường cong y = f(x) tại điểm M0(x0, f(x0): y = f′(x0) (x – x0) + f(x0) Hoành độ giao điểm của tiếp tuyến với trục hoành là f() x0 x1= x 0 − f'( x0 ) f′′(x0) > 0 , f(b) > 0 Tiến hành lặp theo công thức f() xi xi+ 1= x i − (5.11) f'( xi ) đến khi nghiệm gần đúng đạt được độ chính xác đã cho. ¾ Sự hội tụ và sai số: - Sử dụng các giả thiết đã cho về hàm f(x) trên [a, b] và nhận xét rằng dãy {xi} trong (5.11) là dãy đơn điệu nên tồn tai giới hạn lim x i = ξ . i→∞ Rõ ràng ξ là nghiệm của (5.1), do tính duy nhất nghiệm trên (a, b) nên ξ ≡ x*. - Sử dụng các giả thiết đã cho về hàm f(x) và công thức khai triển Taylo ta nhận được đánh giá độ chính xác M 2 2 xi+ 1 − x * ≤ xi+ 1− x i (5.12) 2m1 trong đó 0 < m1≤ |f′(x)| và |f′′(x)| ≤ M2 với mọi x ε [a, b] Ví dụ: Tìm nghiệm gần đúng của phương trình x3 – 0,2x2 – 0,2x - 1,2 = 0 trên (1,1;1,4) bằng phương pháp tiếp tuyến với sai số ε = 0,01 Giải: do f(1,4) cùng dấu với f′′(x) nên ta chọn x0 = 1,4 rồi tiến hành lặp theo (5.11) với m1 = f(1,1) = 2,99 và M2 = |f′′(1,4)| = 8: i xi |xi-xi-1| εi f(xi) f′(xi) f(xi)/f′(xi) 0 1,4 0,872 5,12 0,17031 1 1,22969 0,17031 0,03880 0,11109 3,84454 0,02890 2 1,20079 0,02890 0,00112 67
8 2 Dừng ở bước lặp thứ hai ta có ε =|x − x * | ≤ 1,, 22969− 1 20078≤ 0,, 00112 0) printf("f(a) va f(b) cung dau\n"); else { if(f(a)*D2f(a) > 0) x0 = a; else x0 = b; sbl = 1 ; printf("CAC KET QUA TRUNG GIAN\n"); do { x = x0 - f(x0) / Df(x0); printf("x[%d] = %f\n" , sbl , x); ss = (M / (2 * m)) * (x - x0) * (x - x0); x0 = x ; sbl++ ; } while (ss > epsilon) ; printf("NGHIEM PHUONG TRINH x = %f\n " , x); } Nhận xét: - Độ phức tạp thời gian của thuật toán là O(n). - Phương pháp tiếp tuyến là trường hợp riêng của phương pháp lặp đơn. 68
¾ Sử dụng GOAL SEEK trong EXCEL tìm nghiệm của phương trình f(x) = 0 Chọn ô chứa phương trình ràng buộc cho hộp Set cell và ô chứa nghiệm, cho hộp By changing cell rồi nhập giá trị cho vế phải của phương trình 5.3 GIẢI HỆ PHƯƠNG TRÌNH PHI TUYẾN 5.3.1 Lập công thức: Cho hệ phi tuyến ⎧f1() x 1 , ,x n = 0 ⎪ ⎪f2() x 1 , ,x n = 0 ⎨ (5.13) ⎪ ⎪ ⎩fn() x 1 , ,x n = 0 trong đó các hàm fi(x1, , xn), i = 1, , n có các đạo hàm riêng cấp 2 liên tục và bị chặn. (0) (0) (0) Giả sử X = (x1 , , xn ) là nghiệm xấp xỉ ban đầu của hệ (5.13). Khai triển Taylor n hàm fi(x1, , xn), i = 1, , n tại điểm X(x1, x2, , xn) thuộc lân cận điểm X(0) và bỏ qua các số hạng bậc hai trở đi, ta có hệ: n (0) (0) ⎧ (0) (0) ∂f1 (x 1 , ,xn ) (0) ⎪ f1 (x , ,xn )≈ f 1 (x 1 , ,xn ) + (xj− x j ) = 0 ∑ ∂x ⎪ j= 1 j ⎨ ⎪ n ∂f (x(0) , ,x(0) ) f (x , ,x )≈ f (x(0) , ,x(0) ) + n 1 n (x− x(0) ) = 0 ⎪ n 1 n n 1 n ∑ ∂x j j ⎩⎪ j= 1 j 69
n (0) (0) ⎧ ∂f1 (x 1 , ,xn ) (0) (0) (0) ⎪ (xj− x j ) = - f1 (x 1 , ,xn ) ∑ ∂x ⎪ j= 1 j ⇔ ⎨ (5.14) ⎪ n ∂f (x(0) , ,x(0) ) n 1 n (x− x(0) ) = - f (x(0) , ,x(0) ) ⎪ ∑ ∂x j j n 1 n ⎩⎪ j= 1 j (0) (0) (0) (0) -1 ⎡ ∂f1 (x 1 , ,xn ) ∂f1 (x 1 , ,xn ) ⎤ (0) ⎢ ⎥ ⎡ (0) (0) ⎤ ⎡x1− x 1 ⎤ ∂x ∂x - f1 (x 1 , ,xn ) ⎢ ⎥ ⎢ 1 n ⎥ ⎢ ⎥ ⇔ = ⎢ ⎥ ⎢ ⎥ (0) (0) (0) (0) ⎢ ⎥ ⎢ (0) ⎥ ⎢∂fn (x 1 , ,xn ) ∂fn (x 1 , ,xn )⎥ ⎢ (0) (0) ⎥ ⎣xn− x n ⎦ ⎣- fn (x 1 , ,xn )⎦ ⎢ ∂x ∂x ⎥ ⎣ 1 n ⎦ (0) (0) (0) (0) -1 ⎡ ∂f1 (x 1 , ,x n ) ∂f1 (x 1 , ,xn ) ⎤ ⎡ (0) ⎤ ⎢ ⎥ ⎡ (0) (0) ⎤ ⎡x1 ⎤ x1 ∂x ∂x f1 (x 1 , ,xn ) ⎢ ⎥ ⎢ ⎥ ⎢ 1 n ⎥ ⎢ ⎥ ⇔ = − ⎢ ⎥ (5.15) ⎢ ⎥ ⎢ ⎥ (0) (0) (0) (0) ⎢ ⎥ ⎢ ⎥ ⎢ (0) ⎥ ⎢∂fn (x 1 , ,xn ) ∂fn (x 1 , ,xn )⎥ ⎢ (0) (0) ⎥ ⎣x n ⎦ ⎣x n ⎦ ⎣ fn (x 1 , ,xn )⎦ ⎢ ∂x ∂x ⎥ ⎣ 1 n ⎦ LÆp theo (5.15) ta cã ((k− 1) (k− 1) (k− 1) (k− 1) -1 ⎡∂f1 (x 1 , ,xn ) ∂f1 (x 1 , ,xn ) ⎤ ⎡ (k) ⎤ ⎡ (k− 1) ⎤ ⎢ ⎥ ⎡ (k− 1) (k− 1) ⎤ x1 x1 ∂x ∂x f1 (x 1 , ,xn ) ⎢ ⎥ ⎢ ⎥ ⎢ 1 n ⎥ ⎢ ⎥ = − ⎢ ⎥ (5.16) ⎢ ⎥ ⎢ ⎥ (k− 1) (k− 1) (k− 1) (k− 1) ⎢ ⎥ (k) (k− 1) (k− 1) (k− 1) ⎢ ⎥ ⎢ ⎥ ⎢∂fn (x 1 , ,xn ) ∂fn (x 1 , ,xn )⎥ ⎢ ⎥ ⎣x n ⎦ ⎣x n ⎦ ⎣ fn (x 1 , ,xn )⎦ ⎢ ∂x ∂x ⎥ ⎣ 1 n ⎦ - Người ta đã chứng minh rằng nếu nghiệm gần đúng X(0) chọn đủ gần nghiệm đúng X* của hệ (5.13) thì phép lặp (5.16) sẽ hội tụ về X* khi k→ ∞ 5.3.2 Các bước giải hệ phi tuyến bằng phương pháp lặp Newton-Raphson Bước 1: Tính ma trận các đạo hàm riêng ⎡ ∂f ∂f ⎤ 1 1 ⎢∂x ∂x ⎥ ⎢ 1 n ⎥ JX()= ⎢ ⎥ (5.17) ⎢ ⎥ ∂f ∂f ⎢ n n ⎥ ⎣⎢∂x1 ∂x n ⎦⎥ rồi thế X = X(0) để tính J(X(0)): (0) (0) (0) (0) (0) Bước 2: Tính vectơ F(X ) = [f1(x1 , , xn ), , fn(x1 , , xn )] Bước 3: Thế X = X(0).vào (5.15), giải hệ (5.14) J(X(0)) Y = –F(X (0) ) bằng một trong các phương pháp đã biết, tìm được véc tơ Y(1) = X(1) - X(0) (1) (1) (1) Bước 4: Suy ra nghiệm xấp xỉ thứ nhất X = (x1 , , xn của hệ (5.13) X(1) = X(0) + Y(1) Tiếp tục lặp từ bước 1 đến bước 4 theo công thức 70
J(X(k-1)) [[X(k) – X(k-1)] = –F(X(k-1)) (5.17’) cho đến khi n ()k (k− 1 ) ∑ xi − xi 0 cho trước (5.18) i=1 Ví dụ: Giải hệ phương trình ⎧ x2 + x2 = 1 1 2 ⎨ 3 ⎩x1 − x2 = 0 với xấp xỉ ban đầu X(0) = [0,9; 0,5] và ε = 0,1. Bước 1: Tính ma trận J(X(0)) theo công thức (5.17) ⎡ ∂f1 ∂f1 ⎤ ⎢ ⎥ ∂x1 ∂x2 ⎡2x1 2 x2 ⎤ JX()= ⎢ ⎥ = ⎢ ⎥ ⎢∂f ∂f ⎥ 2 2 2 ⎣⎢3x1 − 1 ⎦⎥ ⎢ ⎥ ⎣∂x1 ∂x2 ⎦ rồi thế X = X(0) để tính J(X(0)): ⎛⎡2.0,9 2.0,5 ⎤⎞ ⎛⎡1,8 1 ⎤⎞ ()JX()()0 = ⎜ ⎟ = ⎜ ⎟ ⎜⎢ 2 ⎥⎟ ⎜⎢ ⎥⎟ ⎝⎣3.0,9− 1 ⎦⎠ ⎝⎣2,43− 1 ⎦⎠ Bước 2: Tính vectơ ⎡f() 0,9; 0,5 ⎤ ⎡0,92 + 0,52 − 1⎤ ⎡ 0,06 ⎤ FX ()0 = 1 = = ()⎢ ⎥ ⎢ 3 ⎥ ⎢ ⎥ ⎣f2 () 0,9; 0,5 ⎦ ⎣0,9− 0,5 ⎦ ⎣0,229⎦ Bước 3: Giải Gauss ⎡ 1, 8 1 − 0, 06 ⎤ ⎛ 2, 43− 1 − 0, 229⎞ ⎡− 0, 0683⎤ (0) ⎜ ⎟ ()JXY()= ⎢ ⎥ ⇔ ⎜ ⎟ ⇒Y = ⎢ ⎥ ⎣2, 43− 1 − o, 229⎦ ⎝ 0 1, 711 0, 11 ⎠ ⎣ 0, 0630 ⎦ (1) (1) (1) Bước 4: Tìm xấp xỉ thứ nhất X = (x1 x2 ) của nghiệm theo công thức (5.16) ()1 ( 0 ) ⎡0, 9⎤ ⎡− 0, 0683⎤ ⎡0, 831679⎤ XXY= + = ⎢ ⎥ + ⎢ ⎥ = ⎢ ⎥ ⎣0, 5⎦ ⎣0, 0630 ⎦ ⎣0, 562979⎦ Tính sai số theo công thức (5.18) 2 ∑ yi = −0,0683 + 0,0630 = 0,1313 > 0,1 i= 1 lặp lại từ bước 1 đến bước 4 với X(1). Bước 1: Tính ma trận ()1 ⎡1,663 1,126⎤ ()JX()= ⎢ ⎥ ⎣2,075− 1 ⎦ Bước 2: Tính vectơ 71
⎡f() 0,831679; 0,562979 ⎤ ⎡0,8316792 + 0,5629792 −1⎤ ⎡ 0, 008634⎤ FX(1) = 1 = = ()⎢ ⎥ ⎢ 3 ⎥ ⎢ ⎥ ⎣f2 () 0,831679; 0,562979 ⎦ ⎣⎢0,831679 − 0,562979 ⎦⎥ ⎣0, 012284 ⎦ Bước 3: Giải Gauss ()1 ⎡1,663 1,126 − 0,008634⎤ ⎡−2,075− 1 0,012⎤ ⎡− 0,0056⎤ ()JXY()= ⎢ ⎥ ⇔ ⎢ ⎥ ⇒Y = ⎢ ⎥ B ⎣2,075− 1 − 0,012284⎦ ⎣ 0 1,928 0,001 ⎦ ⎣ 0,0006 ⎦ (2) (2) (2) ước 4: Tìm xấp xỉ thứ hai X = (x1 x2 ) của nghiệm theo công thức (5.16) ()2 ()1 ⎡ 0,831679⎤ ⎡0,008634⎤ ⎡0,826062⎤ XXY= + = ⎢ ⎥ + ⎢ ⎥ = ⎢ ⎥ ⎣0,562979 ⎦ ⎣0,012284⎦ ⎣0,563608⎦ Tính sai số theo công thức (5.18) 2 (2) ⎡0,826062⎤ ∑ yi = − 0056 + 0,0006 = 0,0062 < 0,1 dừng ở bước 2 với X = ⎢ ⎥ là i= 1 ⎣0,563608⎦ nghiệm gần đúng. Nhận xét: - Khi số phương trình lớn hơn 2, để việc giải lặp hệ (5.17’) J(X(k-1)) [[X(k) – X(k-1)] = –F(X (k-1) ) bằng máy tính được thuận tiện, người ta viết riêng các chương trình con tính ma trận J(X), véc tơ F(X), và giải hệ phương trình tuyến tính (5.17’) bằng phương pháp khử Gauss – Jordan rồi gọi các chương trình này vào mỗi bước lặp. - Sự hội tu của phép lặp về nghiệm đúng X* phụ thuộc vào việc chọn nghiệm X(0) ban đầu. 72
5.3.3 Sơ đồ khối và chương trình Bắt đầu Nhập X(0), ε Tính J(X(0)), F(X(0)) Giải hệ J(X(0))Y= –F(X(0)) X = X(0) + Y n )k( (k− 1) ss=∑ xi − xi i= 1 X(0) = X ss<ε Dừng #include <stdio.h) void main () { float epsilon , ss; int i, m, sbl, L ; // so buoc lap L toi da mang1 F, x, x0, y; matran J; clrscr(); printf("NHAP SO PHUONG TRINH CUA HE m = "); scanf("%d",&m); printf("NHAP VECTO NGHIEM BAN DAU XO\n"); for (i = 1; i <= m; i++) { printf("x0[%d] = ", i); scanf("%f",&x0[i]); } printf("SO BUOC LAP = "); scanf("%d",&L); printf("SAI SO epsilon = "); scanf("%f",&epsilon); sbl = 0; do { Nhap_vecto_F(m, x0, F); 73
Tao_ma_tran_J(m, x0, F, J); Gauss(m, J, y); printf("\n"); for(i = 1; i epsilon) && (sbl < L)); printf("\nNGHIEM PHUONG TRINH SAU %d LAN LAP VOI SAI SO %7.6f:\n",sbl, ss); for(i = 1; i <= m; i++) printf("x[%d] = %10.6f ",i,x[i]); } 5.4 PHƯƠNG PHÁP LẶP SEIDEL Bước 1: Đưa hệ (5.13) về hệ tương đương ⎧x1= g 1() x 1 , ,x n ⎪ ⎨ (5.19) ⎪ ⎩xn= g n() x 1 , ,x n Bước 2: Kiểm tra điều kiện hội tụ n ∂g ∑ i <1, i = 1,n (5.20) j= 1 ∂x j (0) (0) (0) Bước 3: Thay xấp xỉ ban đầu X = (x1 , ,xn ) của nghiệm vào vế phải của (5.19) (1) (1) (1) để tìm xấp xỉ thứ nhất X = (x1 , ,xn ) ⎧x()1 = g() x ()0 ,x ()0 , ,x (0) ⎪ 1 1 2 n ⎪x()1 = g() x ()1 ,x ()0 , ,x ()0 ⎨ 2 2 1 2 n (5.21) ⎪ ⎪ ()1 ()1 ()1 ()1 (0 ) ⎩xn = gn () x1 ,x2 , ,xn− 1 ,x n (1) (1) Như vậy các x1 , ,xi đã tìm được ở các hàng trên được thay vào các hàng dưới để (1) (1) tính xi+1 , ,xn , i = 1, , n-1. (0) (1) (1) (2) Bước 4: Lặp lại bước 2 bằng cách thay xi bởi xi ,thay xi bởi xi i = 1, , n tương (2) (2) (2) ứng để tìm xấp xỉ thứ hai X =(x1 , ,xn ) cho đến khi thoả mãn độ chính xác n ()l (l− 1 ) ∑ xi − xi < ε cho truoc < ε (5.22) i= 1 hoặc số bước lặp l = M cho trước. Ví dụ Giải hệ 74
3 2 2 ⎪⎧ x+ x − 10x − 2x + 1 = 0 1 1 1 2 ⎨ 2 2 ⎩⎪2x1 + 3x2 − 20x2 − 4 = 0 với ε = 0,03 ; M = 10 và X(0) = (0 0). Bước 1: Đưa hệ đã cho về hệ tương đương theo công thức (5.19) 3 2 2 ⎪⎧x= 0,1x + 0,1x − 0,2x + 0,1 1 1 1 2 ⎨ 2 2 ⎩⎪x1= 0,1x 1 + 0,15x2 − 0,2 Bước 2: Kiểm tra điều kiện (5.20) với 0 0,03 Vì vậy ta chuyển sang bước 4. Bước 4: Thay xấp xỉ thứ nhất X(1) của nghiệm để tìm xấp xỉ thứ hai (2) (2) (2) X = (x1 ,x2 ) theo công thức (5.21): (2) 3 2 2 x1 = 0,1.0,1 + 0,1.0,1 – 0,2.( –0,199) + 0,1 = 0,09318 (1) 2 2 x2 = 0,1.0,093 + 0,15.( –0,199) – 0,2 = –0,1854109 Tính sai số theo công thức (5.22) 2 ()2 ()1 ∑ xi − xi = 0,093 − 0,1 + − 0,185 −() − 0,199 i= 1 = 0,021 < 0,03 Vậy ta dừng ở xấp xỉ thứ hai của nghiệm. 75
5.5 Sử dụng Solver trong EXCEL giải hệ phương trình phi tuyến Nhập các giá trị đoán nhận ban đầu cho x1, x2, , xn 1) Nhập các vế trái của các phương trình f1, f2, , fn vào các ô riêng được liên kết với giá trị của các ẩn 2) Chọn công cụ Solver (Data Tab Æ Solver) Đặt “Target Cell” cho ô có công thức = 0 Chọn “Equal to Value of 0” Đối với hộp “By Changing Cells,” trỏ tới các ô chứa giá trị ban đầu của x1, x2, , xn Thêm các ràng buộc (các phương trình) Nhấn nút “Solve” 76
5.6 BÀI TẬP CHƯƠNG 5 1> Muốn dùng phương pháp chia đôi giải phương trình xex – 2 = 0 trên đoạn [0; 1] với độ chính xác ε = 0,001 thì cần phải lặp bao nhiêu lần? 2> Đa thức f(x) = x3 – 2x – 5 có một nghiệm trong khoảng (2; 3). (a) Áp dụng bốn bước của phương pháp chia đôi để thu nhỏ khoảng chứa nghiệm với độ rộng còn 1/16. (b) Tính x3 và x4 bằng phương pháp dây cung, bắt đầu với xl= 3 và x2 = 2. Đánh giá sai số của các kết quả (c) Tính x2, x3, và x4 bằng phương pháp Newton với x1 = 2. Đánh giá sai số của các kết quả. (d) Tính x2, x3, và x4 bằng phương pháp lặp đơn với x1 = 2. Đánh giá sai số của các kết quả. 3> Giải hệ phương trình sau bằng phương pháp Newton f(x,y) = x2 + xy3 – 9 = 0 g(x,y) = 3x2y – y3 – 4 = 0. 77
6 CH ƯƠNG 6: CÁC PHƯƠNG PHÁP SỐ GIẢI PHƯƠNG TRINH VI PHÂN 6.1 ĐẶT VẤN ĐỀ Nhiều bài toán của khoa học kĩ thuật dẫn đến việc tìm nghiệm của phương trình vi phân (PTVP) thoả mãn một số điều kiện nào đó. Những phương pháp giải đúng chỉ áp dụng được cho một lớp rất hẹp các PTVP. Phần lớn các PTVP mô tả các hệ thống cơ học, vật lí, thuỷ lực, sinh học, đều không có lời giải đúng. Chương này sẽ nghiên cứu các phương pháp số giải gần đúng hai bài toán sau đây của PTVP. 6.1.1 Bài toán Cauchy (bài toán giá trị đầu) - Cho PTVP cấp một sau đây: y' = f(x,y), (x0 ≤ x ≤ x*, hàm f đã biết) (6.1a) Tìm hàm y(x) thoã mãn các điều kiện y(x0) = y0 (vế phải đã biết) (6.1b) - Tương tự như trên, bài toán Cauchy đối với PTVP cấp n là: Tìm hàm y(x) thoã mãn các điều kiện (n) (n-1) y = f(x,y, y′, , y ), (x0 ≤ x ≤ x*, hàm f đã biết) (6.2a) (n-1) (n-1) y(x0) = y0, y′(x0) = y′0 , , y (x0) = y0 (vế phải đã biết) (6.2b) - Đối với hệ PTVP cấp một, bài toán Cauchy là: Tìm các hàm y1(x), y2(x), , yn(x) thỏa mãn các điều kiện: ⎧ dy 1 = f (x, y , y , , y ) ⎪ dx 1 1 2 n ⎪ dy ⎪ 2 = f (x, y , y , , y ) ⎨ dx 2 1 2 n (6.3a) ⎪ ⎪ dy n ⎪ = fn (x, y 1 , y 2 , , y n ) ⎩ dx y1(x0) = a1, y2(x0) = a2, , yn(x0) = an, (6.3b) (n-1) Nhận xét: Bằng cách đặt y1(x) = y′, y2 = y′′, , yn-1 = y ta đưa được PTVP cấp n (6.2a) về hệ PTVP cấp một tương đương: Tìm các hàm y(x), y1(x), y2(x), , yn-1(x) ⎧dy = y(,) x y ⎪dx 1 ⎪ dy1 ⎪ = y2 (,) x y ⎪ dx thoả mãn hệ PTVP cấp một: ⎨ (x0 ≤ x ≤ x*) (6.4a) dy ⎪ n− 2 = y(,) x y ⎪ dx n− 1 ⎪ dy ⎪ n− 1 = f( x , y , y , , y ) ⎩⎪ dx 1 n− 1 thoả mãn các điều kiện đầu: y(x0) = a1, y1(x0) = a2, , yn-1(x0) = an (6.4b) Vì vậy ta chỉ cần giải các bài toán Caushy (6.1a), (6.1b) và (6.3a), (6.3b) là đủ 78
6.1.2 Bài toán biên hai điểm tuyến tính đối với PTVP cấp hai: Tìm hàim y(x) thoả mãn các điều kiện y′′ + p(x)y′ + q(x)y = f(x), (a ≤ x ≤ b) (6.5a) α0y(a) + α1y′(a) = A (6.5b) β0y(b) + β1y′(b) = B (6.5c) 6.1.3 Các phương pháp số giải bài toán Cauchy Trong các phương pháp số giải bài toán Cauchy (6.1a), (6.1b), người ta tìm nghiệm tại các điểm x0 0, thay việc sử dụng khai triển Taylor hàm y(x0+h) bằng việc tính giá trị của f(x, y) tại một dãy điểm (xi, yi): 6.2 CÁC PHƯƠNG PHÁP SỐ GIẢI BÀI TOÁN CAUCHY Tìm nghiệm bằng số của hàm y(x) thoã mãn các điều kiện y' = f(x,y), x0 ≤ x ≤ x*, hàm f đã biết) (6.7a) y(x0) = y0 (vế phải đã biết) (6.7b) 6.2.1 Phương pháp Euler ¾ Lập công thức Giả sử cho các mốc xk = x0 + kh, k = 0,1, 2, , n; hằng số h > 0 được chọn sao cho xn = x*. Kí hiệu giá trị gần đúng của nghiệm y(x) tại các mốc xk là yk ≈ y(xk) Khai triển Taylor hàm y(x) đến cấp một tai lân cận điểm xk ()x− x 2 y()()()() x= y x + x − x y' x + k y'' ( c ), c∈ ( x , x ) k k k 2! k Thay x = xk+1 = xk + h và y′(xk) = f(xk, y(xk)) vào đẳng thức trên ta có h 2 y( x )= y ( x ) + hf ( x , y ( x ))+ y'' ( c ) (6.8) k+ 1 k k k 2 Bỏ qua số hạng cuối cùng ở vế phải và thay y(xk)≈ yk vào đẳng thức trên ta được yk+1 ≈ yk + hf(xk,yk) (6.9) Như vậy xuất phát từ y0 ta có thể lần lượt tìm gia trị gần đúng y1, y2,, , yn ¾ Sai số Từ (6.8) người ta chứng minh được rằng sai số địa phương M ε =y(x ) − y ≤h2 = O() h 2 (6.10) k k k 2 79
trong đó |y′′(c)| ≤ M., c (xk, xk+1) và sai số toàn phần ε = O(h) ¾ Ví dụ Dùng phương pháp Euler tìm giá trị gần đúng y(1) từ phương trình vi phân sau đây y' = xy / 2 thoả mãn điều kiện đầu y(0) = 1 với sai số địa phương cỡ 10-2 Giải: Để đạt được sai số địa phương cỡ 10-2 ta chọn bước h = 0,1 hay chia [0, 1] thành 10 khoảng bằng nhau i xi yi f(xi, yi) hf(xi, yi) Nghiệm đúng 0 0 1 0 0 1 1 0,1 1 0,05 0,00 1,03 2 0,2 1,01 0,10 0,01 1,01 3 0,3 1,02 0,15 0,02 1,02 4 0,4 1,03 0,21 0.02 1,04 5 0,5 1,05 0,26 0,03 1,06 6 0,6 1,07 0,32 0,03 1,09 7 0,7 1,11 0,39 0,04 1,13 8 0,8 1,15 0,46 0,05 1,17 9 0,9 1,19 0,58 0,05 1,22 10 1 1,25 1,28 ¾ Sơ đồ khối và chương trình Bắt đầu Nhập x0, y0, x, n h = (x – x0) / n k = 1, 2, , n yk = yk-1 + hf(xk-1, yk-1) xk = xk-1 + h Dừng 80
printf("\nSo diem chia n = "); scanf("%d", &n); printf("Gia tri diem dau x0 = "); scanf("%lf", &x[0]); printf("Gia tri diem cuoi x* = "); scanf("%lf", &x[n]); printf("Gia tri ban dau y0 = "); scanf("%lf", &y[0]); h = (x[n]-x[0]) / n; for (i=1; i<=n; i++) { y[i] = y[i-1] + h*f(x[i-1], y[i-1]); x[i] = x[i-1] + h; } Nhận xét: - Phương pháp Euler đơn giản, dễ lập trình nhưng có độ chính xác thấp. - Có thể áp dụng phương pháp Euler để giải bài toán Cauchy cho hệ phương trình vi phân cấp một, chẳng hạn với hệ hai phương trình ⎧y' = f (x, y,z) ⎨ ' ⎩z= g(x, y,z) và điều kiện đầu y(x0) = y0; z(x0) = z0 công thức (6.3a) có dạng sau ⎧yk+ 1= y k + hf(xk , y k ,z k ) ⎨ k= 1, 2, , n (6.11) ⎩zk+ 1= z k + hg(xk , y k ,z k ) ¾ Ví dụ: Giải hệ phương trình sau trên [0; 0,6] ⎧y' = (z − y)x ⎨ ' ⎩z= (z + y)x với điều kiện đầu y(0) = z(0) = 1,0000.và h = 0,1 Giải: Áp dụng công thức (6.11) ta tính được bảng sau: i xi yi zi f(xi, yi, zi) hf(xi, yi, zi)g(xi, yi, zi) hg(xi, yi, zi) 0 0 1 1 0 0 0,0 0 1 0,1 1 1 0 0 0,2000 0,0200 2 0,2 1 1,0200 0,0040 0,0004 0,4040 0,0404 3 0,3 1,0004 1,0604 0,0180 0,0018 0,6182 0,0618 4 0,4 1,0022 1,1222 0,0480 0,0048 0,8498 0,0850 5 0,5 1,0070 1,2072 0,1001 0,0100 1,1071 0,1107 6 0,6 1,0170 1,3179 6.2.2 Phương pháp Euler cải tiến ¾ Lập công thức Áp dụng công thức Newton – Lebnitz ta có 81
xk+ 1 ' y()()() xk+ 1 − y xk = ∫ y x dx xk rồi tính gần đúng tích phân xác định ở vế phải theo công thức hình thang ta nhận được h h 3 y( x )− y ( x ) = [y '( x )+ y '( x )]− y '''( c ) k+ 1 k 2 k k+ 1 12 Bỏ qua số hạng cuối ở vế phải và thay y(xk) ≈ yk, y(xk+1) ≈ yk+1, y′(xk) ≈ f(xk, yk), y′(xk+1) ≈ f(xk+1, yk+1) ta nhận được h y− y = [f ( x , y )+ f ( x , y )] ⇔ k+ 1 k 2 k k k+ 1 k + 1 h y= y + [f ( x , y )+ f ( x , y )], i = 0,n (6.12) k+ 1 k 2 k k k+ 1 k + 1 ¾ Sai số Người ta chứng minh được rằng phương pháp Euler cải tiến có sai số toàn phần O(h3) Nhận xét: Nhược điểm cơ bản của phương pháp Euler cải tiến là muốn tính yk+1 ta phải giải phương trình (6.12) đối với yk+1. Để khắc phục điều này, người ta sử dụng phối hợp hai phương pháp trên như sau: h y= y + [f(x ,y )+ f(x + h,y + hf(x ,y ))] (6.13) k+ 1 k 2 k k k k k k ¾ Ví dụ: Dùng phương pháp Euler cải tiến giải lại bài toán Cauchy trong ví dụ trước ta theo công thức (6.12) ta nhận được bảng sau i xi yi f(xi, yi) hf(xi, yi) zi=yi+hf(xi, yi) hf(xi+1, zi) Nghiệm đúng 0 0 1 0 0 1 1 0,1 1,0025 0,0511 0,0051 1 0,0050 1,0025 2 0,2 1,0100 0,1010 0,0101 1,0075 0,0101 1,0101 3 0,3 1,0227 0,1534 0,0153 1,0201 0,0153 1,0228 4 0,4 1,0408 0,2082 0.0208 1,0381 0,0208 1,0408 5 0,5 1,0646 0,2661 0,0266 1,0616 0,0265 1,0645 6 0,6 1,0943 0,3283 0,0328 1,0911 0,0327 1,0942 7 0,7 1,1305 0,3957 0,0396 1,1270 0,0394 1,1337 8 0,8 1,1738 0,4695 0,0470 1,1698 0,0468 1,1735 9 0,9 1,2248 0,5512 0,0551 1,2204 0,0549 1,2245 10 1 1,2812 1,2840 82