Bài giảng Phương pháp số - Chương 3: Phép nội suy và hồi quy

Chương 3: Phép nội suy và hồi quy 
CHƯƠNG 3 
PHÉP NỘI SUY VÀ HỒI QUY 
MỤC ĐÍCH, YÊU CẦU 
Sau khi học xong chương 3, yêu cầu sinh viên: 
1. Hiểu được thế nào là bài toán nội suy và hồi quy. 
2. Nắm được các phương pháp nội suy đa thức, biết cách tìm các đa thức nội suy theo các 
phương pháp đó. 
3. Biết được khớp đường cong - Nội suy Spline là gì? 
4. Nắm và giải được các bài toán bằng phương pháp bình phương tối thiểu 
5. Biết cách đánh giá sai số của từng phương pháp. 
3.1. MỞ ĐẦU 
Thông thường trong một số lĩnh vực như kinh tế chẳng hạn, các đại lượng khảo sát thường 
không được cho dưới dạng hàm liên tục, mà là bảng các giá trị rời rạc. Các phương pháp giải tích 
toán học thường tính toán với các hàm cho bởi các công thức, do đó không thể áp dụng trực tiếp 
để nghiên cứu các hàm cho dưới dạng rời rạc như thế này. Cũng có khi ta biết rằng đại lượng y là 
một hàm của đại lượng x, tức là y = f(x), nhưng ta không biết biểu thức hàm f(x) mà chỉ biết một 
số giá trị yi tương ứng với các giá trị của x tại các điểm xi như trong bảng sau: 
x x
0
x
1
x
2
. 
. . 
x
n-1
x
n
y y
0
y
1
y
2
. 
. . 
y
n-1
y
n
Thông thường thì x0 < x1 < x2 < . . . < xn và các điểm này có thể phân bố cách đều hoặc 
không. Mặc dầu ta chỉ biết các giá trị của y tại các điểm mốc xi, nhưng trong nhiều trường hợp ta 
cần tính toán với các giá trị y tại các vị trí khác của x. Một câu hỏi đặt ra là: cho một điểm x 
không thuộc các điểm xi cho ở trên, làm thế nào chúng ta có thể tính được giá trị y tương ứng với 
nó, sao cho chúng ta có thể tận dụng tối đa các thông tin đã có? 
Bài toán nội suy là bài toán tìm giá trị gần đúng của y tại các điểm nằm giữa các giá trị x 
không có trong bảng trên. Nếu cần tìm các giá trị gần đúng của y tại các điểm x nằm ngoài khoảng 
[x0,xn] thì bài toán được gọi là bài toán ngoại suy. Một bộ n+1 cặp các giá trị đã biết của x và y: 
(x0,y0), (x1,y1), . . . ,(xn,yn) được gọi là một mẫu quan sát, còn x0, x1, ... , xn được gọi là các điểm 
quan sát và y0, y1, ... , yn là các kết quả quan sát. 
 42 
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Chương 3: Phép nội suy và hồi quy 
Vì bài toán của chúng ta không chỉ giải quyết với một giá trị x cụ thể, mà là cả một miềm 
giá trị nào đó của x. Do đó câu hỏi trên cũng tương đương với vấn đề sau: hãy tìm một hàm g(x) 
sao cho miền giá trị của nó chứa các điểm (x0, x1, ..., xn) và hàm này xấp xỉ tốt nhất tập số liệu đã 
có là các cặp (x0,y0), (x1,y1), ..., (xn,yn) theo một nghĩa nào đó. Chúng ta thấy ngay là tập số liệu là 
hữu hạn, còn tập các giá trị cần ước lượng là vô hạn, nên sẽ có vô số hàm g(x) nếu chúng ta không 
đưa ra một số ràng buộc nào đó về g(x). Điều đầu tiên chúng ta quan tâm là nên chọn dạng hàm 
g(x) như thế nào. 
Một cách tự nhiên, ta có thể đặt điều kiện về hàm g(x) như sau: 
• 
• 
• 
g(xi) i =0,1,2,...,n gần các điểm yi nhất theo một nghĩa nào đó. 
g(x) là duy nhất theo một số điều kiện nào đó. 
Hàm g(x) liên tục, không có điểm gấp khúc và ít thay đổi trong từng đoạn [xi,xi+1]. 
Các định lý về xấp xỉ sau đây của Weierstrass sẽ cho chúng ta gợi ý về dạng hàm của g(x). 
Định lý Weierstrass 1 về xấp xỉ hàm. 
Cho f (x) là một hàm thực liên tục xác định trên khoảng [a,b]. Khi đó với mọi ε>0 tồn tại 
một đa thức p(x) bậc m với các hệ số thực sao cho với mọi giá trị x∈[a,b] ta có |f(x) - p(x)|<ε. 
Định lý Weierstrass 2 về xấp xỉ hàm. 
Cho f (x) là một hàm thực liên tục xác định trên khoảng [-π,π] và f(-π) = f(π). Khi đó với 
mọi ε>0 tồn tại một đa thức lượng giác 
qm(x) = 2
0a + ∑ [a
=
m
j 1
j cos(jx) + bj sin(jx)] 
với các hệ số thực sao cho với mọi giá trị x∈[-π,π] ta có |f(x) - q(x)|<ε. 
Từ các định lý trên đây ta thấy rằng chọn đa thức là thích hợp cho dạng hàm g(x). Đa thức 
là hàm quen thuộc và ta đã biết nhiều tính chất của nó. 
Người ta thường dùng các phương pháp xấp xỉ sau để xác định đa thức p(x): 
1. Nếu ta biết rằng các cặp giá trị (x0,y0), (x1,y1), ..., (xn,yn) là thể hiện của một hàm f(x) nào 
đó, tức là ta biết rằng y=f(x) và như vậy tại các điểm xi, i=0,1,...,n yi = f(xi). Trong trường 
hợp này ta đòi hỏi đa thức p(x) phải đi qua các điểm (xi,yi), i=0,1,...,n. 
Bài toán nội suy bây giờ có thể phát biểu cụ thể hơn như sau: 
Cho một mẫu quan sát gồm n+1 cặp các giá trị đã biết của x và y : (x0,y0), (x1,y1), . 
. . ,(xn,yn) . Hãy xây dựng một đa thức bậc m ≤ n 
 pm(x) = a0 + a1x1 + . . . am-1xm-1 + amxm (3.1) 
