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

 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