CHƯƠNG 3: NỘI SUY VÀ XẤP XỈ HÀM
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 ...
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:
( x − x0 )...( x − xi −1 )( x − xi + 1 )...( x − x n )
Li =
( xi − x 0 )...( xi − xi −1 )( x i − x i + 1 )...( x i − x n )
Rõ ràng là Li(x) là một đa thức bậc n và :
⎧1 j=i
L i (x j ) = ⎨
⎩0 j≠i
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
Pn ( x) = ∑ f( x i )L i ( x)
i =0
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 y 1
Đa thức nội suy sẽ là :
P1(x) = yoL0(x) + y1L1(x1)
x − x1 x − x0
L0 = L 1 =
x 0 − x1 x1 − x 0
210
x − x1 x − x0
nên P1 ( x) = y 0 + y1
x 0 − x1 x1 − x 0
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 y 0 y1 y2
Đa thức nội suy sẽ là :
P2(x) = yoL0(x) + y1L1(x1) + y2L2(x2)
( x − x1 )( x − x 2 )
L0 =
( x 0 − x1 )( x0 − x 2 )
( x − x0 )( x − x 2 )
L1 =
( x1 − x0 )( x1 − x 2 )
( x − x 0 )( x − x1 )
L2 =
( x 2 − x 0 )( x 2 − x1 )
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
211
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à :
yi − y j
y[x i , x j ] =
xi − x j
Tỉ hiệu cấp hai của y tại xi, xj, xk là :
y[x i , x j ] − y[x j , x k ]
y[xi , x j , x k ] =
xi − xk
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 :
P ( x) − Pn ( x0 )
Pn [x , x0 ] = n
x − x0
là một đa thức bậc (n ‐ 1). Tỉ hiệu cấp 2 tại x, x0, x1 :
P [x , x0 ] − Pn [x0 , x1 ]
Pn [x , x 0 , x1 ] = n
x − x1
là một đa thức bậc (n‐2) v.v và tới tỉ hiệu cấp (n + 1) thì :
212
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ó:
∆y 0
y[x 0 , x 1 ] =
h
∆2 y 0
y[x 0 , x 1 , x 2 ] =
2h 2
...........
213
∆n y 0
y[x 0 , x 1 , x 2 ,..., x n ] =
n! h n
Bây giờ đặt x = x0 + ht trong đa thức Newton tiến ta được:
t( t − 1) 2 t( t − 1) ⋅ ⋅ ⋅ ( t − n + 1) n
Pn ( x 0 + ht) = y 0 + t∆y 0 + ∆ y0 + ⋅ ⋅ ⋅ + ∆ y0
2! n!
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ó:
t( t − 1) 2
Pn ( x 0 + ht) = y 0 + t∆y 0 + ∆ y0
2!
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:
t( t + 1) 2 t( t + 1) ⋅ ⋅ ⋅ ( t + n − 1) n
Pn ( x 0 + ht) = y n + t∇y n + ∇ yn + ⋅ ⋅ ⋅ + ∇ yn
2! n!
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
214
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:
y0 x0 − x
y x1 − x
P01 ( x) = 1
x1 − x 0
Đây là một đa thức bậc 1:
x − x1 x − x0
P01 ( x) = y 0 + y1
x 0 − x1 x1 − x 0
Khi x = x0 thì:
y0 x0 − x0
y1 x1 − x 0
P01 ( x 0 ) = = y0
x1 − x 0
Khi x = x1 thì:
y 0 x 0 − x1
y x1 − x1
P01 ( x1 ) = 1 = y1
x1 − x 0
Đa thức nội suy Lagrange của f(x) qua 3 điểm x0, x1, x2 có dạng:
215
P01 ( x) x0 − x
P ( x) x 2 − x
P012 ( x) = 12
x2 − x0
và là một đa thức bậc 2:
( x − x1 )( x − x 2 ) ( x − x 0 )( x − x 2 ) ( x − x 0 )( x − x1 )
P012 ( x) = y 0 + y1 + y2
( x 0 − x1 )( x 0 − x 2 ) ( x1 − x 0 )( x1 − x 2 ) ( x 2 − x 0 )( x 2 − x1 )
Khi x = x0 thì:
y0 x0 − x0
P ( x) x 2 − x 0
P012 ( x0 ) = 12 = y0
x 2 − x0
Khi x = x1 thì:
y 1 x 0 − x1
y x 2 − x1
P012 ( x1 ) = 1 = y1
x2 − x0
Khi x = x2 thì:
P01 ( x 2 ) x0 − x 2
y2 x2 − x2
P012 ( x 2 ) = = y2
x2 − x0
Tổng quát đa thức nội suy Lagrange qua n điểm là:
P01..( n −1) ( x) x 0 − x
P12..n ( x) x n − x
P012..n ( x) =
x2 − x0
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));
216
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:
1
f31(x) =
1 + 8x 2
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)
217
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ʹ)
y
và nhận được kết quả. fi‐1,i fi,i+1
Để 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 yi‐1 yi yi+1
đường cong này được ghép trơn tại các x
điểm nối. Ta chọn các đường cong này là xi‐1 xi xi+1
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)
fi ,i+1 (xi+1 ) = di h i + c i h i + bi h i + a i = y i+1
3 2
(4)
fi′,i+1 (xi ) = bi (5)
fi′,i+1 (xi+1 ) = 3di h i2 + 2c i h i + bi (6)
fi′′ +1 (x i ) = 2c i = y′′
,i i (7)
fi′′ +1 (xi+1 ) = 6di h i + 2c i = y′′+1
,i i (8)
Muốn nối trơn ta cần có đạo hàm bậc nhất liên tục và do đó:
fi′′ 1,i (x i ) = fi′′ +1 (x i ) = k i
− ,i
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 fi′′ +1 (xi ) . Sử
,i
dụng nội suy Lagrange cho hai điểm ta có:
fi′′ +1 (x i ) = k i L i (x) + k i+1L i+1 (x)
,i
Trong đó:
218
x − x i +1 x − xi
Li (x) = Li+1 (x) =
x i − x i +1 x i +1 − x i
Do vậy:
k i (x − x i+1 ) − k i+1 (x − x i )
fi′′ +1 (x i ) =
x i − x i +1
,i
Tích phân biểu thức trên hai lần theo x ta có:
k i (x − xi+1 )3 − k i+1 (x − xi )3
fi ,i+1 (xi ) = + A(x − xi+1 ) − B(x − xi )
6(xi − xi+1 )
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ó:
k i (xi − xi+1 )3
+ A(xi − xi+1 ) = y i
6(xi − x i+1 )
nên:
yi k (x − x i+1 )
A= − i i
x i − x i +1 6
Tương tự, điều kiện fi,i+1(xi+1) = yi+1 cho ta:
y i +1 k (x − xi+1 )
B= − i +1 i
x i − x i +1 6
Kết quả là:
k ⎡ (x − xi+1 )3 ⎤
fi ,i+1 (xi ) = i ⎢ − (x − xi+1 )(xi − xi+1 ) ⎥
6 ⎣ x i − x i +1 ⎦
k i+1 ⎡ (x − xi )3 ⎤
− − (x − x i )(xi − xi+1 )⎥
6 ⎢ x i − x i +1
⎣ ⎦
y i (x − xi+1 ) − y i+1 (x − xi )
+
x i − x i +1
Đạo hàm cấp 2 ki tại các nút bên trong được tính từ điều kiện:
fi′−1,i (x i ) = fi′,i+1 (x i )
Sau khi biến đổi ta có phương trình:
k i−1 (xi−1 − xi ) + 2k i (xi−1 − xi+1 ) + k i+1 (xi − xi+1 )
⎛ y − y i y i − y i +1 ⎞
= 6 ⎜ i −1 − ⎟
⎝ x i −1 − x i x i − x i + 1 ⎠
Khi các điểm chia cách đều (xi+1 ‐ xi) = h ta có:
219
6
k i−1 + 4k i + k i+1 = ( yi−1 − 2yi + yi+1 ) i = 2, 3,…, n ‐ 1
h2
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
220
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 ‐1 x′1 1
đoạn nội suy [‐1, 1] ta có:
2n + 1 − 2k
x′k = cos π k = 1, 2,…,n (1)
2(n + 1)
Với đoạn nội suy [a, b] bất kì:
b−a b+a b−a 2n + 1 − 2k a+b
xk = x′k + = cos π+ k = 1, 2,…,n (2)
2 2 2 2(n + 1) 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:
1
f(x) =
1 + 8x 2
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
221
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ó:
Tn +1 (x′) = cos(arccos(x′) + narccos(x′))
= cos(arccos(x′))cos(narccos(x′) − sin(arccos(x′))sin(narccos(x′))
= x′T n(x′) + 0.5 ⎡cos((n + 1)arccos(x′) − cos((n − 1)arccos(x′)⎤
⎣ ⎦
= x′T n(x′) + 0.5T n +1(x′) − 0.5T n −1(x′)
nên:
Tn +1(x′) = 2xT n(x′) − T n −1(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
222
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
f(x) = ∑ d m Tm (x′) x′= 2 ⎛ a+b ⎞
⎜ x− ⎟
(6)
m =0 b −a ⎝ 2 ⎠
Trong đó:
1 n 1 n
d0 = ∑ f(xk )T0 (x′k ) = n + 1 ∑ f(xk )
n + 1 k =0 k =0
(7)
2 n
dm = ∑ f(xk )Tm (x′k )
n + 1 k =0
(8)
2 n m(2n + 1 − 2k)
= ∑ f(xk )cos 2(n + 1) π m = 1,2,...,n
n + 1 k =0
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)
223
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
1
Để tìm đa thức Chebyshev dùng xấp xỉ hàm f(x) = ta dùng chương
1 + 8x 2
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ỉ:
Q (x − x 0 )
Pm ,n (x − x 0 ) = m
D n (x − x 0 )
q 0 + q 1 (x − x0 ) + q 2 (x − x0 )2 + L + q m (x − x0 )m
= (1)
1 + d1 (x − x0 ) + d 2 (x − x0 )2 + L + d n (x − x0 )n
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).
224
f(x) ≈ Tm + n (x) = f(x0 ) + f ′(x0 )(x − x0 )
f′′(x0 ) f (m + n) (x0 )
+ (x − x0 ) + L +
2
(x − x0 )m + n (2)
2! (m + n)!
= a 0 + a1(x − x0 ) + a 2(x − x0 )2 + L + a m + n(x − x0 )m + n
Để đơ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:
Q (x)
Tm + n (x) − m = 0 hay Tm+n(x)Dn(n) ‐ Qm(x) = 0
Dn (x)
nghĩa là:
(a 0 + a1x + L + a m + n x m + n )(1 + d1x + L + d n x n ) = (q 0 + q1x + L + q m x m ) (3)
Cân bằng các số hạng cùng bậc ở hai vế ta có:
⎧a 0 = q 0
⎪a + a d = q
⎪ 1
⎪
0 1 1
⎨a 2 + a1d1 + a 0d 2 = q 2 (4)
⎪L
⎪
⎪a m + a m −1d1 + a m −2d 2 + L + a m −nd n = q m
⎩
⎧a m +1 + a md1 + a m −1d 2 + L + a m −n+1d n = 0
⎪a
⎪ m + 2 + a m +1d1 + a md 2 + L + a m −n+ 2d n = 0
⎨ (5)
⎪L
⎪a m + n + a m + n −1d1 + a m + n−2d 2 + L + a md n = 0
⎩
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
225
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)]; %pt.(4)
end
num = q(M + 1:‐1:1)/d(N); den = [d(N:‐1:1)ʹ 1]/d(N); %giam dan
if nargout == 0 % ve ham thuc, khai trien taylor va ham Pade
if nargin §7. NỘI SUY BẰNG ĐA THỨC HERMIT
Trong một số trường hợp, ta cần tìm hàm đa thức không những đi qua
các điểm cho trước mà còn phải thoả mãn điều kiện về đạo hàm tại các điểm
đó. Ta gọi đa thức như vậy là đa thức nội suy Hermit. Để đơn giản, ta khảo
sát một đa thức bậc 3:
h(x) = H 3 x 3 + H 2 x 2 + H1x + H0 (1)
đi qua hai điểm (x0, y0), (x1, y1) và có các đạo hàm là y′ , y′ . Ta tìm các hệ số
0 1
Hi bằng cách giải hệ phương trình:
⎧h(x 0 ) = H 3 x0 + H 2 x0 + H1x0 + H0 = y 0
3 2
⎪
⎪h(x1 ) = H 3 x1 + H 2 x1 + H1x1 + H0 = y1
3 2
⎨ (2)
⎪ h′(x0 ) = 3H 3 x 0 + 2H 2 x0 + H1 = y′
2
0
⎪h′(x ) = 3H x 2 + 2H x + H = y′
⎩ 1 3 1 2 1 1 1
Các đạo hàm bậc nhất được tính gần đúng bằng:
h(x0 + ε) − h(x0 ) y 2 − y 0
y′ = =
ε ε
0
(3)
h(x1 ) − h(x1 − ε) y1 − y 3
y′ = =
ε ε
1
Bây giờ ta tìm đa thưc nội suy Lagrange hay Newton đi qua 4 điểm:
′
(x0, y0), (x 2 = x0 + ε ,y2 = y0 + y0 ε) , (x 3 = x1 − ε , y3 = y1 − y′ ε) , (x1, y1)
1
Hàm hermit() tạo nên phương trình (2):
function H = hermit(x0, y0, dy0, x1, y1, dy1)
A = [x0^3 x0^2 x0 1; x1^3 x1^2 x1 1;
3*x0^2 2*x0 1 0; 3*x1^2 2*x1 1 0];
b = [y0 y1 dy0 dy1]’; %Pt.(2)
H = (A\b)’;
Hàm hermits() dùng hàm hermit() để tính các hệ số của đa thức Hermit trên
nhiều đoạn và giá trị nội suy:
function [H,yi] = hermits(x, y, dy, xi)
% Tim cac he so cua c da thuc Hermite tren c doan
clc
for n = 1:length(x)‐1
H(n,:) = hermit(0, y(n), dy(n), x(n + 1)‐x(n), y(n + 1), dy(n + 1));
227
end
yi = ppval(mkpp(x, H),xi)
Để nội suy ta dùng chương trình cthermite.m:
clear all, clc
x = [0 1 2 3];
y = [1 2 4 5];
dy = [0 2 4 6];
[h, y] = hermits(x, y, dy, 1.5)
§8. BIẾN ĐỔI FOURIER
1. Biến đổi Fourrier: Tín hiệu thực tế thường bao gồm các thành phần có tần
số khác nhau. Chuỗi Fourier và phép bíến đổi Fourier là công cụ toán học
dùng để phân tích đặc tính tần số của tín hiệu. Có 4 định nghĩa tương tự nhau
về chuỗi và phép biến đổi Fourier, gồm: chuỗi Fourier liên tục theo t(CFS),
phép biến đổi Fourier liên tục theo t(CFT), chuỗi Fourier gián đoạn theo
t(DFS) và phép biến đổi Fourier gián đoạn theo t(DFT). Trong các công cụ
này, DFT dễ dàng lập trình trên máy tính nên trong phần này ta sẽ chú ý đến
nó.
Giả sử chuỗi số liệu { x[n] = x(nT), n = 0 : M ‐ 1} với T là chu kì lấy mẫu
có được bằng cách lấy mẫu một tín hiệu liên tục x(t) T lần trong một giây. N
cặp điểm DFT và iDFT được định nghĩa bằng:
N −1
DFT: X(k) = ∑ x[n]e ‐j2πnk/N (1a)
n =0
N −1
1
iDFT: x[n] = ∑ X(k)e j2πnk/N
N n =0
(1b)
Nói chung hệ số DFT của X(k) là một số phức và nó xác định biên độ và pha
của thành phần tín hiệu có tần số số Ωk = kΩ0(rad), tương ứng với tần số
tương tự ωk = kω0 = kΩ0/T = 2πk/NT (rad/s). Ta gọi Ω0 = 2π/N và ω0 = 2π/NT là
các tần số cơ bản số và tương tự (tần số phân giải) vì đây là hiệu tần số có thể
phân biệt bởi N điểm DFT.
DFT và DFS có cùng bản chất nhưng khác nhau về phạm vi thời
gian/tần số. Cụ thể là tín hiệu x[n] và DFT X[k] của nó kéo dài hữu hạn trên
phạm vi thời gian/tần số {0 ≤ n ≤ N‐1} và {0 ≤ k ≤ N‐1}. Tín hiệu x[n] được
228
phân tích bởi DFS và DFS của nó X(k) là chu kì tín hiệu với chu kì N trên toàn
bộ tập số nguyên.
Biến đổi Fourier nhanh FFT là thuật toán hiệu quả để tính DFT và iDFT
được xây dựng bằng cách dùng tính chu kì và tính đối xứng cuả nhân tử
ei2πnk/N để giảm bớt số nhân tử phức từ N2 thành (N/2)log2N )N thể hiện kích
thước của DFT. Hàm MATLAB fft() và ifft() thực hiện thuật toán đối với N =
2l (l là số nguyên không âm). Nếu độ dài M của chuỗi số liệu ban đầu không
phải là bội số của 2, có thể mở rộng bằng cách đệm thêm số 0 vào cuối chuỗi
và gọi là đệm zero.
Ta xem xét hiệu qủa này bằng cách thực hiện đoạn lệnh trong
ctcompdftfft.m.
%So sanh phep bien doi Fourier nhanh va roi rac
clear, clf
N = 2^10;
n = [0:N ‐ 1];
x = cos(2*pi*200/N*n)+ 0.5*sin(2*pi*300/N*n);
tic %ngung dong ho
for k = 0:N ‐ 1
X(k+1) = x*exp(‐j*2*pi*k*n/N).ʹ;
end %DFT
k = [0:N ‐ 1];
for n = 0:N ‐ 1
xr(n + 1) = X*exp(j*2*pi*k*n/N).ʹ;
end %IDFT
time_dft = toc
plot(k,abs(X))
pause, hold on
tic
X1 = fft(x); %FFT
xr1 = ifft(X1); %IFFT
time_fft = toc %dua ra thoi gian thuc hien
clf, plot(k,abs(X1),ʹrʹ) %pho bien do
Chạy đoạn lệnh và so sánh thời gian thực hiện 1024 điểm tính DFT/iDFT và
FFT/iFFT.
229
