Đồ án Phương pháp bình phương tối thiểu lập công thức từ thực nghiệm

MỤC LỤC Mục lục…………………………………………………………………...1 Đặt vấn đề…………………………………………………………3 Đặt vấn đề…………………………………………………..3 Bài toán đặt ra………………………………………………4 Sai số trung bình phương và phương pháp bình phương tối thiểu tìm xấp xỉ tốt nhất với một hàm…………………………………...5 Sai số trung bình phương…………………………………...5 Định nghĩa………………………………………………….5 Ý nghĩa của sai số trung bình phương……………………...5 Xấp xỉ hàm theo nghĩa trung bình phương…………………7 Xấp xỉ hàm trong thực nghiệm bằng đa thức suy rộng……………9 Định nghĩa………………………………………………….9 Nội dung……………………………………………………9 Sai số của phương pháp…………………………………...11 Mở rộng trên hệ trực giao để đơn giản hóa kết quả……….13 Định nghĩa………………………………………….13 Tiếp cận lời giải…………………………………….13 Sai số của phương pháp…………………………….14 Chú ý……………………………………………….14 Xấp xỉ hàm trong thực nghiệm bằng đa thức đại số……………..16 Đặt vấn đề…………………………………………………16 Tiếp cận lời giải……………………………………….......16 Sai số trung bình…………………………………………..17 Áp dụng……………………………………………….......17 Chú ý………………………………………………………19 Xấp xỉ hàm trong thực nghiệm bằng đa thức trực giao………….24 Định nghĩa hệ hàm trực giao………………………..…….24 Đặt vấn đề…………………………………………………24 Nội dung của phương pháp………………………….…….25 Sai số của phương pháp……………………………..…….34 Thiết kế chương trình…………………………………………………..36 Code chương trình……………………………..……………………….38 Kết luận…………………………………………………………….......50 Tài liệu tham khảo………………………………………………….......51 PHƯƠNG PHÁP BÌNH PHƯƠNG TỐI THIỂU LẬP CÔNG THỨC TỪ THỰC NGHIỆM I. ĐẶT VẤN ĐỀ 1.1 Đặt vấn đề Có rất nhiều phương pháp khác nhau để lập những đa thức từ thực nghiệm mà ta đã biết đến như phép nội suy để lập đa thức cấp n:(x) (đại số hoặc lượng giác) xấp xỉ hàm số y=f(x) mà ta đã biết các giá trị của hàm này là y=y tại các điểm x=x. Phương pháp nội suy nói trên khi sử dụng trong thực tiễn thì có những điều cần cân nhắc là: Trong các đa thức nội suy (x) ta đòi hỏi) = y.Tuy nhiên sự đòi hỏi này không có ý nghĩa nhiều trong thực tế. Thật vậy các số ylà giá trị của hàm y=f(x) tai các điểm x=x mà trong thực tế chúng ta cho dưới dạng bảng và thường thu được từ những kết quả đo đạc hoặc tính toán trong thực hành. Bởi vậy những số y này nói chung chỉ xấp xỉ với các giá trị đúng f(x) của hàm y=f(x) tại x=x.Sai số mắc phải =y - f (x) nói chung khác không. Nếu buộc )=y thì thực chất đã đem vào bài toán các sai số của các số liệu ban đầu nói trên(chứ không phải là làm cho giá trị của hàm nội suy và hàm f(x) trùng nhau tại các điểm x=x) Để cho đa thức nội suy biểu diễn xấp xỉ hàm f(x)một cách sát thực đương nhiên cần tăng số mốc nội suy x(nghĩa là làm giảm sai số của công thức nội suy). Nhưng điều này lại kéo theo cấp của đa thức nội suy tăng lên do đó những đa thức nội suy thu được khá cồng kềnh gây khó khăn cho việc thiết lập cũng như dựa vào đó để tính giá trị gần đúng hoặc khảo sát hàm f(x). 1.2 Bài toán đặt ra. Chính vì những lý do đó mà phương pháp tìm hàm xấp xỉ có thể sẽ sát thực hơn thông qua hai bài toán: Bài toán 1(tìm hàm xấp xỉ) Giả sử đã biết giá trị y (i=1,2,…n) của hàm y=f(x) tại các điểm tương ứng x=x. Tìm hàm (x) Xấp xỉ với hàm f(x) trong đó (x) = với là những hàm đã biết, a là những hệ số hằng số. Trong khi giải quyết bài toán này cần chọn hàm sao cho quá trình tính toán đơn giản đồng thời nhưng sai số có tính chất ngẫu nhiên (xuất hiện khi thu được các số liệu y) cần phải được chỉnh lý trong quá trình tính toán. Trong bài toán tìm hàm xấp xỉ trên việc chọn dạng của hàm xấp xỉ là tùy thuộc ý nghĩa thực tiễn của hàm f(x). Bài toán 2 (tìm các tham số của một hàm có dạng đã biết). Giả sử đã biết dạng tổng quát của hàm Y = f(x, a, a, …, a). (1 – 2) Trong đó : a(i= 0,1,…,m) là những hằng số. Giả sử qua thực nghiệm ta thu được n giá trị của hàm y=y(i=1, 2, …, m) ứng với các giá trị x=x của đối. Vấn đề là từ những số liệu thực nghiệm thu được cần xác định các giá trị của tham số a, a, …, a để tìm được dạng cụ thể của biểu thức (1 – 2): y=f(x) về sự phụ thuộc hàm số giữa y và x. II. SAI SỐ TRUNG BÌNH PHƯƠNG VÀ PHƯƠNG PHÁP BÌNH PHƯƠNG TỐI THIỂU TÌM XẤP XỈ TỐT NHẤT VỚI MỘT HÀM. 2.1 Sai số trung bình phương: Những hàm trong thực nghiệm thu được thường mắc phải những sai số có tính chất ngẫu nhiên. Những sai số này xuất hiện do sự tác động của những yếu tố ngẫu nhiên vào kết quả thực nghiệm để thu được các giá trị của hàm. Chính vì lý do trên, để đánh giá sự sai khác giữa hai hàm trong thực nghiệm ta cần đưa ra khái niệm về sai số (hoặc độ lệch) sao cho một mặt nó chấp nhận được trong thực tế, một mặt lại san bằng những sai số ngẫu nhiên (nghĩa là gạt bỏ được những yếu tố ngẫu nhiên tác động vào kết quả của thực nghiệm). Cụ thể nếu hai hàm thực chất khá gần nhau thì sai số chúng ta đưa ra phải khá bé. Khái niệm về sai số nói trên không chú ý tới những kết quả có tính chất cá biệt nên được gọi là sai số trung bình phương. 2.2 Định nghĩa. Theo định nghĩa ta sẽ gọi là sai số (hoặc độ lệch) trung bình phương của hai hàm f(x) và trên tập X = {x,x,…,x}, nếu = . (2 – 1) 2.3 Ý nghĩa của sai số trung bình phương. Để tìm hiểu ý nghĩa của sai số trung bình phương ta giả thiết f(x), (x) là những hàm liên tục trên đoạn [a, b] và X = { x, x, …,x} là tập hợp các điểm cách đều trên [a, b]: a = x < x < … < x = b. Theo định nghĩa fích phân xác định ta có lim = , (2 – 2) Trong đó: = . (2 – 3) Giả sử │f(x) - │ có trên [a, b] một số hữu hạn cực trị và là một số dương nào đó cho trước. Khi đó trên [a, b] sẽ có k đoạn riêng biệt [a,b] ( i = 1, 2, …,k ) sao cho: │f(x) -│ ( với x [a, b], i = 1, 2, …, k) Gọi là tổng các độ dài của k đoạn nói trên. Với n đủ lớn và đủ bé, từ (2 – 2) ta suy ra < ( bé tùy ý). Từ (2 – 3) suy ra > Do đó : < ( b – a ) nghĩa là tổng đọ dài của các đoạn [a, b] sẽ bé tùy ý. Tóm lại: với đủ bé (n khá lớn) thì trên đoạn [a, b] ( trừ tại những điểm của những đoạn [a, b] mà có tổng độ dài bé tùy ý), ta có │f(x) - │ < , trong đó là một số dương tùy ý cho trước. Từ nhận xét trên ta rút ra nhận ý nghĩa thực tiễn của sai số trung bình phương như sau: Nếu sai số trung bình phương của hai hàm f(x) và trên tập hợp n điểm X [a, b] ( n đủ lớn) mà khá bé thì với tuyệt đại đa số giá trị của x trên [a, b] cho sai số tuyệt đối giữa f(x) và khá bé. 2.4 Xấp xỉ hàm theo nghĩa trung bình phương. Từ ý nghĩa của sai số trung bình phương nói trên Ta nhận thấy nếu các giá trị y ( i = 1, 2, …, n ) của hàm f(x) tại các điểm x và nếu sai số trung bình phương = khá bé thì hàm sẽ xấp xỉ khá tốt với hàm f(x). Cách xấp xỉ một hàm số lấy sai số trung bình phương làm tiêu chuẩn đánh giá như trên gọi là xấp xỉ hàm theo nghĩa trung bình phương. Rõ ràng: nếu hàm f(x) thu được bằng thực nghiệm (nghĩa là y f(x)) thì cách xấp xỉ nói trên đã san bằng những sai lạc tại từng điểm (nảy sinh do những sai số ngẫu nhiên của thực nghiệm). Đó là lý do giải thích lý do vì sao phương pháp xấp xỉ theo nghĩa trung bình phương được sử dụng rộng rãi trong thực tiễn. Ta xét trường hợp là phụ thuộc các tham số a, a, . . . ,a: = ( x; a, a, …,a ). (2 – 4) Trong số những hàm có dạng (2 – 4) ta sẽ gọi hàm (x) = ( x ; , , …, ) (2 – 5) là xấp xỉ tốt nhất theo nghĩa trung bình phương với hàm f(x) nếu sai số trung bình phương (x) với f(x) là bé nhất. Cụ thể là: (, , … , ) = min trong đó (a, a, …, a) = (2 – 6) Từ (2 – 6) ta nhận thấy (2 – 5) tương đương với đẳng thức: = min . (2 – 7) Từ đó việc tìm hàm xấp xỉ tốt nhất (trong số những hàm dạng (2 – 4) với hàm f(x)) sẽ đưa về tìm cực tiểu của tổng bình phương trong đó = y - Bởi vậy phương pháp tìm xấp xỉ tốt nhất theo nghĩa trung bình còn gọi là phương pháp bình phương tối thiểu để xấp xỉ hàm trong thực nghiệm. III. XẤP XỈ HÀM TRONG THỰC NGHIỆM BẰNG ĐA THỨC SUY RỘNG 3.1. Định nghĩa: Giả sử cho hệ hàm: , ,…,,…Ta sẽ gọi hàm là đa thức suy rộng cấp m nếu có dạng = (3 – 1) trong đó a, a, … ,a là các hệ số hằng số. Hệ hàm {} đã cho gọi là hệ cơ bản. 3.2. Nội dung. Theo phần trên về tìm hàm xấp xỉ giả sử đã biết n giá trị thực nghiệm y (i=1,2,…,n) của hàm y=f(x) tại các điểm tương ứng x. Khi đó việc tìm một đa thức suy rộng có dạng (3 – 1) mà xấp xỉ với hàm f(x) nói trên [a, b]{x,x,…,x} sẽ chuyển về việc tìm m+1 hệ số a trong (3 – 1). Để quá trình tính toán được đơn giản ta xét đa thức suy rộng với cấp m không lớn lắm. Tuy nhiên ta vẫn phải chọn n đủ lớn do đó có thể giả thiết n m+1. Khác với bài toán nội suy ở đây ta không cần xác định m+1 giá trị a từ n phương trình: y = ( i=1,2,…,n) (vì số phương trình thường nhiều hơn số ẩn). Ta sẽ áp dụng phương pháp bình phương tối thiểu để tìm đa thức suy rộng = xấp xỉ tốt nhất với hàm f(x) trên [a, b]. Trong (2 – 7) ta coi (x; a,a,…,a) = = Từ đó ta suy ra : (,,…,) là điểm cực tiểu của hàm m+1 biến F(a,a,….,a) = (3 – 2) Do đó (a ,a ,…,a) là nghiệm của hệ phương trình = 0 ; = 0 ; ……; = 0. Hoặc dạng tương đương với nó (3 – 3) . . . . . . . . . . . . . . . . . . Gọi là véc tơ n chiều với thành phần thứ i là Gọi y là véc tơ n chiều với thành phần thứ i là y Theo định nghĩa tích vô hướng các véc tơ ta có [y, ] = ; [ = (3 – 4) Do đó (3 – 3) được chuyển về dạng (3 – 5) . . . . . . . . . . . . . . . . Ta nhận thấy (3 – 5) là hệ (m + 1) phương trình đại số tuyến tính dùng để xác định m + 1 hệ số : ,,… trong đa thức xấp xỉ . Ma trận của hệ phương trình tuyến tính (3 – 5) có các phần tử là , do đó là một ma trận đối xứng (dựa vào tính chất giao hoán của tích vô hướng ).Ta sẽ gọi hệ phương trình (3 – 5) là hệ phương trình chuẩn. Định thức của hệ phương trình chuẩn có dạng G( = (3 – 6) Ta gọi định thức G=( là định thức Gram của hệ véc tơ . Mà ta đã biết: Nếu hàm cơ sở là hệ hàm độc lập tuyến tính trên [a, b] X = {x thì trong số những đa thức suy rộng cấp m có dạng (3 – 1) luôn tồn tại một đa thức suy rộng (3 – 1’) là xấp xỉ tốt nhất theo nghĩa trung bình phương đối với hàm f(x). Ngoài ra còn có thể chứng minh khi hệ cơ sở là những độc lập tuyến tính trên [a, b] thì G(> 0. Nghĩa là trong trường hợp này hệ phương trình chuẩn (3 – 5) có và duy nhất nghiệm ứng với các hệ số của đa thức (3 – 1’) xấp xỉ tốt nhất với hàm f(x) (theo nghĩa trung bình phương). Do vậy ta có thể cho rằng hệ hàm cơ sở nghĩa là hệ hàm đôc lập tuyến tính trên đoạn [a, b]. 3.3 Sai số của phương pháp. Cùng với việc tìm hàm xấp xỉ cho hàm f(x) ta cần đánh giá sai số hoặc độ lệch của nó đối với hàm f(x). Sai số ở đây hiểu theo nghĩa trung bình phương. Cụ thể là ta đi tìm đại lượng (3 – 7) Từ (3 – 1’) ta có = [ (3 – 8) Mặt khác = (3 – 9) Kết hợp (3 – 9) với (3 – 5) ta có: Thay kết quả trên vào (3 – 8) ta có: (3 – 10) Thay (3 – 10) vào (3 – 7) ta có: (3 – 11) trong đó là nghiệm của hệ phương trình chuẩn (3 – 5). 3.4 Mở rộng trên hệ trực giao (để đơn giản hóa kết quả). 3.4.1 Định nghĩa: Tuy nhiên để đơn giản hóa kết quả trên thì ta định nghĩa về hệ hàm trực giao như sau: Hệ hàm , , …, gọi là hệ trực giao trên tập X={, , …, x} nếu: (3 – 12) Số mà gọi là chuẩn của hàm trên tập hợp X. Trong trường hợp hệ hàm , ,…, trực giao mà (r = 0, 1, 2, …, m) thì hệ hàm được gọi là hệ trực chuẩn trên tập hợp X. 3.4.2 Tiếp cận lời giải. Từ một hệ cơ sở bất kỳ , , …, bao giờ cũng lập được một hệ trực chuẩn tương ứng , , … , sao cho mỗi hàm của hệ trực chuẩn là một tổ hợp tuyến tính của các hàm trong hệ cơ sở đã cho: (r = 0, 1, … , m) (3 – 13) Từ (3 – 5) và (3 – 12) ta nhận thấy rằng: Nếu , , … , là hệ trực giao thì đa thức xấp xỉ tốt nhất (3 – 1’) của f(x) có các hệ số cho bởi công thức: (i = 0, 1, … , m) hay (i = 0, 1, … , m) (3 – 14) Từ đó ta có: 3.4.3 Sai số của phương pháp. Dựa trên (3 – 11) ta suy ra sai số trung bình phương của đa thức xấp xỉ là: (3 – 15) Vì nên tổng: là một đại lượng đơn điệu tăng theo m. Do đó từ (3 – 15) ta suy ra sai số trung bình phương sẽ giảm khi m tăng. Tóm lại nếu cấp m của đa thức xấp xỉ (3 – 1’) (với hệ cơ sở , , … , là trực giao) càng lớn thì đa thức xấp xỉ f(x) càng tốt. 3.4.4. Chú ý. Một đặc điểm chú ý ở đây là: Trong trường hợp chung khi cần thay đổi cấp m của đa thức xấp xỉ (3 – 1’) thì hệ phương trình chuẩn (3 – 5) dùng để xác định các hệ số ,,… của đa thức hoàn toàn thay đổi. Do đó quá trình tình toán (giải hệ phương trình chuẩn) cần làm lại từ đầu. Tuy nhiên khi hệ hàm cơ sở là trực giao thì muốn thay đổi cấp m của đa thức xấp xỉ (3 – 1’) (chẳng hạn tăng từ m lên m+1) ta chỉ cần thêm số từ công thức (3 – 14). Còn các hệ số ,,… đã thu được cho đa thức vẫn dùng được cho đa thức: Nhận xét trên rất bổ ích về mặt thực hành tính toán vì khi muốn xấp xỉ một hàm thực nghiệm bằng một đa thức suy rộng cấp m (3 – 1’): do khuôn khổ của sự tính toán ta không có khả năng chọn ngay từ đầu số m đủ lớn. Khi đó nếu hệ hàm cơ sở , ,…, , … là một hệ trực giao thì khi xuất phát ta có thể chọn số m nhỏ (chẳng hạn m = 1 hoặc 2). Sau khi thực hành tính toán nếu thấy sai số trung bình phương tương ứng chưa đủ bé (so với yêu cầu) thì ta có thể tăng dần số m lên và tính thêm các hệ số bổ xung (từ công thức (3 – 14)). IV. XẤP XỈ HÀM TRONG THỰC NGHIỆM BẰNG ĐA THỨC ĐẠI SỐ. 4.1 Đặt vấn đề. Giả sử biết n giá trị thực nghiệm y(i = 1, 2, … , n) của hàm f(x) tại các điểm x tương ứng. Ta đặt vấn đề xấp xỉ hàm f(x) bởi một đa thức cấp m có dạng : P=a + ax + … + ax (4 – 1) 4.2 Tiếp cận lời giải. Để giải bài toán này ta áp dụng những kết quả tổng quát ở phần III, trong đó hệ hàm cơ sở có dạng: , , … , (4 – 2) khi đó từ (3 – 4) ta có: và (4 – 3) Dựa vào (3 – 5) ta suy ra các hệ số a của đa thức xấp xỉ (4 – 1) là nghiệm của hệ phương trình chuẩn có dạng sau: (4 – 4) 4.3 Sai số trung bình. Từ (3 – 7) và (3 – 11) ta suy ra sai số trung bình của đa thức xấp xỉ có dạng (4 – 4) là: (4 – 5) Về mặt thực hành, để tìm các hệ số của phương trình chuẩn (4 – 4) ta làm theo lược đồ trong bảng 1. Các hệ số vế trái của phương trình đầu tiên cho bởi các tổng ô lần lượt từ cột (1) đến cột (m), của phương trình thứ 2 cho bởi các tổng lần lượt từ cột 2 đến cột (m+1), …., còn các vế phải của (4 – 4) cho bởi các tổng ở lần lượt từ cột (2m+2) đến cột cuối cùng (3m+2). Bảng 1 x x x … x y xy xy … xy (1) (2) (3) (2m+1) (2m+2) (2m+3) (2m+4) (3m+2) 1 1 … 1 x x … x x x … x … … … … x x … x y y … y xy xy … xy xy xy … xy … … … … xy xy … xy n … … 4.4 Áp dụng. Ví dụ 1: Áp dụng phương pháp binh phương tối thiểu tìm đa thức bậc 2: P = a + ax + ax xấp xỉ với hàm cho bởi bảng sau: x 0.78 1.56 2.34 3.12 3.81 y 2.50 1.20 1.12 2.25 4.28 Ở đây m=2, n=5 và từ bảng 1 ta thu được bảng 2 để tính các hệ số của phương trình chuẩn. ( Quá trình tính toán thực hiện với 3 chữ số sau dấu phẩy). Bảng 2 x x x x x y xy xy 1 1 1 1 1 0.78 1.56 2.34 3.12 3.81 0.608 20434 5.476 9.734 14.516 0.475 3.796 12.813 30.371 55.306 0.370 5.922 29.982 94.759 210.717 2.50 1.20 1.12 2.25 4.28 1.950 1.872 2.621 7.020 16.307 1.520 2.921 6.133 21.902 62.128 5 11.61 32.768 102.761 341.750 11.35 29.770 94.604 Từ đó suy ra các hệ số : a, a , a của đa thức xấp xỉ P(x) cho từ hệ phương trình: (4 – 6) Giải hệ phương trình (4 – 6), ta có: a=5,045 ; a= - 4,043 ; a=1.009 Do đó đa thức xấp xỉ cần tìm có dạng P(x) = 5,045 - 4,043x + 1,009x Để so sánh các y với P và chuẩn bị tính sai số trung bình phương ta thực hiện tính toán trên bảng 3. Bảng 3 x Y P P-y [P-y] 0,78 1,56 2,34 3,12 3,81 2,50 1,20 1,12 2,25 4,28 2,505 1,194 1,110 2,252 4,288 0,005 -0,006 -0,010 0,002 0,008 0,000025 0,000036 0,000100 0,000004 0,000064 Suy ra =0,000229 Và sai số trung bình phương ==0,007 4.5 Chú ý Đối với trường hợp các điểm x cách đều nhau: x- x = h (i = 0, 1, … , n – 1) thì quá trình tính toán sẽ đơn giản hơn rất nhiều. Dưới đây ta sẽ trình bày kết quả trong trường hợp này. Trường hợp 1: Nếu n là số lẻ (n = 2k+1). Đặt u = hay x = x u.h Do đó khi x nhận các giá trị x, x, … , x , … , x thì u nhận các giá trị nguyên sau: -k, - k + 1, … , - 1, 0, 1, … , k – 1, k. Sau phép đổi biến (4 – 8) thì đa thức (4 – 1) cũng có bậc m và có dạng: Q= b+ bu + … + bu (4 – 9) Tương tự như (4 – 4) các hệ số b của (4 – 9) thu được từ hệ phương trình: (4 – 10) Hệ phương trình (4 – 10) so với hệ (4 – 4) đơn giản hơn rất nhiều vì các tổng những lũy thừa lẻ của u bằng 0: (4 – 11) Trường hợp 2: n chẵn (n = 2k) Ta đặt u = hoặc x = (4 – 12) Khi đó nếu x nhận các giá trị x, x, … , x thì u nhận các giá trị nguyên sau đây: -2k + 1, - 2k + 3, … , - 3, - 1, 1, 3, … , 2k – 3, 2k – 1 và trong hệ (4 – 10) cũng vắng mặt những tổng các lũy thừa lẻ của u: (4 – 13) Tóm lại, trong mội trường hợp (n lẻ hoặc n chẵn) vế trái của (4 – 10) đều vắng mặt các hệ số có dạng S (j là số lẻ). Ngoài ra các hệ số còn lại của vế trái (có dạng S, j chẵn) chỉ phụ thuộc vào n (vì unhận các giá trị nguyên). Do đó có thể lập những bảng tính sẵn các hệ số này (tùy thuộc vào n). Cuối cùng, sau việc giải phương trình (4 – 10) ta thu được Q dưới dạng (4 – 9). Để trở lại P dưới dạng (4 – 1) ta cần làm phép đổi biến ngược lại để chuyển biến u về biến x ban đầu. Cụ thể trong Q thu được ta sẽ dùng công thức đổi biến (4 – 8) nếu n lẻ, dùng công thức (4 – 12) nếu n chẵn. Dưới đây ta giải thích tinh thần của việc giải hệ (4 – 10) trong các trường hợp m = 1, m = 2. Trường hợp m = 1, nghĩa là (4 – 9) có dạng: Q(u) = b + bu. Đồng thời (4 – 10 ) có dạng (4 – 4) Từ đó suy ra: (4 – 14) Trường hợp m = 2 nghĩa là (4 – 9) có dạng Q(u) = b + bu + bu Và khi đó (4 – 10) có dạng: giải hệ 3 phương trình trên ta được: (4 – 15) Nếu ta gọi: (4 – 16) Khi đó các kết quả (4 – 14) và (4 – 15) có thể tóm tắt trong bảng 4. Ngoài ra từ (4 – 16) ta nhận thấy các số theo những giá trị lẻ của n từ 3 đến 21. Trong phần dưới của bảng 5 cho các số theo những giá trị chẵn của n từ 4 đến 22. m Các hệ số của Q(u) b b b 1 2 ( Để đơn giản trong phần này ta hiểu là ) Với n lẻ n 3 5 7 9 11 13 15 17 19 21 333333.10 200000.10 142857.10 111111.10 909091.10 769231.10 666667.10 588235.10 526316.10 476190.10 500000.10 100000.10 357143.10 166667.10 909091.10 549451.10 357143.10 245098.10 175439.10 129870.10 100000.10 485714.10 333333.10 255411.10 207459.10 174825.10 151131.10 133127.10 118973.10 107551.10 100000.10 142857.10 476190.10 216450.10 116550.10 699301.10 452489.10 309598.10 221141.10 163452.10 150000.10 714286.10 119048.10 324675.10 116550.10 499500.10 242405.10 128999.10 737137.10 445778.10 Với n chẵn n 4 6 8 10 12 14 16 18 20 22 250000.10 166667.10 125000.10 100000.10 833333.10 714286.10 625000.10 555556.10 500000.10 454545.10 500000.10 142857.10 595283.10 303030.10 174825.10 109890.10 735294.10 515996.10 375940.10 282326.10 640625.10 394531.10 289062.10 228906.10 189732.10 162109.10 141555.10 125651.10 112973.10 102628.10 781250.10 195312.10 781250.10 390625.10 223214.10 139509.10 930060.10 651042.10 473485.10 355114.10 156250.10 167411.10 372024.10 118371.10 468282.10 214629.10 109419.10 604683.10 356004.10 220567.10 V. XẤP XỈ HÀM TRONG THỰC NGHIỆM BẰNG ĐA THỨC TRỰC GIAO. 5.1 Định nghĩa hệ hàm trực giao. Xét hệ đa thức: R(x), R(x), …, R(x) (5 – 1) Trong đó R(x) = 1, R(x) = x + , R(x) = x+ x +, … Tổng quát R(x) = x + x + … + x + (5 – 2) Theo định nghĩa ta sẽ gọi (5 – 1) là hệ đa thức trực giao trên tập hợp X={x, x, …, x}, nếu (5 – 1) là hệ hàm trực giao trên tập X. Cụ thể là: (5 – 3) 5.2 Đặt vấn đề. Từ định nghĩa ta nhận thấy hệ đa thức trực giao là trường hợp đặc biệt của hệ hàm trực giao. Do đó ta áp dụng kết quả ở phần III với (x) = R(x) (r = 0, 1, …, m). Cụ thể là: khi cho u các giá trị thực nghiệm y(i = 1, 2, …, n) của hàm f(x) tại các điểm x (i = 1, 2, …, n) ta đặt vấn đề xấp xỉ hàm f(x) bởi một đa thức suy rộng cấp m (với hệ cơ sở (5 – 1)) có dạng: M(x) = . (5 – 4) Từ (3 – 14) ta suy ra các hệ số a của (5 – 4) có thể thu được từ công thức: a. (5 - 5) Từ (3 - 11) ta suy ra sai số trung bình phương của đa thức xấp xỉ M(x) là: (5 - 6) Ở đây M(x) (có dạng (5 – 4)) là một tổ hợp tuyến tínhcủa những đa thức đại số cấp từ 0 đến m, do đó M(x) thực chất cũng là một đa thức cấp m (như P(x) cho bởi (4 – 1)). Nghĩa là hàm xấp xỉ cũng là một đa thức đại số thông thường như đã thu được trong phần IV. Tuy nhiên do tính trực giao của hàm cơ sở (5 – 1) nên khác với phần IV ở đây ta không cần giải hệ phương trình chuẩn mà tìm các hệ số của đa thức (5 – 4) trực tiếp từ công thức (5 – 5) đã chỉ ra ở trên. Ngoài ra do những đặc điểm của hệ hàm trực giao ta có thể tăng dần cấp của M(x) mà không cần phải làm lại từ đầu quá trình tính toán. Đó chính là ưu điểm của phương pháp xấp xỉ hàm ở đây so với những kết quả thu được trong phần IV. 5.3 Nội dung của phương pháp. Nội dung chủ yếu của việc tìm đa thức xấp xỉ (5 – 4) thực chất là tìm hệ thức trực giao (5 – 1). Để làm được điều này ta tìm công thức truy hồi để xác định lần lượt các đa thức trực giao của hệ (5 – 1). Trước hết ta đi tìm những hàm đầu tiên: R(x), R(x) của hệ (5 – 1). Theo địng nghĩa thì R(x) = 1. Ngoài ra, từ (5 – 2) ta thấy R(x) có dạng: R(x) = x + Để xác định hệ số trong (5 – 8) ta sử dụng điều kiện đầu tiên trong (5 – 3) với r = 1 và s = 0: Từ đó suy ra: thay kết quả này vào (5 – 8) ta có: R(x) = x - (5 – 9) Để xác định những đa thức trực giao của hệ còn lại của hệ (5 – 1): R(x), R(x), …, R(x) ta sẽ chứng minh bổ đề sau đây: BỔ ĐỀ1: Mọi đa thức trực giao cấp r +1 (r 1): R(x) của hệ (5 – 1) được xác định theo các đa thức R(x) và R(x) từ công thức truy hồi sau: R(x) = (x + )R(x) + R(x) ( 5–10 ) Trong đó: CHỨNG MINH: Từ (5 – 3) ta có: (5 – 13 ) Từ (5 – 10) ta lại có == = (5 – 14) Nhưng Vậy từ (5 – 14) suy ra: (5 – 15) Kết hợp (5 – 13) và (5 – 15) ta có: hay (*) Từ (5 – 3) ta cũng có: (5 – 16) Từ (5 – 10) ta cũng có = = = (5 – 17) Nhưng Nên từ (5 – 17) suy ra: (5 – 18) Kết hợp (5 – 16) và (5 – 18) ta có: hay (* *) Từ (*) và (* *) thì bổ đề 1 được chứng minh hoàn toàn. Tuy nhiên để đơn giản các tử số và mẫu số của các công thức và ta sẽ chứng minh bổ đề sau: BỔ ĐỀ 2: Các tử và mẫu số của các công thức (5 – 11) và (5 – 12) có thể khai triển thành tổng những lũy thừa có dạng: S. (5 – 19) Cụ thể là: + = = (5 – 20) + (5 – 21) + (5 – 22) trong đó , , … , là các hệ số của đa thức R(x) cho dưới dạng (5 – 2). CHỨNG MINH: Từ (5 – 2) ta có thể viết lần lượt các đa thức R(x), R(x), … , R(x), …, R(x) dưới dạng: (5 – 23) Ta có x = R(x) - = R(x) - R(x) = R(x) + aR(x). x = R(x) - x - = R(x) - ( R(x) - R(x)) - = = R(x) - R(x) + ( - ) R(x) = = R(x) + a R(x) + a R(x). Tổng quát : x = R(x) + a R(x) + … + a R(x) + a R(x). Hay ta có thể thu được: (5 – 24) Với k < r và dựa trên (5 – 24) ta có: = = = = = = = + … + Từ đó dựa trên (5 – 3) ta có (với k < r): = 0 + .0 + …+ .0 = 0 (5 – 25) Giả sử (x) là đa thức bậc s < r : (x) = Khi đó: = = = Từ đó dựa trên (5 – 25) ta có: (5 – 26) Ngoài ra dựa trên (5 – 2) ta có: = = = Kết hợp với (5 – 25) ta có; = (5 – 27) Dựa trên (5 – 2) ta có: = = (5 – 28) Kết hợp (5 – 27) và (5 – 28 ) ta suy ra (5 – 20): Mặt khác vì R(x) là đa thức bậc r – 1 với hệ số của x là 1 nên x. R(x) là đa thức bậc r với hệ số của x là 1. Do đó ta có thể viêt: x. R(x) = x + (x) (5 – 29) Trong đó (x) là đa thức bậc r – 1. Từ (5 – 29) suy ra: = = (5 – 30) Nhưng dựa trên (5 – 26) ta có: = 0 Do đó từ (5 – 30) ta thu được (5 – 21): Cuối cùng từ (5 – 2) ta có: R(x) = x + (x) (5 – 31) Trong đó (x) là một đa thức cấp r – 2 Từ (5 – 31) suy ra xR(x) = x + x + (x). (5 – 32) Trong đó (x) là đa thức bậc r – 1. Từ (5 – 32) ta có x = { x + x + (x)} R(x). Do đó suy ra = = = (5 – 33) Nhưng dựa trên (5 – 26) ta có (5 – 34) Mặt khác dựa trên (5 – 2) ta có = = = (5 – 35) Thay (5 – 34), (5 – 35) vào (5 – 33) ta thu được (5 – 22) tức là Từ đó bổ đề 2 hoàn toàn được chứng minh. Từ bổ đề 1 và bổ đề 2 ta nhận thấy rằng: để thu được các đa thức trực giao của hệ (5 – 1), từ các công thức (5 – 10) và (5 – 12) ta cần tính tất cả các tổng những lũy thừa có dạng S. ( = 1, 2, … , 2m – 1) Ngoài ra khi áp dụng công thức (5 – 5) để tìm các hệ số a, a, …, a của (5 – 4) lại cần tính các tổng , …, ở mẫu số của công thức. Nghĩa là dựa trên (5 – 20) cần tính các tổng những lũy thừa: S ( = 1, 2,… , 2m) Còn các tử số của công thức (5 – 5) lần lượt là [y,R], [y,R], …, [y,R]. Để tính mỗi tử số này ta dựa vào (5 – 2) và dựa vào khai triển (r = 0, 1, …, m) = = (5 – 36) Tóm lại: Để tìm hàm xấp xỉ M(x) ta cần tìm 2m tổng S( = 1, 2, …, 2m) và m tổng (r = 1, 2, …, m). Khi đó trong bảng tính toán giải bài toán mỗi tổng nói trên được lập theo một cột. 5.4 Sai số của phương pháp. Cuối cùng để tính sai số trung bình phương một cách thuận lợi ta dùng công thức (5 – 6) nên khi tính toán ta cần tính thêm tổng . Khi sai số trung bình phương tìm được chưa đủ bé (nghĩa là m chưa đủ lớn) ta cần tăng dần cấp m của hàm xấp xỉ M(x). Khi đó trong bảng tính cũ cần bổ xung những cột tính S và mới nhưng kết quả cũ vẫn được sử dụng. THIẾT KẾ CHƯƠNG TRÌNH Giao diện chương trình: Cho phép người sử dụng lựa chọn các mục chính trong chương trình để nhập các ssố liệu cần thiết. Sauđó ta chọn mục “Nhập hàm số dưới dạng bảng” và ta nhập dữ liệu thì kết quả thu được như sau: Bảng tính hệ số tức là các thông số trung gian để xây dựng đa thức: Kết quả thu được như sau: CODE CỦA CHƯƠNG TRÌNH #include #include #include #include int n,m; double x[50],y[50],g[50]; double xx[50],yy[50]; double A[50][50],B[50],a[50]; FILE *f1; int t; double e; char xau[256]; /*--------------------------------------------*/ void nhap_ham(void) { clrscr(); int i; f1=fopen("kq","w"); printf("Nhap so mau quan sat n= "); scanf("%d",&n); fprintf(f1,"So mau quan sat n= %d\n",n); /*Xoa bo nho dem ban phim*/ fflush(stdin); /*Nhap ham so duoi dang bang*/ printf("Nhap so chu so phan thap phan t= "); scanf("%d",&t); printf("\nDoc cac diem quan sat x[i] (i=1..n)\n"); for (i=1;i<=n;i++) { printf("x[%d]= ",i); scanf("%lf",&x[i]); } printf("Doc ket qua quan sat y[i] (i=1..n)\n"); for (i=1;i<=n;i++) { printf("y[%d]= ",i); scanf("%lf",&y[i]); } /*In ra man hinh ca mau quan sat*/ printf("\nCac mau quan sat la:\n"); printf("X="); fprintf(f1,"\nCac mau quan sat la:\n"); fprintf(f1,"X="); for (i=1;i<=n;i++) { printf("%10.*lf",t,x[i]); fprintf(f1,"%10.*lf",t,x[i]); } printf("\n"); printf("Y="); fprintf(f1,"\nY="); for (i=1;i<=n;i++) { printf("%10.*lf",t,y[i]); fprintf(f1,"%10.*lf",t,y[i]); } printf("\n"); fprintf(f1,"\n"); fclose(f1); } /*----------------------------------------*/ void bac_da_thuc(void) { int i; clrscr(); f1=fopen("kq","a"); printf("Nhap vao bac da thuc xap xi m= "); scanf("%d",&m); printf("\nHam so se xap xi da thuc co dang:\n"); fprintf(f1,"\nHam so xap xi da thuc co dang:\n"); for (i=0;i<=m;i++) { printf(" a%d*x^%d",i,m-i); if (i>=m) continue; printf("+"); fprintf(f1,"+"); } fclose(f1); } /*------------------------------------------*/ /*Xay dung ham mu de lap bang tinh he so*/ double mu(double c,int k) { double tg; int i; if (k == 0) tg = 1; if (k == 1) tg = c; if (k > 1) { tg = 1; for (i=1;i<=k;i++) tg = tg*c; } return(tg); } /*Lap bang tinh he so*/ void bang_he_so(void) { clrscr(); int i,j; f1=fopen("kq","a"); xx[0]=n; for (j=1;j<=2*m;j++) xx[j]=0; for (j=1;j<=2*m;j++) { for(i=1;i<=n;i++) xx[j]=xx[j]+mu(x[i],j); } yy[0]=0; for (i=1;i<=n;i++) yy[0]=yy[0]+y[i]; for (j=1;j<=m;j++) yy[j]=0; for (j=1;j<=m;j++) { for (i=1;i<=n;i++) yy[j]=yy[j]+mu(x[i],j)*y[i]; } printf("Bang tinh he so gom:\n"); fprintf(f1,"\n\nBang tinh he so gom:\n"); printf("\n1.Tong cac x^i (i=0..2*m)"); fprintf(f1,"\n1.Tong cac x^i (i=0..2*m)"); for (j=0;j<=2*m;j++) { printf("\nxx[%d]= %10.*lf",j,t,xx[j]); fprintf(f1,"\nxx[%d]= %10.*lf",j,t,xx[j]); } printf("\n\n2.Tong cac (x^i)*y (i=0..m)"); fprintf(f1,"\n\n2.Tong cac (x^i)*y (i=0..m)"); for (j=0;j<=m;j++) { printf("\nyy[%d]= %10.*lf",j,t,yy[j]); fprintf(f1,"\nyy[%d]= %10.*lf",j,t,yy[j]); } fclose(f1); /*Gan bang he so cho ma tran he phuong trinh*/ for (i=0;i<=m;i++) for (j=0;j<=m;j++) A[i][j]=xx[i+j]; for (i=0;i<=m;i++) B[i]=yy[i]; } /*------------------------------------------*/ void gauss(void) { int i=0,done=0; int j,p,q; double s,s1; double max,c; g[i]=0; f1=fopen("kq","a"); /*Hien thi he phuong trinh duoi dang ma tran*/ printf("He phuong trinh duoi dang ma tran "); fprintf(f1,"\n\nHe phuong trinh duoi dang ma tran"); for (q=0;q<=m;q++) { printf("\n"); fprintf(f1,"\n"); for (j=0;j<=m;j++) { printf("%12.*lf",t,A[q][j]); fprintf(f1,"%12.*lf",t,A[q][j]); } printf(" =%12.*lf",t,B[q]); fprintf(f1," =%12.*lf",t,B[q]); } printf("\n\nGiai he phuong trinh bang phuong phap Gauss\n"); fprintf(f1,"\n\nGiai he phuong trinh bang phuong phap Gauss\n"); /*--------------------------------------------------------*/ while (!done) { if (A[i][i]==0) { max=0; p=i; for (q=i+1;q<=m;q++) if (max < fabs(A[q][i])) { p=q; max=fabs(A[q][i]); } if (p !=i) { for (j=i;j<=m;j++) { c=A[i][j]; A[i][j]=A[p][j]; A[p][j]=c; } c=B[i]; B[i]=B[p]; B[p]=c; } if (p==i) done =1; } if (A[i][i] !=0) { c=1/A[i][i]; for (j=i;j<=m;j++) A[i][j]=A[i][j]*c; B[i]=B[i]*c; for (q=i+1;q<=m;q++) { for (j=i+1;j<=m;j++) A[q][j]=A[q][j]-A[i][j]*A[q][i]; B[q]=B[q]-B[i]*A[q][i]; A[q][i]=0; } } i++; if (i>m) done=1; } if (i>m) { a[m]=B[m]/A[m][m]; for (i=m-1;i>=0;i--) { a[i]=0; for (j=m;j>i;j--) B[i]=B[i]-A[i][j]*a[j]; a[i]=B[i]/A[i][j]; } } printf("He so da thuc la:\n"); fprintf(f1,"\nHe so da thuc la:\n"); for (i=0;i<=m;i++) { printf("a[%d]= %12.*lf\n",i,t,a[i]); fprintf(f1,"%12.*lf\n",t,a[i]); } /*Hien dang da thuc thu duoc*/ printf("\nDa thuc thu duoc la:\n\n"); fprintf(f1,"\nDa thuc thu duoc la:\n\n"); printf("\n "); for(i=0;i<=m;i++) { cprintf("%8.*f*x^%d",t,a[i],m-i); fprintf(f1,"%8.*f*x^%d",t,a[i],m-i); if(i>=m) continue; if (a[i]>0) { printf(" "); fprintf(f1," "); } else printf("+ "); fprintf(f1," +"); } printf("\n"); fprintf(f1,"\n"); for(i=0;i<=m;i++) for(j=0;j<=m;j++) { g[i]=g[i]+a[i]*mu(x[j],m-i); } printf("g[%d]=%12.*lf\n",i,t,g[i]); s1=g[i]; for(i=0;i<=m;i++) { s1=s1+ g[i]; } printf("%12.lf\n",s1); s=s1; for (i=0;i<=m;i++) { s=s-y[i]; } s=sqrt((1/n)*s); printf("%12.lf\n",s); for(i=0;i<=m;i++) { for (j=1;j<=m;j++) { printf("g[%d]=%12.*lf\n",i,t,g[i]); fprintf(f1,"%12.*lf",t,g[i]); } } fclose(f1); } /*-----------------------------------------*/ void luu_tep(void) { char ten_tep[40]; int i=1; fflush(stdin); printf("Ban hay go vao ten tep in ket qua (nho them duoi .txt).\n "); printf("Ten tep la: "); gets(ten_tep); unlink(ten_tep); f1=fopen("kq","r"); rename("kq",ten_tep); printf("%s","Noi dung cua file"); while (!feof(f1) ) { fgets(xau,256,f1); i++; if (i==20) { clrscr(); i=1; } gotoxy(2,i); printf("%s",xau); getch(); } fclose(f1); } void hienkt(char *dkt,int h,int c,int dai,int mn,int mc) { textbackground(mn); textcolor(mc); window(c,h,c+dai-1,h); clrscr(); /*Cat bot ki tu neu can*/ if (strlen(dkt) >= dai) dkt[dai] =0; puts(dkt); gotoxy(1,1); } /*Ket thuc ham hienkt*/ void menu(int n,int hd,int cd,int dai,char nd[][40],int mnc, int mcc,int mnr,int mcr,int *ch) { int i,chonm,mn,mc; int ch1,ch2,ch4,chon; chon=*ch; /*Chinh li bien chon va kiem tra*/ if (n= 20) return; if (chon n) chon =1; /*Trinh bay*/ chonm=chon; for (i=1;i<=n;++i) hienkt(nd[i],hd+i-1,cd,dai,mnc,mcc); hienkt(nd[chonm],hd+chonm-1,cd,dai,mnr,mcr); /*Bat phim chon chuc nang*/ ch1=0; while (ch1 !=13) { ch1=getch(); if (ch1 == 0) ch2=getch(); if (ch1 == 0) if (ch2 == 72) { chonm=chon-1; if (chonm <=0 ) chonm=n; } else if (ch2 == 80) { chonm=chon+1; if (chonm>n) chonm=1; } if (chonm != chon) { hienkt(nd[chon],hd+chon-1,cd,dai,mnc,mcc); hienkt(nd[chonm],hd+chonm-1,cd,dai,mnr,mcr); chon=chonm; } } /*Ket thuc vong lap while*/ *ch = chon; } /*Ket thtuc ham menudh*/ void main() { char nd[7][40]; int chon; strcpy(nd[1],"1.Nhap ham so duoi dang bang"); strcpy(nd[2],"2.Nhap dang cua da thuc xap xi"); strcpy(nd[3],"3.Lap bang tinh he so"); strcpy(nd[4],"4.Ket qua"); strcpy(nd[5],"5.Hien thi tep"); strcpy(nd[6],"6.Ket thuc"); textmode(C80); clrscr(); chon = 1; while (chon !=6) { textbackground(BLUE); clrscr(); hienkt(" Toan tin 2 K48",21,25,18,BLUE,WHITE); menu(6,9,25,30,nd,CYAN,YELLOW,BROWN,WHITE,&chon); if (chon == 1) { window(1,1,80,25); textbackground(BLACK); clrscr(); nhap_ham(); getch(); } if (chon == 2) { window(1,1,80,25); textbackground(BLACK); clrscr(); bac_da_thuc(); getch(); } if (chon == 3) { window(1,1,80,25); textbackground(BLACK); clrscr(); bang_he_so(); getch(); } if (chon == 4) { window(1,1,80,25); textbackground(BLACK); clrscr(); gauss(); getch(); } if (chon == 5) { window(1,1,80,25); textbackground(BLACK); clrscr(); luu_tep(); getch(); } } } KẾT LUẬN Phương pháp bình phương tối thiểu lập công thức từ thực nghiệm là một bài toán cơ bản hay gặp trong toán học cũng như trong thực tế đời sống. Thông qua phương pháp này ta có thể tìm ra được gần như chính xác đa thức ban đầu. Thông qua chương trình cụ thể viết trên ngôn ngữ lập trinh C thì ta có thể thấy phần nào tính ưu việt của phương pháp này. Tuy nhiên vì thời gian có hạn nên bài toán mà em tìm hiểu không tránh khỏi những thiếu sót. Em rất mong sự chỉ bảo, đóng góp ý kiến của thầy cô và các bạn để chương trình của em hoàn thiện hơn. Em xin chân thành cảm ơn PGS – TS Lê Trọng Vinh đã nhiệt tình chỉ bảo và giúp đỡ em trong quá trình làm đồ án môn học này. Hà Nội ngày 16 tháng 6 năm 2007 Sinh viên : Bùi Văn Bằng TÀI LIỆU THAM KHẢO Cơ sở phương pháp tính. Tác giả: Phan Văn Hạp – Nguyễn Quý Hỷ - Hồ Thuần – Nguyễn Công Thúy. Nhà xuất bản Đại học và trung học chuyên nghiệp. Giáo trình phương pháp tính. Tác giả: Dương Thủy Vỹ. Nhà xuất bản khoa học kỹ thuật – 2005. Giải tích số. Tác giả: Lê Trọng Vinh. Nhà xuất bản khoa học và kỹ thuật – 2000. Giải tích số. Tác giả: Phạm Kỳ Anh. Nhà xuất bản đại học quốc gia – 1996. Phương pháp tinh. Tác giả: Tạ Văn Đĩnh – Lê Trọng Vinh. Nhà xuất bản đại học và trung học chuyên nghiệp – 1983.