Bài giảng Nội suy và xấp xỉ hàm

210 CHƯƠNG 3: NỘI SUY VÀ XẤP XỈ HÀM §1. NỘI SUY LAGRANGE  Trong thực tế nhiều khi ta cần tính giá trị của hàm y = f(x) tại một giá trị  x  trong một đoạn  [a, b] nào đó mà chỉ biết một số nhất định các giá  trị của  hàm  tại một  số  điểm  cho  trước.  Các  giá  trị  này  được  cung  cấp  qua  thực  nghiệm hay tính toán. Vì vậy nảy sinh vấn đề toán học là trên đoạn a ≤ x ≤ b  cho một loạt các điểm xi ( i = 0, 1, 2...) và tại các điểm xi này giá trị của hàm là  yi = f(xi) đã biết và ta cần tìm y = f(x) dựa trên các giá trị đã biết đó. Lúc đó ta  cần tìm đa thức :    Pn(x) = aoxn + a1xn‐1  + …+an‐1x  + an   sao cho Pn(xi) = f(xi) = yi. Đa thức Pn(x) được gọi  là đa thức nội suy của hàm  y=f(x). Ta chọn đa  thức để nội suy hàm y =  f(x) vì đa  thức  là  loại hàm đơn  giản, luôn có đạo hàm và nguyên hàm. Việc tính giá trị của nó theo thuật toán  Horner cũng đơn giản.    Bây giờ ta xây dựng đa thức nội suy kiểu Lagrange. Gọi Li là đa thức:  )xx)...(xx)(xx)...(xx( )xx)...(xx)(xx)...(xx(L ni1ii1ii0i n1i1i0 i −−−− −−−−= +− +−     Rõ ràng là Li(x) là một đa thức bậc n và :  ⎩⎨ ⎧ ≠ == ij0 ij1 )x(L ji   Ta gọi đa thức này là đa thức Lagrange cơ bản.  Bây giờ ta xét biểu thức :    ∑ = = n 0i iin )x(L)x(f)x(P   Ta  thấy Pn(x)  là một đa  thức bậc n vì các Li(x)  là các đa  thức bậc n và  thoả mãn điều kiện Pn(xi) = f(xi) = yi. Ta gọi nó là đa thức nội suy Lagrange.  Với n = 1 ta có bảng   x  x0  x1  y  y0  y1  Đa thức nội suy sẽ là :    P1(x) = yoL0(x) + y1L1(x1)  10 1 0 xx xxL − −=         01 0 1 xx xxL − −=   211 nên  01 0 1 10 1 01 xx xxy xx xxy)x(P − −+− −=   Như vậy P1(x) là một đa thức bậc nhất đối với x  Với n = 2 ta có bảng   x  x0  x1  x2  y  y0  y1  y2  Đa thức nội suy sẽ là :    P2(x) = yoL0(x) + y1L1(x1) + y2L2(x2)  )xx)(xx( )xx)(xx(L 2010 21 0 −− −−=   )xx)(xx( )xx)(xx(L 2101 20 1 −− −−=   )xx)(xx( )xx)(xx(L 1202 10 2 −− −−=   Như vậy P1(x) là một đa thức bậc hai đối với x.   Ta xây dựng hàm  lagrange() để  thực hiện việc nội suy hàm  theo  thuật  toán  Lagrange:  function [l, L] = lagrange(x, y)  %Dua vao : x = [x0 x1 ... xn], y = [y0 y1 ... yn]  %ket qua: l = He so cua da thuc Lagrange bac n  % L = Da thuc Lagrange   n = length(x) ‐ 1; %bac cua da thucl  l = 0;  for m = 1:n + 1      p = 1;      for k = 1:n + 1          if k ~= m              p = conv(p, [1 ‐x(k)])/(x(m) ‐ x(k));           end      end      L(m, :) = p; %da thuc Lagrange       l = l + y(m)*p;   end  212 Cho hàm dưới dạng bảng:  x  ‐2  ‐1  1  2  y  ‐6  0  0  6  và tìm y(2.5) ta dùng chương trình ctlagrange.m:    clear all, clc  x = [‐2 ‐1 1 2];  y = [‐6 0 0 6];  l = lagrange(x, y);  yx = polyval(l, 2.5)  §2. NỘI SUY NEWTON  Bây giờ ta xét một cách khác để xây dựng đa thức nội suy gọi là phương  pháp Newton. Trước hết ta đưa vào một khái niệm mới là tỉ hiệu     Giả sử hàm y = y(x) có giá trị cho trong bảng sau:  x  x0  x1  x2  …  xn‐1  xn  y  y0  y1  y2  …  yn‐1  yn  Tỉ hiệu cấp 1 của y tại xi, xj là :  ji ji ji xx yy ]x,x[y − −=     Tỉ hiệu cấp hai của y tại xi, xj, xk là :  ki kjji kji xx ]x,x[y]x,x[y ]x,x,x[y − −=   v.v.      Với y(x) = Pn(x) là một đa thức bậc n thì tỉ hiệu cấp 1 tại x, x0 :  0 0nn 0n xx )x(P)x(P]x,x[P − −=   là một đa thức bậc (n ‐ 1). Tỉ hiệu cấp 2 tại x, x0, x1 :  1 10n0n 10n xx ]x,x[P]x,x[P]x,x,x[P − −=   là một đa thức bậc (n‐2) v.v và tới tỉ hiệu cấp (n + 1) thì :  213     Pn[ x, xo,.., xn] =  0  Từ các định nghĩa tỉ hiệu ta suy ra :    Pn(x) = Pn(x0) + ( x‐ x0)Pn[x, xo]    Pn[x, x0] = Pn[x0, x1] + ( x ‐ x1)Pn[x, xo,x1]    Pn[x, xo, x1] = Pn[x0, x1, x2] + ( x ‐ x2)Pn[x, xo, x1, x2]    ............    Pn[x, xo,.., xn‐1] = Pn[x0, x1,.., xn] + ( x ‐ xn)Pn[x, xo,.., xn]  Do   Pn[ x, xo,.., xn] =  0 nên từ đó ta có :  Pn(x) = Pn(x0) + (x ‐ x0)Pn[xo, x1] + (x ‐ x0)(x ‐ x1)Pn[x0, x1, x2] +…  +(x ‐ x0)…(x ‐ xn‐1)Pn[x0,…, xn]  Nếu Pn(x) là đa thức nội suy của hàm y = f(x) thì:    Pn(xi) = f(xi) = yi với i = 0 ÷ n     Do đó các tỉ hiệu từ cấp 1 đến cấp n của Pn và của y  là trùng nhau và  như vậy ta có :     Pn(x) = y0 + (x ‐ x0)y[x0, x1] + (x ‐ x0)(x ‐ x1)y[x0, x1, x2] + .. +                         (x ‐ x0)(x ‐ x1)...(x ‐ xn‐1)y[x0,..,xn]  Đa thức này gọi là đa thức nội suy Newton tiến xuất phát từ nút x0 của  hàm y = f(x). Ngoài đa thức tiến còn có đa thức nội suy Newton lùi xuất phát  từ điểm xn có dạng như sau :    Pn(x) = yn + (x ‐ xn)y[xn, xn‐1] + (x ‐ xn)(x ‐ xn‐1)y[xn, xn‐1,xn‐2] +..+                       (x ‐ xn)(x ‐ xn‐1)...(x ‐ x1)y[xn,.., x0]  Trường hợp các nút cách đều thì xi = x0 +  ih với  i = 0, 1,.., n. Ta gọi sai  phân tiến cấp 1 tại i là :      ∆yi = yi+1 ‐ yi  và sai phân tiến cấp hai tại i:      ∆2yi = ∆(∆yi) = yi+2 ‐ 2yi+1 + yi      .........  và sai phân tiến cấp n là :      ∆nyi = ∆(∆n‐1yi)  Khi đó ta có:      [ ] h yx,xy 010 ∆=         [ ] 20 2 210 h2 yx,x,xy ∆=      ...........  214   [ ] n0 n n210 h!n yx,...,x,x,xy ∆=   Bây giờ đặt x = x0 + ht  trong đa thức Newton tiến ta được:    0 n 0 2 000n y!n )1nt()1t(ty !2 )1t(tyty)htx(P ∆+−⋅⋅⋅−+⋅⋅⋅+∆−+∆+=+     thì ta nhận được đa thức Newton tiến xuất phát từ x0 trong trường hợp nút  cách đều. Với n = 1 ta có :    P1(x0 + ht) = y0 + ∆y0  Với n = 2 ta có:  0 2 000n y!2 )1t(tyty)htx(P ∆−+∆+=+   Một cách tương tự ta có khái niệm các sai phân lùi tại i:      ∇yi = yi ‐ yi‐1      ∇2yi = ∇(∇yi) = yi ‐ 2yi‐1 + yi‐2      .........      ∇nyi = ∇(∇n‐1yi)  và đa thức nội suy Newton lùi khi các điểm nội suy cách đều:  n n n 2 nn0n y!n )1nt()1t(ty !2 )1t(tyty)htx(P ∇−+⋅⋅⋅++⋅⋅⋅+∇++∇+=+   Ta xây dựng hàm newton() để nội suy:  function [n,DD] = newton(x,y)  %Dua vao : x = [x0 x1 ... xN]  % y = [y0 y1 ... yN]  %Lay ra: n = he so cua da thuc Newton bac N  N = length(x) ‐ 1;  DD = zeros(N + 1, N + 1);  DD(1:N + 1, 1) = yʹ;  for k = 2:N + 1      for m = 1: N + 2 ‐ k           DD(m,k) = (DD(m + 1, k ‐ 1) ‐ DD(m, k ‐ 1))/(x(m + k ‐ 1) ‐ x(m));      end  end  a = DD(1, :);  n = a(N+1);   for k = N:‐1:1   215     n = [n a(k)] ‐ [0 n*x(k)];   end  Cho hàm dưới dạng bảng:  x  ‐2  ‐1  1  2  4  y  ‐6  0  0  6  60  Ta dùng chương trình ctnewton.m để nội suy:    clear all, clc  x = [‐2 ‐1 1 2 4];  y = [‐6 0 0 6 60];  a = newton(x, y)  yx = polyval(a, 2.5)  §3. NỘI SUY AITKEN ‐ NEVILLE  Một  dạng  khác  của  đa  thức  nội  suy  được  xác  định  bằng  thuật  toán  Aitken  ‐ Neville. Giả sử ta có n điểm đã cho của hàm f(x). Như vậy qua hai  điểm  x0 và  x1  ta  có  đa  thức nội  suy Lagrange  của hàm  f(x)  được viết dưới  dạng:  01 11 00 01 xx xxy xxy )x(P − − − =   Đây là một đa thức bậc 1:  01 0 1 10 1 001 xx xxy xx xxy)x(P − −+− −=   Khi x = x0 thì:  0 01 011 000 001 yxx xxy xxy )x(P =− − − =   Khi x = x1 thì:  1 01 111 100 101 yxx xxy xxy )x(P =− − − =   Đa thức nội suy Lagrange của f(x) qua 3 điểm x0, x1, x2 có dạng:  216 02 212 001 012 xx xx)x(P xx)x(P )x(P − − − =   và là một đa thức bậc 2:  )xx)(xx( )xx)(xx(y )xx)(xx( )xx)(xx(y )xx)(xx( )xx)(xx(y)x(P 1202 10 2 2101 20 1 2010 21 0012 −− −−+−− −−+−− −−=   Khi x = x0 thì:    0 02 0212 000 0012 yxx xx)x(P xxy )x(P =− − − =   Khi x = x1 thì:    1 02 121 101 1012 yxx xxy xxy )x(P =− − − =   Khi x = x2 thì:    2 02 222 20201 2012 yxx xxy xx)x(P )x(P =− − − =   Tổng quát đa thức nội suy Lagrange qua n điểm là:  02 nn..12 0)1n..(01 n..012 xx xx)x(P xx)x(P )x(P − − − = − Như  vậy  ta  có  thể  dùng  phép  lặp  để  xác  định  lần  lượt  các  đa  thức  Lagrange. Sơ đồ tính toán như vậy gọi là sơ đồ Neville ‐ Aitken.  Ta xây dựng hàm aitkenneville() để nội suy:  function a = aitkenneville(xData, yData, x)  % Tra ve gia tri noi suy tai x.  % Cu phap: y = aitkenneville(xData, yData, x)  n = length(xData);  y = yData;  for k = 1:n‐1      y(1:n‐k) = ((x ‐ xData(k+1:n)).*y(1:n‐k)...      + (xData(1:n‐k) ‐ x).*y(2:n‐k+1))...      ./(xData(1:n‐k) ‐ xData(k+1:n));  217 end  a = y(1);  Cho các cặp số (1, 3), (2, 5), (3, 7), (4, 9) và (5, 11), để tìm y tại x = 2.5 ta dùng  chương trình ctaitkennevile.m:  clear all, clc  x = [1  2  3  4];  y = [3  5  7  9];  yx = aitkenneville(x, y, 2.5)  §4. NỘI SUY BẰNG ĐƯỜNG CONG SPLINE BẬC BA    Khi số điểm cho trước dùng khi nội suy tăng, đa thức nội suy có dạng  sóng và sai số tăng. Ta xét hàm thực:    2 1f31(x) 1 8x = +   và nội suy nó bằng thuật toán Newton nhờ chương trình  cttestintp.m  %Noi suy Newton  x1 = [‐1 ‐0.5 0 0.5 1.0];   y1 = f31(x1);  n1 = newton(x1,y1)  x2 = [‐1 ‐0.75 ‐0.5 ‐0.25 0 0.25 0.5 0.75 1.0];   y2 = f31(x2);  n2 = newton(x2,y2)  x3 = [‐1 ‐0.8 ‐0.6 ‐0.4 ‐0.2  0  0.2  0.4  0.6  0.8  1.0];   y3 = f31(x3);  n3 = newton(x3,y3)  xx = [‐1:0.02: 1]; %pham vi noi suy  yy = f31(xx); %ham thuc  yy1 = polyval(n1, xx); %ham xap xi qua 5 diem  yy2 = polyval(n2, xx); %ham xap xi qua 9 diem  yy3 = polyval(n3, xx); %ham xap xi qua 11 diem  subplot(221)  plot(xx, yy, ʹk‐ʹ,  xx, yy1, ʹbʹ)  subplot(224)  218 plot(xx, yy1‐yy, ʹrʹ, xx, yy2‐yy, ʹgʹ, xx, yy3‐yy,ʹbʹ) %do thi sai so  subplot(222)  plot(xx,yy,ʹk‐ʹ,  xx, yy2, ʹbʹ)  subplot(223)  plot(xx, yy, ʹk‐ʹ,  xx, yy3, ʹbʹ)  và nhận được kết quả.    Để tránh hiện tượng sai số lớn khi  số  điểm mốc  tăng  ta  dùng  nội  suy  nối  trơn(spline).  Trên  các  đoạn  nội  suy  ta  thay  hàm  bằng  một  đường  cong.  Các  đường cong này được ghép  trơn  tại các  điểm nối. Ta chọn các đường cong này là   hàm bậc 3  vì  hàm  bậc  1  và bậc hai khó   bảo đảm điều kiện nối trơn.     Cho một loạt giá trị nội suy (x1, y1),…,(xi, yi),…,(xn, yn). Trên mỗi đoạn ta  có một hàm bậc 3. Như vậy giữa nút i và (i +1) ta có hàm fi,i+1(x), nghĩa là ta  dùng (n ‐ 1) hàm bậc 3 f1,2(x), f2,3(x),…, fn‐1,n(x) để thay thế cho hàm thực. Hàm  fi,i+1(x) có dạng:  fi,i+1(x) = ai + bi(x ‐ xi) + ci(x ‐ xi)2 + di(x ‐ xi)3        (1)  Hàm này thoả mãn:  fi,i+1(xi) = ai = yi                  (3)    3 2i ,i 1 i 1 i i i i i i i i 1f (x ) d h c h b h a y+ + += + + + =           (4)    i ,i 1 i if (x ) b+′ =                    (5)  2 i ,i 1 i 1 i i i i if (x ) 3d h 2c h b+ +′ = + +               (6)  i ,i 1 i i if (x ) 2c y+′′ ′′= =                  (7)  i ,i 1 i 1 i i i i 1f (x ) 6d h 2c y+ + +′′ ′′= + =               (8)  Muốn nối trơn ta cần có đạo hàm bậc nhất liên tục và do đó:    i 1,i i i ,i 1 i if (x ) f (x ) k− +′′ ′′= =   Lúc này các giá  trị k chưa biết, ngoại  trừ k1 = kn = 0(ta các các mút  là điểm  uốn). Điểm xuất phát để tính các hệ số của fi,i+1(x) là biểu thức của  i ,i 1 if (x )+′′ . Sử  dụng nội suy Lagrange cho hai điểm ta có:    i ,i 1 i i i i 1 i 1f (x ) k L (x) k L (x)+ + +′′ = +   Trong đó:  y x xi‐1  xi  xi+1  yi‐1  yi+1 yi  fi‐1,i  fi,i+1  219   i 1 ii i 1 i i 1 i 1 i x x x xL (x) L (x) x x x x + + + + − −= =− −   Do vậy:    i i 1 i 1 ii ,i 1 i i i 1 k (x x ) k (x x )f (x ) x x + + + + − − −′′ = −   Tích phân biểu thức trên hai lần theo x ta có:  3 3 i i 1 i 1 i i ,i 1 i i 1 i i i 1 k (x x ) k (x x )f (x ) A(x x ) B(x x ) 6(x x ) + + + + + − − −= + − − −−   Trong đó A và B là các hằng số tích phân  Số hạng cuối trong phương trình trên thường được viết là Cx + D.   Đặt C = A ‐ B và D = ‐Axi+1 + Bxi để dễ dàng tính toán. Từ điều kiện fi,i+1(xi) = yi  ta có:  3 i i i 1 i i 1 i i i 1 k (x x ) A(x x ) y 6(x x ) + + + − + − =−   nên:    i i i i 1 i i 1 y k (x x )A x x 6 + + −= −−   Tương tự, điều kiện fi,i+1(xi+1) = yi+1 cho ta:    i 1 i 1 i i 1 i i 1 y k (x x )B x x 6 + + + + −= −−   Kết quả là:  3 i i 1 i ,i 1 i i 1 i i 1 i i 1 3 i 1 i i i i 1 i i 1 i i 1 i 1 i i i 1 k (x x )f (x ) (x x )(x x ) 6 x x k (x x ) (x x )(x x ) 6 x x y (x x ) y (x x ) x x + + + + + + + + + + + ⎡ ⎤−= − − −⎢ ⎥−⎣ ⎦ ⎡ ⎤−− − − −⎢ ⎥−⎣ ⎦ − − −+ − Đạo hàm cấp 2 ki tại các nút bên trong được tính từ điều kiện:    i 1,i i i ,i 1 if (x ) f (x )− +′ ′=   Sau khi biến đổi ta có phương trình:  i 1 i 1 i i i 1 i 1 i 1 i i 1 i 1 i i i 1 i 1 i i i 1 k (x x ) 2k (x x ) k (x x ) y y y y6 x x x x − − − + + + − + − + − + − + − ⎛ ⎞− −= −⎜ ⎟− −⎝ ⎠ Khi các điểm chia cách đều (xi+1 ‐ xi) = h ta  có:  220   ( )i 1 i i 1 i 1 i i 126k 4k k y 2y yh− + − ++ + = − +   i = 2, 3,…, n ‐ 1  Ta xây dựng hàm cubicspline() để nội suy:  function y = cubicspline(xData, yData, x)  %Ham nay xap xi bang da thuc bac 3 spline  %Cu phap: [yi,f] = cubicspline(xData, yData, x)  n = length(xData);  c = zeros(n‐1, 1); d = ones(n, 1);  e = zeros(n‐1, 1); k = zeros(n, 1);  c(1:n‐2) = xData(1:n‐2) ‐ xData(2:n‐1);  d(2:n‐1) = 2*(xData(1:n‐2) ‐ xData(3:n));  e(2:n‐1) = xData(2:n‐1) ‐ xData(3:n);  k(2:n‐1) = 6*(yData(1:n‐2) ‐ yData(2:n‐1))...  ./(xData(1:n‐2) ‐ xData(2:n‐1))...  ‐ 6*(yData(2:n‐1) ‐ yData(3:n))...  ./(xData(2:n‐1) ‐ xData(3:n));  [c, d, e] = band3(c, d e);  k = band3sol(c, d, e, k);  i = findseg(xData, x);  h = xData(i) ‐ xData(i+1);  y = ((x ‐ xData(i+1))^3/h ‐ (x ‐ xData(i+1))*h)*k(i)/6.0...  ‐ ((x ‐ xData(i))^3/h ‐ (x ‐ xData(i))*h)*k(i+1)/6.0...  + yData(i)*(x ‐ xData(i+1))/h...    ‐ yData(i+1)*(x ‐ xData(i))/h;  Ta có chương trình ctcubicspline.m dùng nội suy:  clear all, clc  x1 = 0:0.1:5;  y1 = (x1+1).^2;  while 1     x = input(ʹx = ʹ);      if isempty(x)         fprintf(ʹKet thucʹ);          break  221     end      y = cubicspline(xData, yData, x)      fprintf(ʹ\nʹ)  end  §5. NỘI SUY BẰNG ĐA THỨC CHEBYSHEV    Khi  nội  suy  bằng  đa  thức Newton  hay  Lagrange,  nghĩa  là  thay  hàm  thực bằng đa thức xấp xỉ, có khoảng cách cách đều thì sai số giữa đa thức nội  suy và hàm thực có xu hướng tăng tại hai mút nội suy. Ta thấy rõ điều này  khi chạy chương trình cttestintp.m.   Do vậy ta nên chọn các điểm mốc nội suy ở  hai mút dày hơn ở giữa. Một  trong những cách  chọn  phân  bố  các  điểm mốc  là  hình  chiếu  lên  trục  x  của  các  điểm  cách  đều  trên  đường  tròn  tâm tại điểm giữa của đoạn nội suy. Như vậy với  đoạn nội suy [‐1, 1] ta có:    k 2n 1 2kx cos 2(n 1) + −′ = π+   k = 1, 2,…,n          (1)  Với đoạn nội suy [a, b] bất kì:    k k b a b a b a 2n 1 2k a bx x cos 2 2 2 2(n 1) 2 − + − + − +′= + = π ++    k = 1, 2,…,n  (2)  Các nút nội suy này được gọi là các nút Chebyshev. Đa thức nội suy dựa trên  các nút Chebyschev gọi là đa thức nội suy Chebyshev.   Ta xét hàm thực:    2 1f(x) 1 8x = +   Ta chọn số nút nội suy  lần  lượt  là 5, 9, 11 và xây dựng các đa thức Newton  (hay Lagrange) c4(x), c8(x) và c10(x) đi qua các nút này và vẽ đồ  thị của hàm  thực cũng như sai số khi nội suy bằng chương trình ctcomchebynew.m với các  N khác nhau.  x1 = [‐1 ‐0.5 0 0.5 1.0];   y1 = f31(x1);  n1 = newton(x1,y1);  xx = [‐1:0.02: 1]; %pham vi noi suy  yy1 = polyval(n1,xx); %ham xap xi qua 5 diem  yy = f31(xx); %ham thuc  ‐1 11x′ 222 subplot(221)  plot(xx,yy,ʹk‐ʹ, x, y, ʹoʹ, xx, yy1, ʹbʹ);  title(ʹNewtonʹ)  subplot(223)  plot(xx, yy1‐yy, ʹrʹ) %do thi sai so  N = 4;   k = [0:N];  x = cos((2*N + 1 ‐ 2*k)*pi/2/(N + 1));  y = f31(x);  c = newton(x, y) %da thuc noi suy dua tren cac nut Chebyshev  xx = [‐1:0.02: 1]; %doan noi suy  yy = f31(xx); %do thi ham thuc  yy1 = polyval(c, xx); %do thi ham xap xi  subplot(222)  plot(xx, yy, ʹk‐ʹ,  x, y, ʹoʹ,  xx, yy1, ʹbʹ)  title(ʹChebyshevʹ)  subplot(224)  plot(xx, yy1‐yy, ʹrʹ) %do thi sai so  Khi tăng số điểm mốc, nghĩa là tăng bậc của đa thức Chebyschev, sai số giảm.  Đa thức Chebyshev bậc n được xác định bằng:    Tn+1(xʹ) = cos((n+1)arccos(xʹ))              (3)  và các nút Chebyshev cho bởi (1) là nghiệm của (3).  Ta có:  n 1 n n n 1 n 1 T (x ) cos(arccos(x ) narccos(x )) cos(arccos(x ))cos(narccos(x ) sin(arccos(x ))sin(narccos(x )) x T (x ) 0.5 cos((n 1)arccos(x ) cos((n 1)arccos(x ) x T (x ) 0.5T (x ) 0.5T (x ) + + − ′ ′ ′= + ′ ′ ′ ′= − ′ ′ ′ ′= + + − −⎡ ⎤⎣ ⎦ ′ ′ ′ ′= + − nên:    n 1 n n 1T (x ) 2xT (x ) T (x )+ −′ ′ ′= −     n ≥ 1          (4)  và  T0(xʹ) = 1    T1(xʹ) = cos(arccos(xʹ) = xʹ        (5)  Các đa thức Chebyshev đến bậc 6 là:    T0(x) = 1    T1(xʹ) = xʹ    T2(xʹ) = 2xʹ2 ‐ 1  223   T3(xʹ) = 4xʹ3 ‐ 3xʹ    T4(xʹ) = 8xʹ4 ‐  8ʹx2 + 1       T5(xʹ) = 16xʹ5 ‐ 20ʹx3 + 5xʹ    T6(xʹ) = 32xʹ6 ‐ 48xʹ4 + 18xʹ2 ‐ 1    T7(xʹ) = 64xʹ7 ‐ 112xʹ5 + 56xʹ3 ‐ 7xʹ  Hàm f(x) được xấp xỉ bằng:  N 2 a bm m x x b a 2m 0 f(x) d T (x ) +⎛ ⎞′= −⎜ ⎟− ⎝ ⎠= ′= ∑               (6)  Trong đó:  n n 0 k 0 k k k 0 k 0 1 1d f(x )T (x ) f(x ) n 1 n 1= = ′= =+ +∑ ∑           (7)  n m k m k k 0 n k k 0 2d f(x )T (x ) n 1 2 m(2n 1 2k)f(x )cos m 1,2,...,n n 1 2(n 1) = = ′= + + −= π =+ + ∑ ∑       (8)  Ta xây dựng hàm cheby() để tìm đa thức nội suy Chebyshev:   function [c, x, y] = cheby(f, N, a, b)  %vao : f = ten ham tren doan [a, b]  %Ra: c = Cac he so cua da thuc Newton bac N  % (x,y) = cac nut Chebyshev  if nargin == 2      a = ‐1;       b = 1;   end  k = [0: N];  theta = (2*N + 1 ‐ 2*k)*pi/(2*N + 2);  xn = cos(theta); %pt.(1)  x = (b ‐ a)/2*xn +(a + b)/2; %pt.(2)  y = feval(f,x);  d(1) = y*ones(N + 1,1)/(N+1);  for m = 2: N + 1    cos_mth = cos((m‐1)*theta);      d(m) = y*cos_mthʹ*2/(N + 1); %pt.(7)  end  xn = [2 ‐(a + b)]/(b ‐ a); %nghich dao cua t. (2)  224 T_0 = 1; T_1 = xn; %pt.(5)  c = d(1)*[0 T_0] +d(2)*T_1; %pt.(6)  for m = 3: N + 1      tmp = T_1;      T_1 = 2*conv(xn,T_1) ‐[0 0 T_0]; %pt.(4)      T_0 = tmp;      c = [0 c] + d(m)*T_1; %pt.(6)  end  Để  tìm  đa  thức Chebyshev dùng xấp xỉ hàm  2 1f(x) 1 8x = +   ta dùng  chương  trình ctcheby.m:  clear all, clc  N = 2;   a = ‐2;   b = 2;  [c, x1, y1] = cheby(ʹf31ʹ, N, a, b) %da thuc Chebyshev  %so sanh voi da thuc Lagrange/Newton   k = [0:N];   xn = cos((2*N + 1 ‐ 2*k)*pi/2/(N + 1));%pt.(1):nut Chebyshev  x = ((b‐a)*xn +a + b)/2; %pt.(2)  y = f31(x);   n = newton(x, y)   l = lagrange(x, y)  §6. XẤP XỈ HÀM BẰNG PHÂN THỨC HỮU TỈ    Xấp xỉ Padé dùng để xấp xỉ hàm f(x) tại x0 bằng hàm hữu tỉ:    m 0m,n 0 n 0 Q (x x )P (x x ) D (x x ) −− = −      2 m 0 1 0 2 0 m 0 2 n 1 0 2 0 n 0 q q (x x ) q (x x ) q (x x ) 1 d (x x ) d (x x ) d (x x ) + − + − + + −= + − + − + + − L L   (1)  với m = n hay m = n + 1  Trong đó f(x0), fʹ(x0),…, f(m+n)(x0) đã cho  Trước hết ta khai triển Taylor hàm f(x) tại x = x0 đến bậc (m + n).  225 m n 0 0 0 (m n) 2 m n0 0 0 0 2 m n 0 1 0 2 0 m n 0 f(x) T (x) f(x ) f (x )(x x ) f (x ) f (x )(x x ) (x x ) 2! (m n)! a a (x x ) a (x x ) a (x x ) + + + + + ′≈ = + − ′′+ − + + −+ = + − + − + + − L L       (2)  Để đơn giản ta coi x0 = 0. Ta cần tính các hệ số của Dn(x) và Qm(x) sao cho:    mm n n Q (x)T (x) 0 D (x)+ − =  hay Tm+n(x)Dn(n) ‐ Qm(x) = 0  nghĩa là:  m n n m 0 1 m n 1 n 0 1 m(a a x a x )(1 d x d x ) (q q x q x ) + ++ + + + + + = + + +L L L (3)  Cân bằng các số hạng cùng bậc ở hai vế ta có:  0 0 1 0 1 1 2 1 1 0 2 2 m m 1 1 m 2 2 m n n m a q a a d q a a d a d q a a d a d a d q− − − =⎧⎪ + =⎪⎪ + + =⎨⎪⎪ + + + + =⎪⎩ L L           (4)  m 1 m 1 m 1 2 m n 1 n m 2 m 1 1 m 2 m n 2 n m n m n 1 1 m n 2 2 m n a a d a d a d 0 a a d a d a d 0 a a d a d a d 0 + − − + + + − + + + − + − + + + + =⎧⎪ + + + + =⎪⎨⎪⎪ + + + + =⎩ L L L L          (5)  Trước hết ta giải (5) để tìm di và sau đó thay vào (4) để tìm qi.   Ta xây dựng hàm padeapp() để tính xấp xỉ:  function [num, den] = padeapp(f, xo, M, N, x0, xf)  %Vao : f = Ham can xap xi trong doan [xo, xf]  %Ra: num = Cac he so cua tu so  % den = Cac he so cua mau so  a(1) = feval(f, xo);  h = .01;   tmp = 1;  for i = 1:M + N      tmp = tmp*i*h; %i!h^i      dix = difapx(i, [‐i i])*feval(f, xo + [‐i:i]*h)ʹ; %dao ham      a(i + 1) = dix/tmp; %he so chuoi Taylor  226 end  for m = 1:N  n = 1:N;       A(m, n) = a(M + 1 + m ‐ n);      b(m) = ‐a(M + 1 + m);  end  d = A\bʹ; %pt.(5)  for m = 1: M + 1      mm = min(m ‐ 1,N);      q(m) = a(m:‐1:m ‐ mm)*[1; d(1:mm)];