sao cho pm(xi) = yi , i = 0, 1,..., n (3.2) 
Người ta gọi bài toán trên đây là bài toán nội suy đa thức, và đa thức pm(x) được gọi 
là đa thức nội suy. 
Trong một số ứng dụng vật lý ta gặp các hiện tượng có tính chất tuần hoàn. Khi đó đa 
thức lượng giác tỏ ra thích hợp hơn trong bài toán nội suy. Và trong bài toán trên đa thức 
pm(x) được thay bằng đa thức lượng giác 
 43
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Chương 3: Phép nội suy và hồi quy 
 qm(x) = 2
0a + ∑ [a
=
m
j 1
j cos(jx) + bj sin(jx)] 
2. Nội suy trong trường hợp số đo không hoàn toàn chính xác: 
Trong thực tế các giá trị yi tại các điểm quan sát lại thường chỉ là các giá trị gần đúng 
của các giá trị thật. Nói cách khác thực ra ta chỉ có yi ≈ f(xi) mà thôi. Trong trường hợp này 
nếu ta áp đặt điều kiện về đa thức nội suy phải thỏa mãn pm(xi) = yi thì không hợp lý. Thay 
vì tìm một đa thức thỏa mãn điều kiện này, ta tìm đa thức pm(x) = a0 + a1x1 + . . . am-1xm-1 
+ amxm, tức là xác định các hệ số a0, a1, . . ,am sao cho tổng bình phương sai số là bé nhất, 
tức là 
 e = ∑ (y
=
n
i 0
i -∑ a
=
m
j 0
j xij )2
là bé nhất. Phương pháp nội suy theo tiêu chuẩn này được gọi là phương pháp bình 
phương bé nhất hay là phương pháp bình phương cực tiểu. 
Ngoài hai phương pháp thông dụng trên, người ta còn dùng phương pháp xấp xỉ Csebisev 
dựa trên tiêu chuẩn: 
ni≤≤0
max |yi - p(xi)| 
cực tiểu. 
3.2. NỘI SUY ĐA THỨC 
3.2.1. Sự duy nhất của đa thức nội suy 
Ta có định lý sau đây: 
Định lý. Có duy nhất một đa thức có bậc không quá n và đi qua n+1 điểm cho trước (x0,y0), 
(x1,y1), . . . ,(xn,yn) . 
Chứng minh. Ta xét đa thức có dạng (3.1) trên đây và thỏa mãn (3.2). Kết hợp (3.1) và 
(3.2) ta có 
⎟⎟
⎟⎟
⎟
⎠
⎞
⎜⎜
⎜⎜
⎜
⎝
⎛
ny
y
y
.
1
0
= (3.3) 
⎥⎥
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢⎢
⎢
⎣
⎡
n
nnnn
n
n
xxxx
xxx
xxx
22
1
2
11
0
2
00
1
.......
...1
...1
⎟⎟
⎟⎟
⎟
⎠
⎞
⎜⎜
⎜⎜
⎜
⎝
⎛
na
a
a
.
1
0
Hay có thể biểu diễn gọn hơn dưới dạng ma trận 
 Y = V a 
 Trong đó 
V = 
⎥⎥
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢⎢
⎢
⎣
⎡
n
nnnn
n
n
xxxx
xxx
xxx
22
1
2
11
0
2
00
1
.......
...1
...1
 44 
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Chương 3: Phép nội suy và hồi quy 
chính là ma trận Vandermon, ta có 
 det V = (x∏
≤<≤ nji0
j - xi) 
Vì ta đã giả thiết các điểm xi và xj khác nhau, do đó ma trận này khác 0 nên hệ phương trình 
(3.3) có nghiệm duy nhất cho các ai, và như vậy đa thức pn(x) được xác định duy nhất. (Nếu khi 
giải phương trình (3.3) mà ta nhận được an ≠ 0 thì đa thức này có bậc là n). 
3.2.2. Tính giá trị đa thức bằng phương pháp Horner 
Trong bài toán nội suy đa thức ta sẽ phải thường xuyên tính giá trị của đa thức pm(x)= a0 + 
a1x1 + . . . am-1xm-1 + amxm tại điểm x. Nếu tính trực tiếp ta phải thực hiện khá nhiều phép tính, và 
khi tính các giá trị mũ của x có thể ta gặp các giá trị lớn, cho dù trong thực tế các thành phần của 
đa thức triệt tiêu lẫn nhau và giá trị của đa thức không lớn. Horner đưa ra cách tính sau loại trừ 
được các nhược điểm trên. 
Ta viết lại đa thức pm(x) dưới một dạng khác: 
pm(x) = amxm + am-1xm-1 + . . . + a1x1 + a0 = (...((amx + am-1)x + am-2)x+...+a1)x+a0 
Từ đây ta có cách tính pm(x) trên máy tính như sau: 
Đặt Pm = am
 Pm-1 = Pmx+am-1 
 ... 
 Pi = Pi-1x + ai
