Bài giảng Phương pháp số - Bài 4: Nội suy bằng đa thức và làm khớp dữ liệu - Nguyễn Thị Vinh

BÀI 4 NỘI SUY BẰNG ĐA THỨC VÀ LÀM KHỚP DỮ LIỆU 0 1 2 3 4 5 6 0 2 4 6 8 10 12 14 NỘI SUY BẰNG ĐA THỨC (1) BÀI TOÁN NỘI SUY: Cho x0, x1, , xn là n + 1 điểm phân biệt trên trục số thực và f(x) là hàm nhận giá trị thực, xác định trên khoảng I = [a, b] chứa các điểm này. Hãy xây dựng một đa thức pn(x) có bậc ≤ n mà tại các điểm x0, x1, , xn pn(xi) = f(xi) i = 0, , n SỰ TỒN TẠI ĐA THỨC NỘI SUY: (Đa thức nội suy Lagrange) Cho hàm f(x) nhận giá trị thực và n + 1 điểm phân biệt x0, x1, , xn, khi đó tồn tại đúng một đa thức bậc ≤ n nôi suy f(x) tại x0, x1, , xn là pn(x) = a0l0(x) + a1l1(x) + + anln(x) với ai = f(xi) và ik i n ki 0i k xx xx (x)l Π      PHƢƠNG PHÁP SỐ - Bài 4 2 NỘI SUY BẰNG ĐA THỨC (2) SỰ DUY NHẤT CỦA ĐA THỨC NỘI SUY Bổ đề: Nếu z1, , zn là các nghiệm phân biệt của đa thức p(x) thì p(x) = (x – z1)(x – z2) (x – zn) r(x) với r(x) là một đa thức. Hệ quả: Nếu p(x) và q(x) là hai đa thức bậc ≤ k có giá trị trùng nhau ở k+1 điểm z0, , zk phân biêt, thì p(x) ≡ q(x)  có nhiều nhất một đa thức bậc ≤ n nội suy f(x) ở n + 1 điểm phân biệt x0, x1, , xn Mặt khác, từ sự tồn tại của đa thức nội suy Lagrange  có ít nhất một đa thức bậc ≤ n nội suy f(x) ở n + 1 điểm phân biệt x0, x1, , xn  Kết luận: có đúng một đa thức bậc ≤ n nội suy f(x) ở n + 1 điểm phân biệt x0, x1, , xn PHƢƠNG PHÁP SỐ - Bài 4 3 NỘI SUY BẰNG ĐA THỨC (3) Ví dụ 1: trường hợp n = 1, nghĩa là cho biết hàm f(x) và hai điểm phân biệt x0, x1. Vậy ta có hai đa thức bậc nhất  Trong ví dụ này, đa thức nội suy Lagrange là đa thức nội suy tuyến tính (n = 1) 01 0 1 10 1 0 xx xx (x)l xx xx (x)l       ) )) 0 01 01 0 10 0110 01 0 1 10 1 01100n x(x xx )f(x)f(x )f(x xx x)(xf(xx)(xf(x xx xx )f(x xx xx )f(x(x))lf(x(x))lf(x(x)p              PHƢƠNG PHÁP SỐ - Bài 4 4 Đây chính là PT đƣờng thẳng đi qua 2 điểm (x0, y0) và (x1, y1) NỘI SUY BẰNG ĐA THỨC (4) Ví dụ 2: Từ bảng các giá trị của tích phân sau, tính giá trị của đa thức nội suy Lagrange K(3.5), biết K(1) = 1.5709, K(4) = 1.5727, và K(6) = 1.5751 Ta có  K(3.5) ≈ (1.5709)(0.08333)+(1.5727)(1.04167)+(1.5751)((–0.12500) = 1.57225    π/2 0 1/222 x]sink)(sin [1 dx K(k) 0 1 2 (3.5 4)(3.5 6) 1.25 l (3.5) 0.08333 (1 4)(1 6) 15 (3.5 1)(3.5 6) 6.25 l (3.5) 1.04167 (4 1)(4 6) 6 (3.5 1)(3.5 4) 1.25 l (3.5) 0.12500 (6 1)(6 4) 10                          PHƢƠNG PHÁP SỐ - Bài 4 5 NỘI SUY BẰNG ĐA THỨC (5) NHƯỢC ĐIỂM CỦA ĐA THỨC NỘI SUY Lagrange: • Việc ƣớc lƣợng nó ở một điểm x cần ít nhất 2(n+1) phép nhân/chia và (2n+1) phép cộng và trừ sau khi tất cả các mẫu số của các đa thức nội suy Lagrange đã đƣợc tính toán và chia vào trong các giá trị hàm tƣơng ứng • Việc dùng đa thức Lagrange dƣờng nhƣ là lãng phí, bởi vì khi tính pk(x) không dùng lại đƣợc các kết quả đã tính ở pk–1(x)  Nên tìm đa thức nội suy dạng khác PHƢƠNG PHÁP SỐ - Bài 4 6 ĐA THỨC NỘI SUY NEWTON (1) Lập công thức pn(x) = A0 + A1 (x – x0) + A2 (x – x0) (x – x1) + + An(x – x0)(x – x1) (x – xn–1) = qk(x) + (x – x0)(x – x1) (x – xk)r(x) với qk(x) = A0 + A1 (x – x0) + A2 (x – x0) (x – x1) + + Ak(x – x0)(x – x1) (x – xk–1) Nhận xét: pn(x) = qk(x) tại các điểm x0, x1, , xk, nghĩa là qk(x) nội suy f(x) tại x0, , xk PHƢƠNG PHÁP SỐ - Bài 4 7 ĐA THỨC NỘI SUY NEWTON (2) Lập công thức: Xây dựng dãy đa thức có tính kế thừa p0(x),p1(x), pk(x) = pk–1(x) + Ak(x – x0)(x – x1) (x – xk–1) Ak = f[x0, , xk] là hệ số của lũy thừa bậc cao nhất xk của pk(x), chỉ phụ thuộc vào các giá trị của f(x) tại các điểm, x0, , xk đƣợc gọi là tỉ sai phân cấp k của f(x) tại x0, , xk PHƢƠNG PHÁP SỐ - Bài 4 8 pn(x) = f[x0] + f[x0, x1](x – x0) + f[x0, x1, x2](x – x0)(x – x1) + + f[x0, x1, , xn](x – x0)(x – x1) (x – xn–1) ĐA THỨC NỘI SUY NEWTON (3) Công thức tính tỉ sai phân Với n = 1, công thức nội suy Newton sẽ là p1(x) = f[x0] + f[x0, x1](x – x0) so sánh với công thức nội suy Lagrange p1(x) = f(x0) + [f(x1) – f(x0)](x – x0) / [x1 – x0]  f[x0] = f(x0) f[x0, x1] = . f[x0, x1, , xk] = 01 01 01 01 xx f[xf[x xx )f(x)f(x      ]] 1 k 0 k 1 k 0 f[x ,...,x ] f[x ,...,x ] x x   PHƢƠNG PHÁP SỐ - Bài 4 9 ĐA THỨC NỘI SUY NEWTON (3) Bảng tỉ sai phân xi f[ ]=f( ) f[,] f[, ,] f[, , ,] f[, , , ,] x0 f(x0) f[x0, x1] x1 f(x1) f[x0, x1, x2] f[x1, x2] f[x0, x1, x2, x3] x2 f(x2) f[x1, x2, x3] f[x0, x1, x2, x3, x4] f[x2, x3] f[x1, x2, x3, x4] x3 f(x3) f[x2, x3, x4] f[x3, x4] x4 f(x4) PHƢƠNG PHÁP SỐ - Bài 4 10 ĐA THỨC NỘI SUY NEWTON (4) Ví dụ: Từ ví dụ 2 về đa thức nội suy Lagrange  p2(x) = 1.5709 + 0.0006 (x – 1) + 0.00012 (x – 1) (x – 4) thế x = 3.5 vào công thức nội suy Newton, ta nhận đƣợc p2(x) = 1.5709 + 0.0006 (2.5) + 0.00012 (2.5) (–0.5) = 1.57225 xi K[ ]=f( ) K[,] K[, ,] 1 1.5709 0.0006 4 1.5727 0.00012 0.0012 6 1.5751 0.0012 64 1.57511.5727 K(4,6) 0.0006 41 1.57271.5709 K(1,4)         0.00012 61 0..00120.0006 K(1,4,6)     PHƢƠNG PHÁP SỐ - Bài 4 11 ĐA THỨC NỘI SUY NEWTON (5) • Thuật toán tính các hệ số của đa thức nội suy Newton Cho n + 1 điểm phân biệt x0, x1, , xn , và f(x0), f(x1), , f(xn) tƣơng ứng, với f(xi) đƣợc chứa trong vector d = (di), i =0, ..., n Với k = 1, , n, lặp công việc (tính d ở lần lặp k) Với i = 0, , n – k, lặp công việc: di = (di + 1 – di )/ (xi+k – xi) Vậy f[x0, , xk] = d0, ở lần lặp thứ k • Sai số của đa thức nội suy pn(x) trên [a, b] chứa x0,x1, ,xn: PHƢƠNG PHÁP SỐ - Bài 4 12 :)b,a()x(ξξ],b,a[x          n 0j j )1n( j n 0j n0 nn )xx( )!1n( )ξ(f )x(x ] x, x..., ,x[f )x(p)x(f)x(e Π ĐA THỨC NỘI SUY NEWTON (6) Chƣơng trình: double noiSuyNewton(vector x, vector y, double x0) { vector d = y; double tong = d[0], tich = 1; for (unsigned k = 1; k <= d.size() – 1; k++) { for (unsigned i = 0; i <= d.size() – k; i++) d[i] = (d[i+1] – d[i])/(x[i+k] – x[i]); tich *= (x0 – x[k – 1]); tong += d[0] * tich; } return tong; } PHƢƠNG PHÁP SỐ - Bài 4 13 ĐA THỨC NỘI SUY NEWTON VỚI MỐC CÁCH ĐỀU (1) Giả sử f(x) đƣợc lập thành bảng đối với xi = a + ih với các giá trị f(xi), i = 0, . . . , N đã biết, dùng phép đổi biến tuyến tính đƣa đa thức bậc n của x về đa thức bậc n của s theo các tâm nguyên. Định nghĩa các sai phân tiến s00 0 fsh)f(x f(x) và shx x(s)x h xx s(x)s             0i fΔfΔ)fΔ(Δ 0i f fΔ s 1i 1s 1i s 1i s s i PHƢƠNG PHÁP SỐ - Bài 4 14 0 i fΔ hi! 1 ]x ..., ,f[x k i iikk  ĐA THỨC NỘI SUY NEWTON VỚI MỐC CÁCH ĐỀU (2) Đa thức Newton bậc ≤ n nội suy f(x) tại xk, . . . , xk+n trở thành Do (x – xk+j) = x0 + sh – [x0 + (k + j) h] = (s – k – j) h nên ta có công thức sai phân tiến ... 1, 0, k ),x(xfΔ hi! 1 (x)p jk 1i 0j k n 0i i in Π       k n k 2 kk 1i 0j k n 0i i 0nn fΔ n! 1) n - k - (s ... k) - (s ... fΔ 2 1)- k - k)(s - (s fk) - (s f 1j jks fΔsh)(xp(x)p Π           PHƢƠNG PHÁP SỐ - Bài 4 15 ĐA THỨC NỘI SUY NEWTON VỚI MỐC CÁCH ĐỀU (3) x k f k ∆fk ∆ 2fk ∆ 3fk ∆ 4fk x 0 f 0 △f 0 x 1 f 1 △ 2f 0 △f 1 △3f 0 x 2 f 2 △ 2f 1 △4f 0 △f 2 △3f 1 x 3 f 3 △ 2f 2 △f 3 x 4 f 4 PHƢƠNG PHÁP SỐ - Bài 4 16 Bảng sai phân tiến ĐA THỨC NỘI SUY NEWTON VỚI MỐC CÁCH ĐỀU (4) Công thức nội suy Newton tiến Phép nội suy hàm f(x) tại x0, . . . , xn khi k = 0, trở thành Các hệ số ∆if0 ở công thức trên chỉ đơn thuần là hiệu số giữa đầu vào phía dƣới bên trái và đầu vào phía trên bên trái. 0 n 0 2 00 0n0 fΔ n! 1) n - (s ... 1) - s(s ... fΔ 2! 1) - (s s fs f sh)(xpsh)f(x f(s)     PHƢƠNG PHÁP SỐ - Bài 4 17 ĐA THỨC NỘI SUY NEWTON VỚI MỐC CÁCH ĐỀU (5) xk fk ∆ fk ∆ 2 fk ∆3 fk ∆4 fk -4 -64 56 -2 -8 -48 8 48 0 0 0 0 8 48 2 8 48 56 4 64 Ví dụ: Cho các giá trị của f(x) tại các điểm cách đều như bảng bên. Tính f(-3). Do x = - 3 ở đầu bảng giá tri xk nên ta chọn x0 = -4.  s = [-3 - ( -4 )]/2 = 0,5 f(-3) = -64 + 56.0,5 -48.0,5.(0,5-1)/2! + 48. 0,5.(0,5-1)(0,5-2)/3! = -27 PHƢƠNG PHÁP SỐ - Bài 4 18 double noiSuyNewtonTien (vector x, vector y, double x0) { unsigned i, k, n = y.size(); vector d = y; double h = x[1] – x[0], s = (x0 – x[0]) / h, mauSo = 1. ; double f = y[0], tich = s; for (k = 1; k <= n – 1; k++) { // Tinh tong thu k for (i = 0; i <= n – k; i++) d[i] = d[i+1] – d[i];// cac sf cap k f += tich * d[0] * mauSo; tich *= (s – k); mauSo /= (k+1); } return f; } PHƢƠNG PHÁP SỐ - Bài 4 19 ĐA THỨC NỘI SUY NEWTON VỚI MỐC CÁCH ĐỀU (6) Công thức nội suy Newton lùi Khi điểm cần tính nội suy x ở cuối bảng, ta lập đa thức nội suy pn(x) = A0 + A1 (x – xn) + A2 (x – xn) (x – xn-1) + + An(x – xn)(x – xn-1) (x – x0) Tƣơng tự nhƣ khi tính toán các hệ số Ai cho công thức nội suy ở đầu bảng, dẫn đến công thức pn(x) = f[xn] + f[xn-1, xn](x – xn) + f[xn-2, xn-1, xn](x – xn)(x – xn-1) + + f[x0, x1, , xn](x – xn)(x – xn-1) (x – x0) Khi các mốc nội suy cách đều nhau một khoảng h, dùng phép biến đổi tuyến tính snn n fsh)f(x f(x) và shx x(s)x h xx s(x)s    PHƢƠNG PHÁP SỐ - Bài 4 20 ĐA THỨC NỘI SUY NEWTON VỚI MỐC CÁCH ĐỀU (7) Công thức nội suy Newton lùi Việc nội suy hàm f(x) tại xn, . . . , x0 trở thành Các hệ số ∆i fn–i ở công thức trên là các sai phân nằm ở cạnh bên dƣới của bảng tam giác sai phân 0 n 2-n 2 1-nn nnn fΔ n! 1) n (s ... 1) s(s ... fΔ 2! 1) (s s fs f sh)(xpsh)f(x f(s)      )x(xfΔ hi! 1 (x)p jn 1i 0j i-n n 0i i in Π      PHƢƠNG PHÁP SỐ - Bài 4 21 ĐA THỨC NỘI SUY NEWTON VỚI MỐC CÁCH ĐỀU (8) xk fk ∆ fk ∆ 2 fk ∆3 fk ∆4 fk -4 -64 56 -2 -8 -48 8 48 0 0 0 0 8 48 2 8 48 56 4 64 Ví dụ: Cho bảng sai phân của f(x) tại các điểm cách đều. Tính f(3). Do x = 3 ở cuối bảng giá tri xk nên ta chọn xn = 4.  s = [3 - 4 ]/2 = -0,5  f(3) = 64 + 56.(-0,5) + 48.(–0,5) (–0,5+1)/2! + 48.(–0,5)(–0,5+1)(–0,5+2)/3! = 27 PHƢƠNG PHÁP SỐ - Bài 4 22 ĐA THỨC NỘI SUY NEWTON VỚI MỐC CÁCH ĐỀU (9) CÁC ĐA THỨC GHÉP TRƠN (1) 1. ĐẶT VẤN ĐỀ • Biểu thức tính sai số của đa thức nội suy Newton gợi ý rằng việc tăng bậc của đa thức nội suy sẽ cung cấp một xấp xỉ tốt hơn. Nhƣng, sai số phụ thuộc vào độ lớn của khoảng chứa các điểm nội suy, làm cho thủ thuật này không thành công PHƢƠNG PHÁP SỐ - Bài 4 23 C0 5 7.5 9.9 12.8 13.2 15.1 16.3 16.8 D 0.0240 0.0437 0.0797 0.1710 0.1990 0.326 0.8460 0.9720 • Xét ví dụ sau: Dữ liệu thực nghiệm về hệ số khuyếch tán D đối với chất dung môi C0 đƣợc cho trong bảng dƣới đây 4 6 8 10 12 14 16 18 -2 -1.5 -1 -0.5 0 0.5 1 Noi suy bang da thuc Newton PHƢƠNG PHÁP SỐ - Bài 4 24 P (x) có nhiều giá trị âm không phù hợp thực tế. Ví du với x0 = 6 thì y0 = -1.565 CÁC ĐA THỨC GHÉP TRƠN (2) 1. ĐẶT VẤN ĐỀ • Tìm các đa thức bậc ba nội suy f(x) từng khúc, nối 2 điểim liên tiếp xi, xi+1, và ghép trơn với nhau, tức là chúng trơn ở các điểm nối Mi(xi,yi) đã cho • Giả sử có hàm f(x) xác định trên đoạn [x0,xN] và biết giá trị của nó tại các điểm x0 < x2 < < xN là f(x0), f(x1) , f(xN) tƣơng ứng. Hãy tìm N đa thức bậc ba ghép trơn trên [xi,xi+1] nội suy f(x) tại các điểm này. Pi(x) =c1,i + c2,i(x–xi) + c3,i(x–xi) 2 + c4,i(x–xi) 3, i = 0, 1, , N-1 PHƢƠNG PHÁP SỐ - Bài 4 25 PHƢƠNG PHÁP SỐ - Bài 4 26 4 6 8 10 12 14 16 18 -2 -1.5 -1 -0.5 0 0.5 1 Noi suy Newton & Noi suy Spline bac 3 Các đa thức ghép trơn Spline bậc 3 phù hợp thực tế. Ví du với x0 = 6 thì y0 = 0.0275 CÁC ĐA THỨC GHÉP TRƠN (3) 2. ĐA THỨC NỘI SUY HERMIT: Do các đa thức Pi(x)nội suy f(x) trên [xi xi+1] nên ta có Pi (xi) = f(xi), Pi (xi+l) = f(xi+1) i = 0, . . . , N-1 và Pi -1 (xi) = Pi(xi), i = 1, . . . , N-1 Các đa thức này trơn tại các điểm trong xi, tức là Pi´(xi) = f´(xi) i = 1, . . . , N-1 (1) Pi´(xi+l) = f´(xi+l) i = 1, . . . , N-1 (2) Từ công thức nội suy Newton bậc ba tại 4 điểm xi, xi, xi+1, xi+1 Pi(x) = f(xi) + f[xi, xi](x – xi) + f[xi, xi, xi+1](x – xi) 2 + f [xi,, xi, xi+1, xi+1](x – xi) 2(x – xi+1) Bởi vì (x – xi+1) = (x – xi) + (xi – xi+1) = (x – xi) - hi, hi = xi+1 – xi nên Pi(x) = f(xi) + f ’(xi) (x – xi) + (f [xi, xi, xi+1] – f [xi, xi, xi+1, xi+1] hi) (x – xi) 2 + f[xi, xi, xi+l, xi+l] (x – xi) 3 PHƢƠNG PHÁP SỐ - Bài 4 27 CÁC ĐA THỨC GHÉP TRƠN (4) 2. ĐA THỨC NỘI SUY HERMIT: SAI SỐ Với x ∈ [xi, xi+1], g3(x) = Pi(x), ở đây Pi (x) nội suy f(x) tại 4 điểm xi, xi, xi+1, xi+1, từ công thức sai số của đa thức nội suy Newton, suy ra với x ∈ [xi, xi+1], f(x) – g3(x) = f[xi, xi, xi+1, xi+l, x] (x – xi) 2 (x – xi+1) 2 với giả thiết f(x) khả vi liên tục cấp bốn. Hơn nữa,     16 4 2 2 12 2 12 1 2 1 )( max iiiii ],xi[xx h hh|)x(x)x|(x i               |)maxmax ! |)()(| (ξ|f)(hxgxf )( kxξkx k k 44 3 1 16 1 4 1 4 i )4( [a,b]ξ 3 384 )hmax( |)(ξf|max|)x(g)x(f| j   PHƢƠNG PHÁP SỐ - Bài 4 28 nên với a ≤ x ≤ b: CÁC ĐA THỨC GHÉP TRƠN (3) 2. ĐA THỨC NỘI SUY HERMIT: Kí hiệu fi = f(xi), si = f’(xi), i = 0 . . . , N ta nhận đƣợc Pi(x) =c1,i + c2,i(x–xi) + c3,i(x–xi) 2 + c4,i(x–xi) 3 với c1,i = fi c2,i = si c3,i = f[xi, xi, xi+1] – f[xi, xi, xi+1, xi+1] hi i,i4 i i1ii hc h s],xf[x     i 1iii1i1ii ,i4 h ],x,xf[x],x,xf[x c    2 i 1iii1i )(h ],xf[x2ss   i1i i1ii1ii1i i hh ],xf[xh],xf[xh s      PHƢƠNG PHÁP SỐ - Bài 4 29 Khi không biết thông tin về các số si = f’(xi), ta sử dụng trung bình trọng số CÁC ĐA THỨC GHÉP TRƠN (5) 3. PHÉP NỘI SUY SPLINE BẬC BA Nhận xét: P’i–1(xi) = P’i(xi) = si = f’(xi), i = 1, . . . , N-1 Từ điều kiện khả vi liên tục hai lần của các đa thức ghép trơn bậc ba P”i–1(xi) = P”i (xi), i = 1 , . . . , N-1  2 c3,i–1 + 4 c4,i–1hi–1 = 2c3,i – 2 c4,ihi , i = 1, . . . , N-1 Hệ PTTT trên cần thêm hai điểm x0, xN để cung cấp giá trị bằng số cho các đạo hàm biên s0, sN của các đa thức ghép trơn. Các điểm này đƣợc chọn để thỏa mãn các điều kiện “đầu mút tự do” P”0(x0) = P”N–1(xN) = 0 1N,,1i)h]x,f[xh]x,3(f[x shs)hh2(sh 1i1iiii1i 1i1iii1i1ii      PHƢƠNG PHÁP SỐ - Bài 4 30 CÁC ĐA THỨC GHÉP TRƠN (6) 3. PHÉP NỘI SUY SPLINE BẬC BA Hệ có tính chéo trội nên luôn có nghiệm duy nhất.                            )h( h2h000 h)h( h2000 00)h( h2h0 00h)h(h2h 000h)h(h2 1N2N-1N- 3N-2N--3N 323 1212 010       Sử dụng nghiệm và đ/k ta tính đƣợc các hệ số của các đa thức cục bộ Pi(x) trong công thức c4,j từ đó tính đƣợc c3,j ][ * 1-N * 2 * 1 s...,,s,s 0ss * N * 0  11hh 1i   N,,i)]x,f[x]x,3(f[x w i1iii1ii  PHƢƠNG PHÁP SỐ - Bài 4 31 Vế phải của hệ là vectơ n-1 chiều w = [w1, w2, , wN-1] T với CÁC ĐA THỨC GHÉP TRƠN (7) 3. PHÉP NỘI SUY SPLINE BẬC BA Ví dụ: Cho trƣớc các giá trị của hàm f(x), tính gần đúng f(6,46) bằng đa thức ghép trơn Spline bậc ba. 0.0031c0.0228c 0.0132s0.0651s 0.0030s0.5630s 2625,0 525,0 675,0 375,3 10200 1830 0261 0016 43 43 21                      Ta có hệ 3 đ/c Do x = 6,46 є [5; 8] = [x3; x4] nên f(6.46) ≈ P3(6.46), c1,3= f3, c2,3= s3 f(6.46) ≈1.2 - 0.0651 (1.46) + 0.0228 (1.46)2 - 0.0032 (1.46)3=1.1438 PHƢƠNG PHÁP SỐ - Bài 4 32 i xi hi fi fi+1 – fi f[xi-1, xi] wi 0 1 1 2 -0.5 -0.5 1 2 2 1.5 -0.25 -0.125 -3.375 2 4 1 1.25 -0.05 -0.05 -0.675 3 5 3 1.2 -0.075 -0.025 -0.525 4 8 2 1.125 -0.025 -0.0125 -0.2625 5 10 1.1 CÁC ĐA THỨC GHÉP TRƠN (7) 3. PHÉP NỘI SUY SPLINE BẬC BA PHƢƠNG PHÁP SỐ - Bài 4 33 0 2 4 6 8 10 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 Noi suy splinre bac 3 cua g(x) = 1+1/x CÁC ĐA THỨC GHÉP TRƠN (8) 3. PHÉP NỘI SUY SPLINE BẬC BA: SAI SỐ Có thể chỉ ra sai số trong nội suy spline bậc ba với a ≤ x ≤ b: 384 hmax5 |)(|fmax(x)|g|f(x) 4 i )4( [a,b]ξ 3 j         :N , . . . 2, i , 60 h |)(f|max|(x)g(x)f| 4 )5( b][a,ξ , 3 '      24 hmax5 |)(f|max|(x)g(x)f| 3 ii(4) b][a,ξ ' 3 '   Biên của sai số này chỉ lớn bằng 5 lần biên của sai số nội suy Hermite, nhƣng nội suy Hermite bậc ba sử dụng các thông tin về hàm f(x) nhiều gấp hai lần,do thêm các giá trị f’(xi),i = 2, , N Các đạo hàm phải xấp xỉ tốt các f’(xi). Có thể chỉ ra rằng)(xg i , 3 PHƢƠNG PHÁP SỐ - Bài 4 34 trong trƣờng hợp dãy điểm cách đều, xi = x0 + ih, với mọi i ta có double spline (vector x, vector f, double x0) { unsigned i, j, n = x.size(); double c3, c4; vector h(n - 1, 0), y(n - 1, 0), u(1, 0), l(1, 0), d(n - 2, 0), w(n - 2, 0), z(n - 2, 0), s(n, 0); // Khởi tạo các vector for (i = 0; i < n - 1; i++) h[i] = x[i+1] - x[i]; // số gia của x for (i = 0; i < n - 1; i++) y[i] = (f[i+1] - f[i]) / h[i]; // tỉ sai phân cấp 1 for (i = 0; i < n - 2; i++) d[i] = 2 * (h[i+1] + h[i]); // ĐC chính l.insert (l.end (), h.begin () + 2, h.end ()); // ĐC dưới u.insert (u.begin (), h.begin () , h.end () – 2); // ĐC trên for (i = 0; i < n - 2; i++) w[i] = 3 * (y[i] * h[i+1] + y[i+1] * h[i]); // vế phải baDuongCheo (u, d, l, w, z); // Giải hệ ba đường chéo tìm nghiệm z for (i = 1; i < n - 1; i++) s[i] = z[i-1]; // Gán các giá tri s1, ... sN-2 i = 0; // Tim khoảng thứ i [xi, xi+1] chứa x0 while (x0 > x[i+1]) i++; c4 = (s[i+1] + s[i] - 2 * y[i]) / h[i] / h[i]; c3 = (y[i] - s[i]) / h[i] - c4 * h[i]; return f[i] + s[i] * (x0 - x[i]) + c3 * (x0 - x[i]) * (x0 - x[i]) + c4 * (x0 - x[i]) * (x0 - x[i]) * (x0 - x[i]); } PHƢƠNG PHÁP SỐ - Bài 4 35 CÁC ĐA THỨC GHÉP TRƠN (9) ĐƢỜNG CONG LÀM KHỚP DỮ LIỆU (1) 1. BÀI TOÁN: Hãy phục hồi quan hệ hàm y = f(x) từ các dữ liệu đo đƣợc yi, tại các xi, i = 1, , n x x1 x2 xn y y1 y2 yn PHƢƠNG PHÁP SỐ - Bài 4 36 Các dữ liệu yi có chứa một thành phần biến đổi chậm là tính xu thế của f(x), hay các thông tin về f(x), và một thành phần biến đổi khá nhanh có biên độ tƣơng đối nhỏ đƣợc gọi là sai số hay nhiễu yi = f(xi) + ei PHƢƠNG PHÁP SỐ - Bài 4 ) x y ĐƢỜNG CONG LÀM KHỚP DỮ LIỆU (2) x x x x x x x x Biến độc lập Biến phụ thuộc y = f1(x) hay f2(x) Mi = (xi, yi) Làm sao biết đƣợc hàm nào khớp nhất ? 1. BÀI TOÁN: ĐƢỜNG CONG LÀM KHỚP DỮ LIỆU (3) 2. PHƢƠNG PHÁP BÌNH PHƢƠNG BÉ NHẤT: Tìm hàm f(x) = f(c1, c2, , ck; x) thể hiện hầu hết các thông tin về f(x) có trong dữ liệu và một phần nhỏ nhiễu. Hàm này phụ thuộc vào các tham số c1, . . . , ck một cách tuyến tính, f(c1, c2, , ck; x) = c1Φ1(x) + c2Φ2(x) + + ckΦk(x) trong đó { Φi } là một tập các hàm đƣợc lựa chọn trƣớc và { ci } là các tham số cần phải xác định. Chọn các tham số phù hợp nhất với các điểm thực nghiệm: 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ố làm cực tiểu hàm E của k biến c1, c2, , ck    n 1i 2 ik21ik21 )] x;c,...,c,(cf- y[)c,...,c,E(c PHƢƠNG PHÁP SỐ - Bài 4 38 PHƢƠNG PHÁP SỐ - Bài 4 x y ĐƢỜNG CONG LÀM KHỚP DỮ LIỆU (4) x x x x x x x x ei = yi - f(c1, .., ck; xi) Quan sát Dự đoán dựa vào x và các tham số c1, . . . , ck Sai số        n i n i ikiik x;c...,,cfyec...,,c 1 1 2 1 2 1 E tiÓu cùc lµm c¸c Tim Bình phƣơng sai số 2. PHƢƠNG PHÁP BÌNH PHƢƠNG BÉ NHẤT: ĐƢỜNG CONG LÀM KHỚP DỮ LIỆU (5) 3. LẬP CÔNG THỨC • Điểm làm cực tiểu hàm E(c1,c2,,ck) phải là điểm dừng, nghĩa là phải thoả mãn hệ k phƣơng trình • Thay x = xi, f(c1,c2,,ck; xi) = c1Ф1(xi) + c2 Ф2(xi) ++ ck Фk(xi) vào hệ trên ta đƣợc hệ PTTT đối xứng k1, j 0, c E n 1ij         j ik21 ik21i c );x,...,c,cf(c )];x,...c,c-f(c[y2           y),(Φ)Φ,(Φ...)Φ,(Φ)Φ,(Φ .................................................................... y),(Φ)Φ,(Φ...)Φ,(Φ)Φ,(Φ y),(Φ)Φ,(Φ...)Φ,(Φ)Φ,(Φ kkk2k1k 2k22212 1k12111 k21 k21 k21 ccc ccc ccc PHƢƠNG PHÁP SỐ - Bài 4 40 ĐƢỜNG CONG LÀM KHỚP DỮ LIỆU (6) 3. LẬP CÔNG THỨC Φi=[Φi(x1), , Φ i(xn)] T, i = 1, , k; y = [f(x1) , f(xn)] T (Φi, Φj) là tích vô hƣớng giữa 2 véc tơ Φi và Φj, cụ thể )(xΦ)(xΦ( kjk n 1k i  ),ΦΦ ji ))f(x(xΦ kk n 1k i  ) ,y (Φ, i PHƢƠNG PHÁP SỐ - Bài 4 41 Một số hàm làm khớp dạng tham số tuyến tính thƣờng gặp f(x) = a+ bx; f(x) = a + bx + cx2; f(x) = a + b/x f(x) = a + bcosx + csinx; hoặc có thể biến đổi về dạng tham số tuyến tính nhƣ f(x) = aebx; (a>0), f(x) = alnx + b f(x) = axb; (a>0, b>0) ĐƢỜNG CONG LÀM KHỚP DỮ LIỆU (7) 4. VÍ DỤ 1: Giả sử ta có 11 điểm thực nghiệm có thể xấp xỉ bởi đƣờng thẳng dạng f(x) = c1 + c2 x. Hãy xác định c1 và c2 bằng phƣơng pháp bình phƣơng bé nhất PHƢƠNG PHÁP SỐ - Bài 4 42 (c1, c2) ĐƢỜNG CONG LÀM KHỚP DỮ LIỆU (8) 4. VÍ DỤ 1: PHƢƠNG PHÁP SỐ - Bài 4 43       n i n i iii )xcc(yec,c 1 1 2 21 2 21 E tiÓu cùc lµm c¸c Tim 1. Xác định Φ1(xi) =1, Φ2(xi) = xi 2. Tính các tích vô hƣớng (Φ 1 , Φ 1 ), (Φ 1 , Φ 2 ),(Φ 2 , Φ 2 ), (Φ 1 , y), (Φ 2 , y) 3. Giải hệ PTTT c1(Φ1,Φ1) + c2 (Φ1,Φ2) = (Φ1,y) c1(Φ1,Φ2) + c2 (Φ2,Φ2) = (Φ2,y) ĐƢỜNG CONG LÀM KHỚP DỮ LIỆU (9) Giải ra ta được c1 = –0.7314 và c2 = 0.7437         328.05 41.04 50666 6611 i x = Φ 2 y Φ 1 (Φ 1 , Φ 1 ) (Φ 1 , Φ 2 ) (Φ 1 ,y) (Φ 2 , Φ 2 ) (Φ 2 ,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.6 11 11 7.17 1 1 11 7.17 121 78.87 Tổng 11 66 41.04 506 328.05 PHƢƠNG PHÁP SỐ - Bài 4 44 ĐƢỜNG CONG LÀM KHỚP DỮ LIỆU (10) 5. SỬ DỤNG HÀM TRENDLINE TRONG VẼ ĐỒ THỊ EXCEL PHƢƠNG PHÁP SỐ - Bài 4 45 ĐƢỜNG CONG LÀM KHỚP DỮ LIỆU (11) 5. SỬ DỤNG HÀM TRENDLINE TRONG VẼ ĐỒ THỊ EXCEL PHƢƠNG PHÁP SỐ - Bài 4 46 y = 0.743x - 0.731 0.00 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 0 5 10 15 y x y Linear (y) y = -0.001x2 + 0.091x + 0.739 0 0.5 1 1.5 2 2.5 3 0 20 40 60 80 y = f(c1, c2, , ck; x) ĐƢỜNG CONG LÀM KHỚP DỮ LIỆU (12) 5. SỬ DỤNG HÀM TRENDLINE TRONG VẼ ĐỒ THỊ EXCEL VÍ DỤ 2: Tìm đƣờng cong làm khớp dữ liệu tại 5 điểm thực nghiệm sau PHƢƠNG PHÁP SỐ - Bài 4 47 x y 5 1.201 15 1.824 30 2.624 50 2.787 70 2.210 BÀI TẬP 1. Bài tập tính toán: a) 2.2-1, 2.2-3, 2.2-5, 2.3-2, 2.6-2, 2.6-4, 6.2-2 b) Cho bảng dữ liệu quan trắc đƣợc của hàm y = f(x) x 1 2 3 5 6 y 4.75 4 5.25 19.75 36 Hãy tính f(4) bằng cách sử dụng - Đa thức nội suy Newton, - Đa thức ghép trơn spline bậc 3 2. Bài tập lập trình: Hãy viết các chƣơng trình - Nội suy Newton với bƣớc không cách đều và cách đều, - Nội suy spline bậc ba PHƢƠNG PHÁP SỐ - Bài 4 48