Giáo trình tính toán khoa học - Chương 5

Chia sẻ: Nguyen Thanhvan | Ngày: | Loại File: PDF | Số trang:24

0
236
lượt xem
82
download

Giáo trình tính toán khoa học - Chương 5

Mô tả tài liệu
  Download Vui lòng tải xuống để xem tài liệu đầy đủ

Nội suy là cơ sở của nhiều khái niệm trong giải tích số. Đó là công cụ để khôi phục các đặc trưng liên tục của một hàm số y=f(x) từ các tập hợp dữ liệu rời rạc do đo đạc hay quan sát được. Khi f(x) là một hàm phức tạp, khó tính toán và khảo sát thì cũng cần được xấp xỉ bởi một đa thức. Nội suy đơn giản nhất là nội suy bằng đa thức.

Chủ đề:
Lưu

Nội dung Text: Giáo trình tính toán khoa học - Chương 5

  1. Chương 5 NỘI SUY VÀ XẤP XỈ HÀM 5.1 NỘI SUY BẰNG ĐA THỨC Nội suy là cơ sở của nhiều khái niệm trong giải tích số. Đó là công cụ để khôi phục các đặc trưng liên tục của một hàm số y=f(x) từ các tập hợp dữ liệu rời rạc do đo đạc hay quan sát được. Khi f(x) là một hàm ph ức tạp, khó tính toán và khảo sát thì cũng cần được xấp xỉ bởi một đa thức. Nội suy đơn giản nhất là nội suy b ằng đa thức. Lý do đa thức là một h àm đơn giản: dễ tính đạo đạo h àm và nguyên hàm… Nội suy bằng đa thức là tìm một đa thức P(xi) b ậc n-1 qua n mốc nội suy xi với i  1, n thỏa mãn P(xi )= f(xi ). Nói cách khác, có thể mô tả tập các điểm dữ liệu rời rạc của h àm y = f(x) dưới dạng bảng: x x1 x2 ... xn y y1 y2 ... yn Sau đó tính các h ệ số của đa thức P(x) bậc n -1 thỏa mãn: yi=P (xi), i  1, n (5.1) Bây giời ta cần xây dựng công thức tính các hệ số của đa thức P(x). Giả sử đa thức P(x) được viết dưới dạng tường minh: P(x) = p1xn-1 + p2xn-2+...+pn-1x+ pn Từ điều kiện (5.1), để tìm các h ệ số pi của đa thức nội suy P(x) ta có thể giải hệ phương trình sau đây:  p x n 1  p2 x1  2  ... n pn 1 x1  pn  y1  11  p1 x2 1  p2 x2  2  ... n n pn 1x1  pn  y2   ... ... ... ...  n 1 p2 xn  2  ... n  p1 xn  pn 1xn  pn  yn 111
  2.  x n 1... x 2 1   p1   y1  x1 1 1      x2 1... x2 n 2 x2 1   p2   y2  hay .    ... ...   ...   ...   ... ...      n 1 xn 1   pn   yn  2  xn ... xn Ma trận hệ số của hệ phương trình chính là ma trận Vandermonde của vector x=( x1,x2,...xn). Do đó, để giải hệ phương trình trên có thể sử dụng một câu lệnh đ ơn giản trong Matlab: >> p =vander(x)\y (5.2) Đa thức nội suy đư ợc tính theo công thức trên đơn giản, nhưng khá cồng kềnh. Tính hệ số của đa thức nội suy dạng tường minh bằng giải hệ phương trình như trên sẽ mất nhiều công sức tính toán khi số nút nội suy lớn. Do đó ta cần phải nghiên cứu một số phương pháp tìm đa thức nội suy khác, đơn giản hơn. Định lý 5.1 (Tính duy nh ất của đa thức nội suy). Đa thức nội suy bậc n-1 thoả mãn (5.1) là duy nhất. Chứng minh. Thật vậy. Giả sử có 2 đa thức nội suy bậc n-1 là P(x) và Q(x) cùng thoả mãn điều kiện (5.1). Nghĩa là yi=P(xi)=Q(x) với i  1, n . Xét đa thức R(x)=P(x)-Q(x). Rõ ràng là: R(xi)=P(xi)-Q(xi) =0 , i  1, n . Như vậy R(x) là đa thức bậc không quá n-1 nhưng có tới n n ghiệm khác nhau x1,x2,...xn. Do đó R(x) =0 với mọi x, hay P(x)  Q(x). Điều đó chứng tỏ rằng đa thức nội suy bậc n-1 thỏa m ãn (5.1) là duy nh ất (đ.p.c.m). Khi số nút nội suy lớn, th ì việc giải hệ phương trình như trên tốn rất nhiều công sức. Sau đây chúng ta sẽ nghiên cứu một số ph ương pháp khác để tìm đa thức nội suy mà không cần giải hệ phương trình. 5.1.1 Đa thức nội suy Lagrange Trước hết ta xây dựng các đa thức cơ bản như sau:  x  x1   x  x2  ...  x  xi 1   x  xi 1  ...  x  xN  với Li  x   i  1, n .  xi  x1   xi  x2  ...  xi  xi 1   xi  xi 1  ...  xi  xN  112
  3. Dễ thấy các đa thức cơ b ản có tính chất:  0 khi i  j  Li x j   1 khi i  j n P (x) =  yi Li ( x) Do đó n ếu đặt (5.3) i 1 thì P(x) là đa th ức bậc không quá n thoả m ãn P(xi)=yi, với i  1, n . Do đó P(x) chính là đa thức nội suy bậc n -1 của hàm số đ ã cho. Đa thức dạng (5.3) còn gọi là đa thức nội suy Lagrange. Nó có d ạng tổng của n đa thức bậc n -1. Thí dụ 1. Xây dựng đa thức nội suy Lagrange bậc 2 từ bảng dữ liệu có dạng sau đây: x x1 x2 x3 y y1 y2 y3 Giải: Ta có các đa thức nội suy cơ bản:  x  x2   x  x3  , L ( x)   x  x1   x  x3  L1 ( x)   x1  x2   x1  x3  2  x2  x1   x2  x3   x  x1   x  x2  . và L3 ( x)   x3  x1   x3  x2   x  x2   x  x3  y   x  x1   x  x3  y   x  x1   x  x2  y . Do đó P( x)   x1  x2   x1  x3  1  x2  x1   x2  x3  2  x3  x1   x3  x2  3 5.1.2 Đa thức nội suy Newton Nội suy bằng đa thức Lagrange là một ph ương pháp khá đơn giản, sử dụng rất ít các kiến thức về đại số, nên dễ nhớ. Tuy nhiên đây lại là một phương pháp kém hiệu quả. Bây giờ xét trường hợp giữ bảng dữ liệu cũ và bổ sung thêm một nút nội suy mới (đ ể h àm số được nội suy chính xác h ơn) th ì tất cả các đa thức nội suy cơ bản lại phải tính toán lại từ đầu. Thay cho công thức nội suy dạng Lagrange ta viết đa thức nội suy P(x) dưới dạng: P(x)= a 1+ a2(x-x1) + a3(x-x1)(x-x2)+…+ aN(x-x1)(x-x2)(x-x3)...(x-xn-1) (5.4) Các hệ số a i của đa thức có thể được tính trong bảng tỉ hiệu (tỉ sai phân) theo công thức qui nạp như sau: 113
  4. y y f1  xi , xi 1   i 1 i : Tỉ hiệu cấp 1 tại xi ; xi 1  xi f [ x , x ]  f1[ xi , xi 1] f 2  xi , xi 1, xi  2   1 i 1 i  2 : Tỉ hiệu cấp 2 tại xi ;… xi  2  xi f [ x , x ,..., xi  k ]  fk 1[ xi , xi 1,..., xi  k 1] fk  xi ,..., xi  k   k 1 i 1 i  2 :Tỉ hiệu cấp k tại xi; xi  k  xi fn  2 [ x2 , x3 ,..., xn ]  f n  2[ x1, x2 ,..., xn 1] … f n 1  x1, x2 ,..., xn   :Tỉ hiệu cấp n -1 tại x1. xn  x1 Tại nút xi ch ỉ phải tính các tỉ hiệu cấp 1 đến cấp n-i. Ta lập một bảng để thu ận tiện khi tính toán các tỉ hiệu. Bảng 5-1 Bảng tỉ hiệu với 6 núy nội suy x y f1 f2 f3 f4 f5 x1 y1 f1[x1,x2] x2 y2 f2[x1,x2,x 3] f1[x2,x3] f3[x 1,x 2 ,x3,x4] x3 y3 f2[x2,x3,x4 ] f4[x1,x2,x 3,x 4,x5] f1[x3,x4] f3[x2,x3,x4,x5 ] f5[x 1,x2,x3,x4,x5,x6] x4 y4 f2[x3,x4,x5 ] f4[x2,x3,x4,x5,x6] f1[x4,x5] f3[x3,x4,x5,x6 ] x5 y5 f2[x4,x5,x6 ] f1[x5,x6] x6 y6 Khi đó các hệ số của đa thức nội suy (5.4) được xác định như sau: a 1=y1 , a2= f1[x1,x2], a3= f2[x1,x2,x3] ... an= fn-1 [x1,x2,x3,...,xn]. Công th ức (5.4) với cách tính các hệ số ai như trên gọi là công thức nội suy Newton tiến xuất phát từ x1. Với công thức này, khi thêm một nút nội suy mới xn+1 thì ta chỉ cần tính thêm một hệ số mới an+1. Khi đó b ảng tỉ hiệu chỉ phải thêm một dòng. Mặt khác, nếu chú ý ta thấy công thức Newton không đòi hỏi sự sắp xếp thứ tự về giá trị của dữ liệu, nên khi đ ảo ngược thứ tự của dữ liệu thì dạng mới của đa thức nội suy là: P(x)=b1+ b2(x-xn) + b 3(x-xn)(x-xn-1)+…+bn(x-xn)(x-xn-1)(x-xn-2)...(x-x2) ( 5.5) Khi đó các hệ số của đa thức nội suy dạng (5.5) được xác định như sau: b1=yn , b 2= f1[xn-1,xn], b3= f2[xn-2,xn-1 ,xn] ... bn= fn-1[x1,x2,x3,...,xn]. 114
  5. Công thức với cách tính các hệ số bi như trên gọi là công thức nội suy Newton lùi xu ất phát từ xn. Do đ ịnh lý 5.1 có thể thấy công thức Lagrange và các công thức Newton tiến hay lùi đ ều xác định cùng một đa thức, chỉ có h ình thức thể hiện là khác nhau và thu ận tiện áp dụng cho các trư ờng hợp khác nhau. Thí dụ 2. Cho h àm y =f(x) d ưới dạng bảng số sau: x 1 2 3 5 6 8 y 5,230 2,092 1,406 -1,202 -1,321 0,015 Hãy lập đa thức nội suy cho hàm F(x) đã cho dưới dạng: a. Tường minh; b. Lagrange; c. Newton tiến xuất phát từ x1 =1 . d. Newton lùi xuất phát từ x6 =8. Giải. a. Để tìm h ệ số tường minh của đa thức P(x) ta sử dụng công thức (5.2): >> x=[1 2 3 5 6 8]; >> y=[ 5.230 ; 2.092; 1.406; -1.202 ; -1.321; 0.015]; >> p=vander(x)\y p= -0.0187 0.4201 -3.4803 13.2919 -24.3719 19.3890 Do đó : P  x   0, 0187 x5  0, 4201x 4  3, 4803 x3 + 13.2919 x 2  24,3719 x  19,3890 b . Lập đa thức nội suy Lagrange: ( x  2)( x  3)( x  5)( x  6)( x  8) ( x  1)( x  3)( x  5)( x  6)( x  8) P( x)  5, 230 2, 092  (1  2)(1  3)(1  5)(1  6)(1  8) (2  1)(2  3)(2  5)(2  6)(2  8) 115
  6. ( x  1)( x  2)( x  5)( x  6)( x  8) ( x  1)( x  2)( x  3)( x  6)( x  8) +1, 406 1, 202 (3  1)(3  2)(3  5)(3  6)(3  8) (5  1)(5  2)(5  3)(5  6)(5  8) ( x  1)( x  2)( x  3)( x  5)( x  8) ( x  1)( x  2)( x  3)( x  5)( x  6) 1,321  0, 015 (6  1)(6  2)(6  3)(6  5)(6  8) (8  1)(8  2)(8  3)(8  5)(8  6) hay P( x)  0,0187( x  2)( x  3)( x  5)( x  6)( x  8)  0,0291( x  1)( x  3)( x  5)( x  6)( x  8) 0,0234( x  1)( x  2)( x  5)( x  6)( x  8)  0,0167 ( x  1)( x  2)( x  3)( x  6)( x  8) 0, 0110( x  1)( x  2)( x  3)( x  5)( x  8)  0, 00001( x  1)( x  2)( x  3)( x  5)( x  6) c. Lập đa thức nội suy Newton tiến xuất phát từ x1 =1. x y f1 f2 f3 f4 f5 1 5,230 -3,1380 1,2260 -0,3580 2 2 ,092 -0,6860 0,1016 -0,2060 0,1502 -0,0187 3 1 ,406 -1,3040 -0,0295 0,3950 -0,0265 5 -1,202 -0,1190 0,2623 6 -1,321 0,6680 8 0 ,015 Từ bảng tỉ hiệu ta có: P  x   5, 230  3,1380  x  1  1, 2260  x  1  x  2  0,3580  x  1  x  2   x  3 0,1016  x  1  x  2   x  3  x  5   0, 0187  x  1  x  2   x  3  x  5   x  6  d. Lập đa thức nội suy Newton lùi xu ất phát từ x6 =8. Từ bản g tỉ hiệu ta có: P  x   0, 015  0, 0668  x  8   0, 2623  x  8   x  6   0, 0265  x  8   x  6   x  5  0, 0295  x  8   x  6   x  5   x  3  0, 0187  x  8   x  6   x  5   x  3  x  2  5.1.3 Sai số nội suy Giả sử P(x) là đa th ức nội suy của h àm f(x) tại n nút nội suy x1,x2, ...,xn; xi[ a,b] và hàm f(x) kh ả vi đến cấp n . Khi đó có thể chứng minh đ ược rằng: ( x ) R(x)= f(x) - P(x) = f   ( c ) n , với c [ a,b]; n! 116
  7. trong đó (x) =(x-x1)(x-x2)...(x-xn-1)(x-xn) là đa th ức bậc n và có n nghiệm tại các nút nội suy x1, x2,...,xn. Do f(n)(c) là hằng số nên dáng điệu của sai số của nội suy R(x) phụ thuộc vào dáng điệu của h àm (x). Xét trường hợp các nút nội suy cách đều. Khi đó h àm (x) có biên độ nhỏ dần ở giữa khoảng nội suy và lớn dần khi đi ra hai biên (xem hình 5.1). Từ đó có một vấn đề nảy sinh là: n ếu được phép chọn các nút nội suy thì nên chọn như th ế nào đ ể sai số nội suy trở thành bé nhất. Điều đó dẫn đến việc giải b ài toán: max|(x )|  min Kết quả gải b ài toán như sau: Nếu a=-1 vàc b=1, thì các nút nội suy “ tối ưu” được chọn như sau: -  2i  1    , i  1, n . (5.6) xi  cos   n 2 Đó chính là các nghiệm của đa thức Chebysev bậc n: 1 Tn  x   cos  n.arccos x  2n 1 Hình 5.1 Đồ thị hàm y= (x) trên lưới đều Điều đó nghĩa là phân bố các nút nội suy “tối ưu” là thưa ở giữa, dày d ần khi tiến sang 2 biên của khoảng nội suy. Khi đó ta có đánh giá sai số : 1   x   Tn  x   n 1 2 117
  8. 2x  a  b - Nếu a -1 hoặc b1 thì tiến hành đổi biến t  đ ể đưa x [a,b] ba  2i  1   về t [-1,1] rồi chọn ti theo công thức (5.4) ti  cos   , i  1, n . Sau đó tiến  n 2  b  a  ti   a  b  . Khi đó sai số nội suy được tính hành đổi biến ngược lại xi  2 theo công thức sau: n b  a  R( x)  Pn  x   f  x   f   ( x) n (5.7) n !2 2 n 1 5.1.4 Một số hàm số tính toán với đa thức Trong Matlab đã có sẵn các hàm nội trú thuận tiện cho thực hiện tính toán các đa thức nội suy nói riêng và tổ hợp các đa thức nói chung.  Hàm POLYFIT Cú pháp: p = polyfit(x,y,N) Giải thích: Hàm POLYFIT tính hệ số của đa thức xấp xỉ hàm cho bởi cặp 2 vector cùng cấp x và y. - Nếu N  length(x) –1, thì hàm tính vector p là vector h ệ số của đa thức nội suy bậc N: P(x)= p 1xN+p2xN-1+...+pNx+pN+1 thỏa mãn P(xi)= yi ,i=1,2,..., length(x). Nếu N < length(x)-1 , hàm sẽ tính vector các hệ số p của đa thức xấp xỉ - tốt nhất đối với dữ liệu theo nghĩa bình phương tối thiểu.  Hàm POLYVAL Cú pháp: y = polyval (p,x) Giải thích. Hàm POLYVAL tính giá trị của đa thức có hệ số cho bởi vector p tại tất cả các giá trị của vector x, nghĩa là : y = p1xn-1 + p2xn-2+...+pn-1x+ pn  Hàm CONV Cú pháp: c = conv (a,b) 118
  9. Giải thích. Hàm CONV thực hiện việc nhân hai đa thức có hệ số a và b thành đa thức hệ số c. Thí dụ 3. Nhân hai đa th ức (3x2+4x+5).(2x3-3x2+2) >> a = [0 3 4 5]; >> b = [2 -3 0 2]; >> c = conv(a,b) c= 0 6 -1 -2 -9 8 10  Hàm DECONV Cú pháp: [q,r] =deconv (a,b) Giải thích. Hàm DECONV chia đa thức có hệ số là vector a cho cho thức có hệ số là vector b, được đa thức có hệ số là vector q và phần dư là đa th ức có hệ số là vector r. Thí dụ 4. Chia đa thức (2x2+3x+6)/(2x+3) >> a = [2 3 6]; >> b = [2 3]; >> [q,r] = deconv (a,b) q= 1 0 r= 0 0 6  Hàm ROOTS Cú pháp: x = roots(p) Giải thích. Hàm ROOTS tính tất cả các nghiệm thực và phức của đa thức có hệ số là vector p trả về kết quả là vector nghiệm x. Thí dụ 5. >> roots([3 4 1]) 119
  10. ans = -1.0000 -0.3333 >> x= roots([1 2 3 2 1]) x= -0.5000 + 0.8660i -0.5000 - 0.8660i -0.5000 + 0.8660i -0.5000 - 0.8660i 5.1.4 Biểu diễn hàm số trên đồ thị  Vẽ đồ thị phẳng bằng lệnh plot plot(x,y,'symbol'): Vẽ đồ thị hàm y đối với x. 2 vector thực x và y cùng cỡ. 'Symbol' là xâu qui định kiểu màu hoặc đường vẽ: y yellow w white * star m magenta g green o circle c cyan b blue x x-mark r red k b lack plot(y,'symbol') : Nếu y là vector thực lệnh vẽ đồ thị của hàm y đối với x là số thứ tự của toạ độ của y. Nếu y là vector phức thì tương đương với lệnh: plot(real(y),imag(y), 'symbol')  Xác định tỉ lệ trên đồ thị bằng lệnh axis axis([xmin xmax ymin ymax]) : Thay đổi lại tỷ lệ của các trục toạ độ. axis auto (mặc định): Matlab tự chọn một tỷ lệ thích hợp nhất cho đồ thị. : Trả về một vectơ hàng mô tả tỷ lệ của đồ thị hiện tại, có dạng V = a xis [xmin xmax ymin ymax]. : Các trục toạ độ có cùng đơn vị. axis equal : Hộp đồ thị hình vuông. axis square Khi xem đồ thị có thể dùng các công cụ “+” và “ –“ để phóng to thu nhỏ vùng đồ thị muốn xem. 120
  11.  Các lệnh liên quan khác liên quan đến đồ thị hold on / off : Xếp chồng đồ thị, mặc định là hold o ff : thay đồ thị cũ. grid on / off : Hiện (on , mặc định ) hay ẩn (off) lư ới. xlabel('text'), ylabel('text') , zlabel('text') : Gắn nh ãn cho các trục toạ độ. : Gắn tiêu đ ề cho đồ thị. tiltle('text')  Minh hoạ đa thức nội suy. Như đa trình bày ở trên, việc nội suy bằng đa thức với các nút nội suy cách đều có thể gây ra sự d ao động không mong muốn cho đa thức nội suy của một hàm số. Thí dụ sau đây sẽ là lời cảnh báo về điều này. Thí dụ 6. Hãy n ội suy hàm số 1 f x  với x [-1,1]. 1  16 x 2 Giải. Ta tạo một tập dữ liệu tại 20 mốc cách đều trên đoạn [-1,1] là: 2(i  1) 1 , i  1, 20 . xi   1, fi  1  16 xi2 19 Từ tập dữ liệu này n ội suy bằng đa thức bậc 19. Kết quả minh hoạ trên đồ thị h ình thị 5.2 cho thấy: m ặc dù hàm P19(x) thoả mãn P 19(xi) = yi , n hưng giá trị xấp xỉ bởi nội suy tại gần 2 đầu mút của khoảng nội suy không tốt. Sau đó, ta tạo lại tập dữ liệu như sau: (i  1) 1 xi  cos , fi  , i  1, 20 1  16 xi2 19 và vẽ lại đồ thị của đa thức nội suy mới. Chương trình tạo file dữ liệu ban đầu của hàm số: % Matlab Code to Produce Data file n=20; x=zeros(n,1); y=x; for t=1:n; x(t)=2*(t-1)/19-1; y(t)=1/(1+16*x(t)^2); end; %% Hoặc save('dulieu.mat', 'x', 'y' save dulieu x y; Chương trình vẽ đồ thị của các hàm nội suy: % Demonstrate the "failure" of Polynomial Interpolation on equidistant grids 121
  12. clear; load dulieu; n=length(x); axis([-1 1 -0.5 4]); plot(x,y,'o'); %% Đánh dấu các điểm dữ liệu bằng hình chữ “ o” hold on; grid on; coef=polyfit(x,y,n-1); xx=[-1:0.01:1]; yy=polyval(coef,xx); plot(xx,yy,'r'); xlabel(' Truc X'); ylabel(' Truc Y'); title( ' Noi suy da thuc '); for t=1:n x(t)=cos((t-1)*pi/19); y(t)=1/(1+16*x(t)^2); end coef1=polyfit(x,y,n-1); yy=polyval(coef1,xx); plot(xx,yy,'g'); hold off; H ình 5.2 “Th ất bại” của nội suy trên lưới đều và ưu điểm của các nút nội suy “tối ưu” Kết quả chạy chương trình được minh họa trong hình 5.2 cho thấy với các nút nội suy mới đồ thị đường nội suy đẹp h ơn rất nhiều. 122
  13. 5.2 NỘI SUY BẰNG HÀM SPLINE B ẬC 3 Nội suy bằng đa thức bằng đa thức có nhược điểm là khi có nhiều mốc nội suy thì bậc của đa thức rất cao, nên không thuận tiện cho việc tính toán. Phương pháp spline là thực hiện ghép nối trơn tru những đa thức bậc thấp, từng khúc để nội suy hàm đã cho. Cho hàm f(x) liên tục trên đoạn [a,b]. Xét một phân hoạch trên [a,b] ={a =x1
  14. mj m j 1  x j 1  x   x  xj  S ''j ( x)  (5.10) hj hj Từ đó suy ra: mj m j 1 3 3  x j 1  x  x  xj      (5.11)   j x j 1  x   j x  x j S j ( x)   6h j 6h j Để tính  j và  j , lần lượt thay x=xj và x=xj+1 vào (5.11) ta sẽ được: mj m j 1 3 3  x j 1  x  x  xj  S j ( x)    6h j 6h j (5.12) m h2  m h2  1   y j  j j  x j 1  x  1  y j 1  j 1 j  x  x j      hj  6 hj  6     Và do đó: m h2  m h2    m j m j 1 1  yj  j j  1 2 2  y j 1  j 1 j   x j 1  x  x  xj  S ' j ( x)    2h j 2h j hj 6  hj 6       Vì S '( x j  0)  S '( x j  0) với j  2 ,n  1 h ay S ' j 1( x j )  S ' j ( x j ) , suy ra: y j  y j 1 y j 1  y j h j 1 h j 1 hj hj m j 1  mj    m j 1  m j  6 3 h j 1 6 3 hj h j 1  h j y j 1  y j y j  y j 1 h j 1 hj (5.13) m j 1  mj  m j 1   6 3 6 hj h j 1 h j 1  h j Nếu chia cả hai vế của (5.11) cho và đặt 6 hj h j 1 j  ,  j  1  j  , h j 1  h j h j 1  h j  y j 1  y j y j  y j 1  6 và (5.14) dj    , j  2, n  1   h j 1  h j 1 h j hj   thì (5.13) trở th ành: (5.15)  j m j 1  2m j   j m j 1  d j , j  2, n  1 Đây là một hệ gồm n-2 phương trình tuyến tính của n ẩn số mj. Để bổ xung thêm 2 (hay m-1) phương trình nữa, có thể sử dụng một trong hai cách chọn 2 điều kiện biên như sau: 124
  15. a) Buộc S '(a) =y'(a) và S'(b) =y'(b) khi đó :  6y y  2m1  m2   2 1  y '1   d1  h1  h1    6 yn  yn 1    mn 1  2 mn  h  y 'n  h   dn n 1  n 1   d1 Buộc S”(a)=y”(a)= b) h ay 2. m1+0.m2 =d1 , 2 dn và buộc S ”(b) =y”(b) = hay 0 .mn-1 + 2mn =dn. 2 Như vậy để tính các hệ số mj với j  1, n , thì cần giải hệ phương trình: 2   m1   d1  1       m2   d 2   2 2 2   m   d  3 2 3  3    3  (5.16)  ...   ...   ...   ... ... ... ... ...   2 n 1   mn 1   d n 1  n 1         2   mn   d n  n   Cần chú ý rằng (5.13) là h ệ phương trình có ma trận hệ số dạng ma trận thưa và chéo trội, nên có th ể hệ giải bằng phương pháp lặp đơn Jacoby hay phương pháp lặp Gauss-Seidel hội tụ. Khi các nút nội suy cách đều thì: 1 ,j  2, n  1 . j  j  2 Có th ể dễ d àng nhận thấy việc xác định các hệ số của h àm spline bậc 3 khá phức tạp và khó nhớ. Do tính ứng dụng lớn của nội suy bằng hàm spline b ậc 3, nên trong Matlab đã cài đặt sẵn một h àm spline b ậc 3.  Hàm SPLINE trong Matlab Cú pháp: yy = spline(x,y,xx) Giải thích. Hàm SPLINE tính giá trị h àm nội suy spline bậc 3 xác định từ tập dữ liệu (x,y). x và y p hải là 2 vector cùng cỡ. Hàm trả về vector yy giá trị của hàm spline bậc 3 tại mọi phần tử của vector xx. 125
  16. Thí dụ 7. Cài đặt ch ương trình nội suy h àm số từ tập dữ liệu (x,y) trong file dulieu.m bằng hàm spline bậc 3. Giải. Soạn file chương trình có nội dung: % Malab code for Cubic Spline Interpolation clear load dulieu; n = length(x); xx = [x(1): (x(n)-x(1))/200:x(n)] ; yy = spline(x,y,xx); plot(x,y,'o'); grid on; hold on; plot(xx,yy,'r'); hold off; xlabel(' Truc X'); ylabel(' Truc Y'); title( ' Noi suy bang Spline bac 3 '); Kết quả chạy ch ương trình được minh họa trong hình 5.3. Đồ thị cho ta thấy việc xấp xỉ bằng spline bậc 3 là khá tốt. H ình 5.3 Nội suy bằng spline bậc 3 126
  17. 5.3 XẤP XỈ HÀM 5.3.1 Phân tích hồi qui tuyến tính Hồi qui tuyến tính là tìm một đ ường thẳng phù hợp với một tập các dữ liệu đã thống kê. Khi nghiên cứu hồi qui tuyến tính có thể sử dụng nhiều phương pháp khác nhau, tuy nhiên trong nghiên cứu thực nghiệm ngư ời ta thường sử dụng phương pháp bình phương tối thiểu. Cho một tập gồm n cặp dữ liệu thống kê {(xi,yi), i  1,n } cần phải tìm các hệ số a và b của một đường thẳng hồi qui có dạng y(x)=ax+b sao cho đường thẳng đó “phù hợp nhất ” với tập dữ liệu đã cho. Đây là một bài toán được áp dụng nhiều trong các ngành khoa học thực nghiệm. Đặt ei = yi - (axi + b) là sai số của đường hồi qui tại điểm xi . Cần phải tìm đường hồi qui sao cho tổng b ình phương của các sai số trở nên bé nhất, nghĩa là: n 2 L(a,b) =   yi   axi  b   min i 1 Để tìm lời giải tối ưu của b ài toán cực trị 2 biến ở trên ta cần giải hệ 2 phương trình của a và b :  L n   2{ yi  (axi  b)}( xi )  0  a i 1    L n   b  2{ yi  (axi  b)}(1)  0   i 1  n n n    xi2  a   xi  b     xi yi  i 1       i 1  i 1 ha y (5.17)  n n      xi  a  n.b  yi    i 1   i 1 Thí dụ 8. Tạo một text file có tên data.m gồm 2 cột dữ liệu: 0 2.38471 1.0e-01 1.82271 2.0e-01 2.03641 3.0e-01 2.41649 4.0e-01 2.17773 5.0e-01 1.98777 127
  18. Sau đây là chương trình xây dựng đư ờng hồi qui tuyến tính: % Matlab code for linear regression clear; load data.m; x = data(:,1); y=data(:,2); p lot(x,y,'*'); a xis([0 0.5 0 4]); hold on; grid on; a11=sum(x.*x) ; a12=sum(x); n =length(x); b1=sum(x.*y) ; b2=sum(y); A=[ a11 a12; a12 n]; B=[b1;b2]; sol = A\B; a =sol(1) ; b=sol(2); x0 = x(1) ; xn=x(n); y0 = a*x0 + b; yn=a*xn+b; p lot([x0 xn],[y0 yn],'r'); hold off; xlabel(' X-Axis'); ylabel('Y-Axis'); title( ' Linear Regession Analysis '); H ình 5.4 Đường thẳng hồi qui 128
  19. 5.3.2 Xấp xỉ hàm bậc cao. Phân tích hồi qui tuyến tính là một trường hợp xấp xỉ hàm số ch ưa biết (hàm thực nghiệm) bởi một hàm bậc nhất. Bây giờ ta xét một trường hợp xấp xỉ hàm thực nghiệm bởi một đa thức bậc 2. Từ đó học viên có thể tự xây dựng công thức tính cho các xấp xỉ h àm b ậc cao hơn và một số dạng khác cũng theo phương pháp bình phương tối thiểu. Để tìm các h ệ số a, b, c của đường parabol y(x) = ax2+bx +c sao cho nó phù hợp nhất theo nghĩa bình phương tối thiểu với một tập gồm n cặp dữ liệu thực nghiệm {(xi,yi)}, i  1, n ta cần giải b ài toán tìm cực trị: n 2   yi   axi2  bxi  c  L(a,b,c) =  min i 1 Từ đó dẫn đến cần giải hệ 3 phương trình tuyến tính của 3 ẩn a, b, c:  L n 2{ yi  (axi2  bxi  c)}( xi2 )  0    a i 1   L n 2{ yi  (axi2  bxi  c)}( xi )  0  (5.19)   b i 1  n  L 2{ yi  (axi2  bxi  c)}(1)  0   c   i 1 Viết lại hệ phương trình dạng ma trận: n 4 n n n 2   xi3 xi2       xi  xi yi   i 1   i 1  i 1 i 1   a    n 3 n n    n   xi2    (5.20) xi   b    xi xi yi     c   i 1  i 1  i 1 i 1    n n  n  x2   yi    xi  n i     i 1   i 1  i 1 Bằng phương pháp tương tự bạn có thể tự xây dựng được hệ phương trình tính hệ số của đường hồi qui bậc 3,4,5... Để tiện cho việc tính toán các hệ số của hệ phương trình (5.20) người ta thường lập bảng tính toán có dạng: 129
  20. Bảng 5-2 Bảng tính toán cho đa thức hồi qui bậc 2 xi2 xi3 xi4 xi2 yi xi yi xi yi x1 y1 … … … … … x2 y2 … … … … … … … … … … … … xn yn … … … … … n n n n n n n xi2 xi3 xi4  xi2 yi  xi  yi     xi yi i 1 i 1 i 1 i 1 i 1 i 1 i 1 Thí dụ 9. Tính hệ số của đa thức bậc 2 xấp xỉ tốt nhất hàm số y=f(x) được cho dưới dạng bảng: x 1 1,5 2 2,5 3 3,5 y 1,2341 3,9242 2,4563 - 0,2224 -1,3215 0,5506 Giải. Đầu tiên ta lập bảng tính toán: xi yi xi2 xi3 xi4 xi2 yi xi yi 1 1,2341 1 ,00 1 ,000 1,0000 1 ,2341 1,2341 1,5 3,9242 2 ,25 3 ,375 5,0625 5 ,8863 8,8294 2 2,4563 4,00 8,000 16,0000 4,9126 9,8252 2,5 -0,2224 6 ,25 15,625 39,0625 -0,5560 -1,3900 3 -1,3215 9,00 27,000 81,0000 -3,9645 -11,8935 3,5 0,5506 12,25 42,875 150,0625 1 ,9271 6,7448 13,5 6,6213 34,75 97,875 292,1875 9 ,4396 13,3501 - Giải hệ phươn g trình:  292,1875 97,875 34, 75   a  13,3501      34, 75 13,5   b    9, 4396   97,875  34, 75 6   c   6, 6213  13,5      ta được các hệ số a= -0,1868, b= -0.4071 và c = 3,1013. 130

CÓ THỂ BẠN MUỐN DOWNLOAD

Đồng bộ tài khoản