Khi tính toán ta không cần lưu trữ tất c các giá trị của Pi, mà chỉ cần lưu trữ các giá trị của 
Pi trong một vị trí bộ nhớ. Thuật toán trên trở thành: 
Đặt P = am 
Cho i chạy từ m-1 đến 0, tức là i=m-1,m-2,...,0 
Đặt P = Px + ai
P cuối cùng tính được chính là giá trị của đa thức tại x. 
3.2.3. Sai số của đa thức nội suy 
Định lý Rolle: Cho f(x) là hàm số thực liên tục trên khoảng đóng [a,b] và khả vi trên 
khoảng mở (a,b) và f(a) = f(b). Khi đó tồn tại điểm ξ∈ (a,b) sao cho f'(ξ) = 0. 
Định lý: Giả sử hàm f(x) có đạo hàm liên tục đến cấp n+1 trên đoạn [a,b] và pm(x) là đa thức 
nội suy, tức là:pm(xi) = f(xi) = yi , i = 0, 1,..., n. Với các mốc nội suy là a = x0 < x1 <... < xn = b. 
Đặt ωn+1(x) = và R(x) = f(x) - p)(
0
i
n
i
xx∏
=
− m(x) 
Khi đó với ∀ x ∈ [a,b] tồn tại η∈ [a,b] (phụ thuộc vào x) sao cho 
 R(x) = 
)!1(
)()1(
+
+
n
f n η ωn+1(x) (3.4) 
 45
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Chương 3: Phép nội suy và hồi quy 
Hệ quả. Gọi M = |f
bxa ≤≤
sup (n+1) (x)| khi đó ta có 
 |R(x)| ≤ | f(x) - pm(x)| ≤ )!1( +n
M | ωn+1(x) | (3.5) 
đây là công thức đánh giá sai số của đa thức nội suy. 
Ta có thể áp dụng hệ quả (3.5) để đánh giá sai số đa thức nội suy. 
Ví dụ. Cho bảng giá trị của hàm số y = sinx 
x 0 
4
π 
2
π 
y 0 0.707 1 
Hãy đánh giá sai số khi dùng đa thức nội suy để tính gần đúng sin 
3
π . 
Giải. Bài ra không đặt vấn đề tính xấp xỉ sin 
3
π mà chỉ yêu cầu tính sai số. 
Ta có n = 2 và như vậy M = |sin
bxa ≤≤
sup (n+1) (x)| =1, do đó 
|R2(
3
π )| ≤ 
!3
1
3
π
12
π
6
π = 0.024 
Sau đây ta sẽ xét một số phương pháp tìm đa thức nội suy dựa vào các điểm mốc cách đều 
và không cách đều. 
3.2.4. Phương pháp nội suy Lagrange 
Giả sử ta có các điểm quan sát x0, x1, ... xn với khoảng chia đều hoặc không đều và một 
dãy các giá trị quan sát y0, y1, ... yn . 
Ý tưởng đơn giản đầu tiên là tìm một đa thức nội suy có bậc n (chính xác hơn là có bậc 
không quá n) sao cho trong đó các cặp (xi,yi) i = 0,1, ..., n có vai trò bình đẳng. Thí dụ ta tìm 
pn(x) có dạng: 
pn(x) = H0(x) + H1(x) + . . . + Hn(x) 
Các hàm Hi(x) đều có bậc không quá n và Hi(xi) = yi, Hki(xj) = 0 khi j≠i. Để Hi(xj) = 0 khi 
j≠i thì Hi(x) có dạng: 
Hi(x) = K(x)(x-x0) (x-x1)... (x-xi-1) (x-xi+1)... (x-xn) 
Từ điều kiện Hi(xi) = yi ta có 
K(x)(x-x0) (x-x1)... (x-xi-1) (x-xi+1)... (x-xn) = yi 
Suy ra 
K(x) =yi )x-(x )...x-)(xx-(x )...x-(x )x-(x
)x-(x )...x-)(xx-(x )...x-(x )x-(x
ni1ii1-ii2i1i
n1i1-i21
+
+ 
 46 
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Chương 3: Phép nội suy và hồi quy 
Nếu ta ký hiệu Li(x) = )x-(x )...x-)(xx-(x )...x-(x )x-(x
)x-(x )...x-)(xx-(x )...x-(x )x-(x
ni1ii1-ii2i1i
n1i1-i21
+
+ 
Ta nhận thấy đa thức Li(x) có tính chất 
 Li(xj) = ⎩⎨
⎧
≠
=
ij
ij
0
1
và đa thức pn(x) có dạng 
pn(x) = y0L0(x) + y1L1(x) + . . . + ynLn(x) (3.6) 
Như vậy ta có 
L0(x) = 
)x-(x )...x-(x )x-(x
)x-(x )...x-(x )x-(x
n02010
n21 
L1(x) = )x-(x )...x-(x )x-(x
)x-(x )...x-(x )x-(x
n12101
n20 
. . . 
Li(x) = )x-(x )...x-)(xx-(x )...x-(x )x-(x
)x-(x )...x-)(xx-(x )...x-(x )x-(x
ni1ii1-ii2i1i
n1i1-i21
+
+ 
. . . 
Ln(x) = )x-(x )...x-(x )x-(x
)x-(x )...x-(x )x-(x
1-nn2n0n
1-n10 
Đa thức nôi suy được xây dựng theo cách trên đây được gọi là đa thức nội suy Lagrange. 
Ví dụ1 :Với hàm số y=sin(x/2) tại các nút giá trị sau: 
I xi yi
0 0 0.000 
1 1.5 0.682 
2 2 0.841 
Hãy xác định đa thức nội suy Lagrange đi qua các điểm trên? Hãy tính giá trị gần đúng của 
hàm số tại điểm x=1? Hãy đánh giá sai số lý thuyết tại x=1 
Theo phần lý thuyết trên đa thức nội suy Lagrange đi qua các điểm (xi,yi) được xác định 
như sau 
P(x)=∑ y
=
n
i 0
iLi
Với Li=( (x-x∏
≠=
n
ij
j 0
j))/(∏ (x
≠=
n
ij
j 0
i-xj)) 
 47
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Chương 3: Phép nội suy và hồi quy 
Ta có: 
L0(x)=(x-1.5)(x-2)/3 
L1(x)=-4/3x(x-2) 
L2(x)=x(x-1.5) 
Vậy P(x)= y0L0(x)+y1L1(x)+y2L2(x)=0- 0.682*4/3x(x-2) + 0.841*x(x-1.5) 
Vậy P(1)=0- 0.682*4/3(1-2) + 0.841(1-1.5) 
 P(1)=0.4888 
* Đánh giá sai số lý thuyết: 
R(x)=(f(3)(η)/3!)*x(x-1.5)(x-2). 
f(x)=sin(x/2). vËy f(3)(x)=(-1/2(3))cos(x/2) 
Ta có |f(3)(x)|≤1/8 
Vậy R(1)≤(f(3)(x)/3!)*(1-1.5)(1-2)≈0.01042. 
Ví dụ 2.Với hàm số y=sin(x/3) tại các nút giá trị sau: 
I xi yi
0 0 0.000 
1 1.5 0.479 
2 2 0.618 
Hãy xác định đa thức nội suy Lagrange đi qua các điểm trên? Hãy tính giá trị gần đúng của 
hàm số tại điểm x=1? Hãy đánh giá sai số lý thuyết tại x=1 
Theo phần lý thuyết trên đa thức nội suy Lagrange đi qua các điểm (xi,yi) được xác định 
như sau 
P(x)=∑ y
=
n
i 0
iLi
Với Li=(∏ (x-x
≠=
n
ij
j 0
j))/(∏ (x
≠=
n
ij
j 0
i-xj)) 
 48 
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Chương 3: Phép nội suy và hồi quy 
Ta có: 
L0(x)=(x-1.5)(x-2)/3 
L1(x)=-4/3x(x-2) 
L2(x)=x(x-1.5) 
Vậy P(x)= y0L0(x)+y1L1(x)+y2L2(x)=0- 0.479*4/3x(x-2) + 0.618*x(x-1.5) 
Vậy P(1)=0- 0.479*4/3(1-2) + 0.618(1-1.5) 
 P(1)≈0.3297 
* Đánh giá sai số lý thuyết: 
R(x)=(f(3)(η)/3!)*x(x-1.5)(x-2). 
f(x)=sin(x/3). vậy f(3)(x)=(-1/3(3))cos(x/3) 
Ta có |f(3)(x)|≤1/27 
Vậy R(1)≤( |f(3)(x)|/3!)*(1-1.5)(1-2)≈((1/27)/6)*0.5=0.00386 
3.2.5. Sai phân 
Cách xây dựng đa thức nội suy Lagrange khá đơn giản về mặt ý tưởng. Tuy nhiên nhược 
điểm của nó là mỗi lần bổ sung thêm một số điểm quan sát mới ta lại phải tính lại từ đầu. Người ta 
tìm cách xây dựng một đa thức nội suy sao cho khi bổ sung các điểm quan sát thì ta không phải 
tính lại phần đa thức đã có. Thí dụ từ các điểm quan sát (x0,y0), (x1,y1),..., (xk,yk) ta tính được đa 
thức pk(x). Khi bổ sung thêm các điểm (xk+1,yk+1),..., (xn,yn) thì đa thức nội suy tương ứng với mẫu 
quan sát (x0,y0),..., (xn,yn) sẽ có dạng pn(x) = pk(x) + u(x). 
Để thực hiện và trình bày điều này một cách rõ ràng, sáng sủa, trước hết ta cần đến khái 
niệm sai phân như sau: 
a. Định nghĩa: 
Cho f(x) là hàm của x và h = Δx là một hằng số không đổi biểu thị cho khoảng thay đổi trên 
biến x và được gọi là số gia của x. Khi đó số gia tương ứng trên f(x): 
Δf(x) = f(x+Δx) - f(x) (3.7) 
được gọi là sai phân tiến cấp một tại điểm x của f(x) tương ứng với h. Gia số được tính bởi 
 ∇f(x) = f(x) - f(x-Δx) (3.8) 
được gọi là sai phân lùi cấp một tại điểm x của f(x) tương ứng với h. 
Vì sai phân tiến g(x) của một hàm lại là một hàm của x do đó ta lại có thể định nghĩa sai 
phân tiến của g(x). Khi đó ta gọi sai phân tiến cấp một của g(x) là sai phân tiến cấp 2 của f(x), và 
cứ như vậy ta có thể định nghĩa sai phân tiến cấp k của một hàm f(x). 
 49
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Chương 3: Phép nội suy và hồi quy 
Với sai phân lùi ta cũng có lập luận và định nghĩa tương tự. 
b. Sai phân tiến 
Giả sử các điểm x0, x1, ... xn thoả mãn điều kiện 
xi+1 - xi = h 
 yi = f(xi), i = 0, 1, ... (3.9) 
Ta có thể thấy rằng sai phân tiến 
Δ2 yi = Δ(Δyi) = Δyi+1 - Δyi = yi+2 - yi+1 - (yi+1 - yi) = yi+2 - 2yi+1 +yi
Tổng quát ta có thể chứng minh rằng 
Δk yi = ∑ (-1)
=
k
j 0
j Ckj yi+(k-j) (3.10) 
Bảng các sai phân tiến 
x y Δy Δ2y Δ3y Δ4y 
x0 y0 Δy0 Δ2y0 Δ3y0 Δ4y0
x1 y1 Δy1 Δ2y1 Δ3y1 Δ4y1
x2 y2 Δy2 Δ2y2 Δ3y2 Δ4y2
x3 y3 Δy3 Δ2y3 Δ3y3 Δ4y3
x4 y4 Δy4 Δ2y4 Δ3y4 Δ4y4
. . . . . . 
 Δ4yn-5
 Δ3yn-4 Δ4yn-4
 Δ2yn-3 Δ3yn-3 
. . Δyn-2 Δ2yn-2 
xn-1 yn-1 Δyn-1 
xn yn 
c. Sai phân lùi 
Với sai phân lùi ta có 
∇2 yi = ∇ (∇yi) = ∇yi - ∇yi-1 = yi - yi-1 - (yi-1 - yi-2) = yi - 2yi-1 +yi-2
Tổng quát ta có thể chứng minh rằng 
 ∇k yi  ... da thuc noi suy Newton, xqs[i] la 
cac diem quan sat*/ 
double polinew(double x) 
 {int i;double mx,s; 
 s=anew[0]; 
 mx=1; 
 for(i=0;i<nqs;i++) 
 {mx*=(x-xqs[i]); 
 s+=anew[i+1]*mx; 
 } 
 return s; 
 } 
//=============================================== 
/*Tra ve gia tri da thuc noi suy tai x; alag[i] la cac he so cua 
 da thuc noi suy Lagrange , xqs[i] la cac diem quan sat*/ 
double polilag(double x) 
 {int i,j;double mx,s; 
 s=0; 
 for(i=0;i<=nqs;i++) 
 {mx=1; 
 for(j=0;j<=nqs;j++) 
 if(j!=i) mx*=(x-xqs[j]); 
 s+=alag[i]*mx; 
 } 
 return s; 
 } 
//=============================================== 
/*Noi suy Newton voi khoang chia khong can deu */ 
void nsnewton(double *a) 
 58 
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Chương 3: Phép nội suy và hồi quy 
 { int h,i,j,k;double tmp;kvecto sp; 
 a[0]=yqs[0]; 
 for(i=0;i<=nqs-1;i++) //Tinh sai phan bac 1; 
 sp[i]=(yqs[i+1]-yqs[i])/(xqs[i+1]-xqs[i]); 
 a[1]=sp[0]; 
 for(k=2;k<=nqs;k++)//Tinh cac sai phan bac k va cac he so a[i] 
 { for(i=0;i<=nqs-k;i++) 
 sp[i]=(sp[i+1]-sp[i])/(xqs[i+k]-xqs[i]); 
 a[k]=sp[0]; 
 } 
 } 
//=============================================== 
/*Noi suy Lagrange voi khoang chia khong can deu 
 Tinh cac he so*/ 
 a[i] = 1/((x[i]-x[0])(x[i]-x[1])...(x[i]-x[i-1])(x[i]-x[i+1])...x[m])*/ 
void nslagrange(double *a) 
 { int h,i,j,k;double tmp; 
 for(i=0;i<=nqs;i++) //Tinh cac he so a[i]; 
 { tmp=1; 
 for(j=0;j<=nqs;j++) if(j!=i) tmp=tmp*(xqs[i]-xqs[j]); 
 a[i]=yqs[i]/tmp; 
 } 
 } 
3.3. KHỚP ĐƯỜNG CONG - NỘI SUY SPLINE 
Trong phần trước ta đã xét bài toán nội suy dùng đa thức và như ta đã thấy, các đa thức nội 
suy thường có bậc là n, trong đó n+1 là số điểm quan sát. Ta có thể nội suy bằng đa thức bậc m 
nhỏ hơn n, nhưng như vậy thì ta cũng chỉ dùng đến mẫu quan sát dựa trên m+1 điểm là (x0,y0), 
(x1,y1), . . . ,(xm,ym) và như thế chỉ nội suy được giá trị của hàm tại các điểm x ∈ [x0,xm] . Điều 
này tỏ ra không được phù hợp với thực tế cho lắm. Thật vậy, giả sử trong thực tế hàm f(x) là 
một đa thức bậc 3 nhưng vì ta không biết điều này nên phải dùng đa thức nội suy. Theo một cách 
tự nhiên, ta nghĩ rằng nếu có càng nhiều thông tin thì ta càng giải quyết bài toán tốt hơn. Nghĩa là 
nếu có càng nhiều điểm quan sát thì kết quả của chúng ta càng gần với thực tế hơn. Tuy nhiên nếu 
dùng đa thức nội suy như kiểu chúng ta vừa khảo sát thì không có được như điều chúng ta mong 
đợi. Mặc dầu dạng thật của đa thức là bậc 3, nhưng nếu dùng 5 điểm quan sát thì ta phải tính các 
hệ số đa thức bậc 4, 10 điểm thì ta phải tính toán với đa thức bậc 9,...nghĩa là càng dùng nhiều 
điểm thì ta càng đi xa thực tế hơn. Phép nội suy đa thức còn có một nhược điểm nữa là số lượng 
 59
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Chương 3: Phép nội suy và hồi quy 
phép tính cần thực hiện phụ thuộc rất nhiều vào cỡ của mẫu quan sát. Trong kỹ thuật truyền thông 
chẳng hạn, việc chuyển đổi một tín hiệu số có hàng ngàn điểm quan sát sang dạng tương tự là vấn 
đề thường gặp. Thế nhưng chỉ cần nội suy đa thức cho 101 điểm quan sát ta đã phải dùng đến đa 
thức bậc 100, và việc dùng đa thức bậc 100 để tính toán cho các điểm còn lại là một việc tiêu tốn 
tài nguyên máy một cách quá lãng phí. Vì vậy có thể nói rằng phép nội suy đa thức chỉ có ý nghĩa 
lý thuyết mà thôi, trong thực tế hầu như người ta không dùng đến. 
Để tìm kiếm một cách nội suy gần với thực tế hơn, chúng ta hãy bắt đầu bằng một thao tác 
đơn giản mà chúng ta hay thực hiện hồi còn học phổ thông. Khi vẽ một đồ thị hàm số nào đó, đầu 
tiên ta vẽ các điểm rời rạc, và vẽ được càng nhiều điểm càng tốt. Sau đó ta dùng bút nối các điểm 
đó với nhau, nhưng ta không nối bằng thước kẻ, mà nối bằng bút và sự quan sát bằng mắt sao cho 
các đoạn nối các điểm thành một đường mịn, không bị gãy khúc. 
Những người chuyên vẽ sơ đồ thiết kế dùng một thiết bị cơ học gọi là spline để vẽ các 
đường cong đẹp, có thẩm mỹ: người vẽ xác định tập hợp các điểm (nút) rồi bẻ cong một giải 
plastic hay thanh gỗ linh hoạt (spline) quanh chúng và lấy vết chúng để tạo thành một đường 
cong. Nội suy spline về mặt toán học tương đương với tiến trình này và cho ra cùng một kết quả. 
3.4. PHƯƠNG PHÁP BÌNH PHƯƠNG CỰC TIỂU 
Trong phương pháp nội suy đa thức đã xét, ta đòi hỏi đa thức phải đi qua các điểm quan sát. 
Điều này đôi khi là điều kiện quá chặt trong thực tế. Sau đây ta sẽ xác định tham số của một hàm 
f(x) có dạng đã biết, sao cho các giá trị f(xi) xấp xỉ tốt nhất các giá trị yi theo một tiêu chuẩn gọi 
là bình phương cực tiểu. Trong bài toán ước lượng bình phương cực tiểu ta phải giả thiết là dạng 
hàm f(x) đã biết và ta chỉ cần ước lượng các tham số. Nói chung đối với dạng bài toán này thì ta 
không thể đặt ra yêu cầu là đồ thị hàm số y = f(x) phải đi qua các điểm quan sát, mà chỉ có thể 
đặt ra yêu cầu là đồ thị gần các điểm quan sát nhất trong tập hợp các dạng hàm đã cho. 
3.4.1. Trường hợp hàm nội suy là đa thức hay y phụ thuộc các tham số một cách tuyến tính 
Bây giờ ta xét một trường hợp hay được áp dụng nhất là hàm f(x) có dạng đa thức bậc m, 
tức là: 
f(x) = a0 + a1x1 + . . . am-1xm-1 + amxm 
Vì giá trị f(xi) nói chung khác giá trị yi, giả sử sai số là ei, ta có 
yi = a0 + a1xi1 + . . . am-1 xi m-1 + am xi m + ei 
 i=0,1,2,...,n 
Như vậy ta có: 
ei = yi - ∑ a
=
m
j 0
j xi j 
Và tổng bình phương các sai số bằng 
S = ∑ e
=
n
i 0
i
2 = ∑ (y
=
n
i 0
i - ∑ a
=
m
j 0
j xi j)2 
Để S đạt giá trị nhỏ nhất thì điều kiện sau phải thỏa mãn 
 60 
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Chương 3: Phép nội suy và hồi quy 
∂ S /∂ak =0, với k=0,1,...,m 
Thực hiện phép lấy đạo hàm riêng từng thành phần của tổng theo ak ta có 
∂ (yi - ∑ a
=
m
j 0
j xi j)2 /∂ak = 2(yi - ∑ a
=
m
j 0
j xi j)(- xik) = 2(-yi xik + ∑ a
=
m
j 0
j xi j+k) 
Như vậy 
∂ S /∂ak = 2∑ (-y
=
n
i 0
i xik + ∑ a
=
m
j 0
j xi j+k) = 0, k=0,1,2,...,m 
Từ đây 
∑
=
m
j 0
aj ∑ x
=
n
i 0
i
 j+k = ∑ y
=
n
i 0
i xik , k = 0,1,2,...,m 
Với k = 0,1,2,..,m 
(n+1)a0 + a1∑ x
=
n
i 0
i + a2∑ x
=
n
i 0
i
2 + . . . + am∑ x
=
n
i 0
i
m = ∑ y
=
n
i 0
i
a0∑ x
=
n
i 0
i + a1∑ x
=
n
i 0
i
2 + a2∑ x
=
n
i 0
i
3 + . . . + am∑
=
n
i 0
xim+1 = ∑ y
=
n
i 0
i xi
a0∑ x
=
n
i 0
i
2 + a1∑ x
=
n
i 0
i
3 + a2∑ x
=
n
i 0
i
4 + . . . + am∑ x
=
n
i 0
i
m+2 = ∑ y
=
n
i 0
i xi2
. . . 
a0∑ x
=
n
i 0
i
m + a1∑ x
=
n
i 0
i
m+1 + a2∑ x
=
n
i 0
i
m+2 + . . . + am∑
=
n
i 0
xim+m = ∑ y
=
n
i 0
i xim 
Đặt 
 ∑
=
n
i 0
xi0 ∑
=
n
i 0
xi ∑
=
n
i 0
xi2 ∑
=
n
i 0
xi3
. . . ∑=
n
i 0
xim
 ∑
=
n
i 0
xi ∑
=
n
i 0
xi2 ∑
=
n
i 0
xi3 ∑
=
n
i 0
xi4
. . . ∑=
n
i 0
xim+1
C = ∑
=
n
i 0
xi2 ∑
=
n
i 0
xi3 ∑
=
n
i 0
xi4 ∑
=
n
i 0
xi5
. . . ∑=
n
i 0
xim+1
 . . . 
 ∑
=
n
i 0
xim ∑
=
n
i 0
xim+1 ∑
=
n
i 0
xim+2 ∑
=
n
i 0
xim+3
. . . ∑=
n
i 0
xi2m
 61
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Chương 3: Phép nội suy và hồi quy 
d = (∑ y
=
n
i 0
i xi0 , ∑ y
=
n
i 0
i xi1 ,. . ., ∑ y
=
n
i 0
i xim) 
Ta có hệ phương trình 
C a = d 
Tuy nhiên, để tiện cho việc tính toán, ta có nhận xét sau đây: 
Đặt y = (y0, y1,..., yn)T, e = (e0, e1,..., en)T , a = (a0, a1,..., an)T 
 1 x0 x02 ... x0m
 1 x1 x12 ... x1m
F = . . ... . 
 1 xn xn2 ... xnm
Hệ phương trình C a = d có thể viết dưới dạng sau: 
 FT F a = FTy 
Ta có thể áp dụng phương pháp Gauss-Jordan để giải hệ phương trình này. 
Ta có thể thấy rằng nếu m ≤ n thì định thức của ma trận C khác không và vì vậy đa thức nội 
suy bình phương cực tiểu được xác định duy nhất. 
Thật vậy, viết chi tiết ma trận C ta có 
C = 
⎥⎥
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢⎢
⎢
⎣
⎡
+++++++++
+++++++++
+++++++++
+++
+++
m
n
mmm
n
mmm
n
mm
m
n
mm
nn
m
n
mm
n
xxxxxxxxx
xxxxxxxxx
xxxxxx
22
1
2
0
11
1
1
010
11
1
1
0
22
1
2
010
1010
............
......
............
.........1...11
Bằng cách tách ra các cột ta thu được (n+1)m+1 ma trận con C1, C2,... , mỗi cột của ma trận 
con chỉ phụ thuộc các số 1, x0, x1,..., xn. Sau khi đã tách được như vậy, bằng cách đặt thừa số 
chung ra ngoài ta lại thu được các ma trận mà mỗi ma trận có m+1 cột, các cột này được lấy từ tổ 
hợp của (n+1) các cột có dạng 
⎟⎟
⎟⎟
⎟⎟
⎠
⎞
⎜⎜
⎜⎜
⎜⎜
⎝
⎛
mx
x
x
0
2
0
0
.
1
, , ..., 
⎟⎟
⎟⎟
⎟⎟
⎠
⎞
⎜⎜
⎜⎜
⎜⎜
⎝
⎛
mx
x
x
1
2
1
1
.
1
⎟⎟
⎟⎟
⎟⎟
⎠
⎞
⎜⎜
⎜⎜
⎜⎜
⎝
⎛
m
n
n
n
x
x
x
.
1
2
Dễ thấy rằng: 
- Nếu m>n thì các ma trận con luôn có 2 cột nào đó trùng nhau nên định thức bằng 0 và 
do đó det C = Σ det Ci = 0. 
- Nếu n ≥ m: ma trận C được tách thành hai loại ma trận: 
 62 
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Chương 3: Phép nội suy và hồi quy 
Loại 1: Các ma trận có 2 cột nào đó cùng phụ thuộc một thành phần trong x0,x1,...,xn 
nên có định thức bằng 0. 
Loại 2: Các ma trận có m+1 cột khác nhau lấy từ (n+1) cột có dạng trên. Các ma trận 
này có dạng Vandermond nên khác không. Do đó det C ≠ 0. 
3.4.2. Trường hợp y phụ thuộc các tham số một cách phi tuyến 
a. y = aebx, a>0 (3.19) 
 Lấy logarit hai vế, ta có 
 lny = lna + bx 
 Đặt Y = lny, A = lna, B = B, X = x (3.18) trở thành 
 Y = A + BX (3.20) 
 Như vậy bằng cách lấy logarit hai vế, ta đã đưa quan hệ phi tuyến (3.19) thành dạng tuyến 
tính (3.20). Ta tính được A và B, từ đây tính được a, b. 
b. y = axb, a>0 (3.21) 
 Lấy logarit hai vế, ta có 
 lny = lna + blnx 
 Đặt Y = lny, A = lnA, B = b, X = lnx (3.21) trở thành 
 Y = A + BX (3.22) 
 Từ đây ta tính được A và B, và suy ra a, b. 
Chương trình cài đặt các đa thức nội suy 
 Sau đây là đoạn chương trình chính thể hiện (mô tả) thuật toán hồi qui bằng bình phương 
cực tiểu 
/*Hoi quy dung da thuc uoc luong theo phuong phap binh phuong cuc tieu*/ 
/*Cho truoc bac m, xac dinh cac he so da thuc thuc nghiem , tra ve tong binh phuong cac sai so*/ 
double regresspoli(double *a,int m) 
 {int i,j,k; 
 kmatran aa; 
 double **f,**ft; 
 f = new double* [nqs+1]; 
 for(I=0;i<=nqs;i++) f[i] = new double [m+1]; 
 ft = new double* [m+1]; 
 for(I=0;i<=m;i++) ft[i] = new double [nqs+1]; 
 /*Tinh ma tran 
 f la ma tran co kieu nhu Vandermon nhung co n hang m cot, 
 ft la ma tran chuyen vi cua f. Nhu vay ft x f se la ma tran aa cua 
 he phuong trinh tuyen tinh 
 */ 
 63
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Chương 3: Phép nội suy và hồi quy 
 for(i=0;i<=nqs;i++) 
 {f[i][0]=1; 
 for(j=1;j<=m;j++) 
 f[i][j]=f[i][j-1]*xqs[i]; 
 } 
 //Tinh ma tran chuyen vi 
 for(i=0;i<=m;i++) 
 for(j=0;j<=nqs;j++) ft[i][j]=f[j][i]; 
 /*Tinh ma tran vuong aa cap m, chi can tinh tinh nua tren roi gan cho 
 nua duoi, vi ma tran aa la doi xung 
 */ 
 for(i=0;i<=m;i++) 
 for(j=i;j<=m;j++) 
 {aa[i][j]=ft[i][0]*f[0][j]; 
 for(k=1;k<=nqs;k++) aa[i][j]+=ft[i][k]*f[k][j]; 
 } 
 //Gan nua duoi 
 for(i=0;i<=m;i++) 
 for(j=0;j<i;j++) aa[i][j]=aa[j][i]; 
 //Tinh ve phai cua he phuong trinh 
 for(i=0;i<=m;i++) 
 {aa[i][m+1]=ft[i][0]*yqs[0]; 
 for(k=1;k<=nqs;k++) aa[i][m+1]+=ft[i][k]*yqs[k]; 
 } 
 for(i=0;i<=nqs;i++) delete [] f[i]; 
 for(i=0;i<=m;i++) delete [] ft[i]; 
 gjordan(aa,a,m); 
 //Tinh tong binh phuong cac sai so 
 double ss,fa,xx; 
 ss=0; 
 for(i=0;i<=nqs;i++) 
 {fa=1;xx=1; 
 for(j=1;j<=m;j++) 
 {xx=xx*xqs[i];fa+=a[j]*xx;} 
 ss+=(yqs[i]-fa)*(yqs[i]-fa); 
 } 
 return ss; 
 } 
 64 
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Chương 3: Phép nội suy và hồi quy 
3.5. BÀI TẬP 
Bài 1. Cho hàm số y = 2x với các giá trị tại xi = 3,50 +0,05i, i=0,1,2,3,4 tuần tự là 33,115;34,813; 
36,598; 38,475; 40,477. Hãy lập đa thức nội suy Newton tiến xuất phát từ nút 3,50. 
Bài 2. Tích phân xác suất 
 φ(x) = (2/sqrt(π))∫ exp(-tx
0
2) dt 
 Không tính được bằng nguyên hàm. Người ta tính gần đúng nó tại xi =1 +0,1i, i=0,1,2,..,10 
được các giá trị xấp xỉ tuần tự là 0,8427; 0,8802; 0,9103; 0,9523; 0,9661; 0,9763; 0,9838; 
0,9891; 0,9928; 0,9953. Hãy tính φ(1,43) 
Bài 3. Cho hàm số y = ex tại x = 0,65 + 0,1i i=0,1,...,5 tuần tự là 1,91554; 2,11700; 2,33965; 
2,58571; 2,85765; 3,15819. Hãy tính ln2. 
Bài 4. Biết rằng đại lượng y là một tam thức bậc hai của x và tại x =0,78; 1,56; 2,34; 3,12; 
3,81 các giá trị của y tuần tự là 2,5; 1,20; 1,12; 4,28. Hãy lập công thức y biểu diễn qua x. 
Bài 5. Hãy đánh giá sai số nhận được khi xấp xỉ hàm số y = sinx bằng đa thức nội suy Lagrange 
bậc 5 L5(x), biết rằng đa thức này trùng với hàm số đã cho tại các giá trị x bằng: 00, 50, 100, 
150, 200, 250. 
Xác định gía trị của sai số khi x = 12030'. 
Bài 6. Cho bảng các giá trị: 
x 
2 4 6 8 10 12 
7
.32 
8
.24 
9
.20 
10
.19 
11
.01 
12
.05 
 Hãy áp dụng phương pháp bình phương cực tiểu xác định các đa thức xấp xỉ có các dạng: 
 y = a + bx, y = a+bx +cx2, y = axb
 và so sánh các sai số để kết luận dạng nào thích hợp nhất. 
Bài 7. Thử lại hoặc viết mới các chương trình cài đặt các thuật toán rồi chạy thử để kiểm tra các 
kết quả trên đây. 
 65
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Chương 3: Phép nội suy và hồi quy 
TÓM TẮT NỘI DUNG CHƯƠNG 3 
Trong chương này chúng ta cần chú ý nhất là các vấn đề sau: 
1. Sai số của đa thức nội suy 
Với M = |f
bxa ≤≤
sup (n+1) (x)| khi đó ta có: |R(x)| ≤ | f(x) - pm(x)| ≤ )!1( +n
M | ωn+1(x) | 
đây là công thức đánh giá sai số của đa thức nội suy. 
2. Phương pháp nội suy Lagrange 
Đa thức pn(x) có dạng như sau, được gọi là đa thức nội suy Lagrange: 
pn(x) = y0L0(x) + y1L1(x) + . . . + ynLn(x) (3.6) 
Với: 
L0(x) = 
)x-(x )...x-(x )x-(x
)x-(x )...x-(x )x-(x
n02010
n21 
L1(x) = )x-(x )...x-(x )x-(x
)x-(x )...x-(x )x-(x
n12101
n20 
. . . 
Li(x) = )x-(x )...x-)(xx-(x )...x-(x )x-(x
)x-(x )...x-)(xx-(x )...x-(x )x-(x
ni1ii1-ii2i1i
n1i1-i21
+
+ 
. . . 
Ln(x) = )x-(x )...x-(x )x-(x
)x-(x )...x-(x )x-(x
1-nn2n0n
1-n10 
3. Phương pháp sai phân Newton 
 Phương pháp sai phân tiến Newton với khoảng chia đều 
Theo phương pháp Newton tiến thì đa thức nội suy có bậc m là pm(x) với khoảng chia 
đều như sau: 
pm(x) = y0 + h
y0Δ (x-x0) + 20
2
2h
yΔ
(x-x0)(x-x1) + . . .+ i
i
hi
y
!
0Δ (x-x0) (x-x1)... (x-xi-1) + .. . + 
m
m
hm
y
!
0Δ (x-x0) (x-x1)... (x-xm-1) 
Hoặc có thể biểu diễn công thức trên dưới một dạng khác bằng phép biến đổi t = 
h
xx 0− 
-> x=x0 + th: 
 66 
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Chương 3: Phép nội suy và hồi quy 
pm(x) = pm(x0 + th) = y0 + (Δy0)t + 
2
0
2 yΔ
t(t-1) + . . . 
+ 
!
0
i
yiΔ
t(t-1)... (t-i+1) + . . . + 
!
0
m
ymΔ
t (t-1)... (t-m+1) 
 Phương pháp sai phân tiến Newton với khoảng chia không đều 
 Theo phương pháp Newton tiến thì đa thức nội suy có bậc m là pm(x) với khoảng chia 
không đều như sau: 
pm(x) = a0 + a1(x-x0) + a2(x-x0)(x-x1) + . . . + ai(x-x0)(x-x1)... (x-xi-1)+ . . . 
+ am(x-x0) (x-x1)... (x-xm-1) (3.16) 
 Trong đó: 
a0 = y0
a1 = 
01
01
xx
yy
−
−
 = y[x1,x0]. 
a2 = )(
],[],[
02
0112
xx
xxyxxy
−
−
 = y[x2, x1, x0] 
............................ 
am = y[xm,..., x1, x0] 
 67
CuuDuongThanCong.com https://fb.com/tailieudientucntt