Giáo trình Fortran - Nguyễn Trung Hiếu

ệ (1): Muốn vậy, nhân phương trình (2) tuần tự với 413121 , , aaa và tuần tự lấy các phương 
trình số 2, 3, 4 trừ đi các tích tương ứng vừa nhận được, ta có ba phương trình: 
ù
ỵ
ù
ý
ỹ
=++
=++
=++
)1(
454
)1(
443
)1(
432
)1(
42
)1(
354
)1(
343
)1(
332
)1(
32
)1(
254
)1(
243
)1(
232
)1(
22
axaxaxa
axaxaxa
axaxaxa
 (3) 
trong đó 
)5,4,3,2;4,3,2(11
)1( ==-= jibaaa jiijij (4) 
Bây giờ chia phương trình thứ nhất của hệ (3) cho phần tử chính 
)1(
22a ta có: 
)1(
254
)1(
243
)1(
232 bxbxbx =++ , (5) 
trong đó 
)543()1(
22
)1(
2)1(
2 , , ja
a
b jj == . 
Bằng cách tương tự như khi loại 1x , bây giờ ta loại 2x khỏi các phương trình thứ ba và thứ tư, ta có: 
ùỵ
ù
ý
ỹ
=+
=+
)2(
454
)2(
443
)2(
43
)2(
354
)2(
343
)2(
33
axaxa
axaxa 
. (6) 
trong đó 
)5,4,3;4,3()1(2
)1(
2
)1()2( ==-= jibaaa jiijij . (7) 
Chia phương trình thứ nhất của hệ (6) cho phần tử chính )2(33a , ta có: 
)2(
354
)2(
343 bxbx =+ , (8) 
293 294 
trong đó 
)5,4()2(
33
)2(
3)2(
3 == ja
a
b jj . 
Sau đó nhờ (8) ta loại 3x khỏi phương trình thứ hai của hệ (6), nhận được: 
)3(
454
)3(
44 axa = 
trong đó 
)5 ,4( )2(3
)2(
43
)2(
4
)3(
4 =-= jbaaa jjj (9) 
Như vậy ta đã đưa hệ (1) về hệ tương đương có ma trận các hệ số là ma trận tam giác 
ù
ù
ỵ
ù
ù
ý
ỹ
=
=+
=++
=+++
)3(
454
)3(
44
)2(
354
)2(
343
)1(
254
)1(
243
)1(
232
154143132121
axa
bxbx
bxbxbx
bxbxbxbx
 (10) 
Từ (10) xác định các ẩn 
ù
ù
ỵ
ù
ù
ý
ỹ
---=
--=
-=
=
 bxbxbxbx
 bxbxbx
 bxbx
 aax
)()()(
)()(
)()(
122133144151
1
233
1
244
1
252
2
344
2
353
3
44
3
454
 (11) 
Vậy thủ tục giải hệ phương trình đại số tuyến tính bậc nhất quy về hai quá trình: 
a) Quá trình thuận: đưa hệ (1) về dạng tam giác (10); 
b) Quá trình nghịch: tìm ẩn theo các công thức (11). 
Nếu phần tử chính của hệ bằng không thì chỉ cần thay đổi chỗ của các phương trình trong hệ tương ứng để làm cho phần tử chính khác không. 
Số phép tính số học N cần thực hiện trong phương pháp Gauss bằng 
295 296 
)1(
3
)2()1(2
-+
++
= nnnnnN . 
Vậy số phép tính số học xấp xỉ tỷ lệ với luỹ thừa bậc ba của số ẩn. 
2. Phương pháp căn bậc giải hệ phương trình đại số tuyến tính trong trường hợp ma trận A là ma trận đối xứng 
Phương pháp này thuận lợi trong trường hợp hệ phương trình 
A x = b (12) 
có ma trận A là ma trận đối xứng, điều thường gặp trong các bài toán kỹ thuật. 
Theo phương pháp này ma trận A được biểu diễn thành tích của hai ma trận tam giác chuyển vị 
TTA Â= (13) 
trong đó 
ữữ
ữ
ữ
ữ
ứ
ử
ỗỗ
ỗ
ỗ
ỗ
ố
ổ
=Â
ữữ
ữ
ữ
ữ
ứ
ử
ỗỗ
ỗ
ỗ
ỗ
ố
ổ
=
nnnnnn
n
n
ttt
tt
t
T,
t
tt
ttt
T
...
............
0...
0...
...00
............
...0
...
21
2212
11
222
11211
Nhân hai ma trận T Â và T và cho tích bằng ma trận A , ta suy ra cá công thức tính các phần tử ijt : 
jikhit
ji
t
tta
t
nitat
j
t
a
tat
ij
ii
i
k
kjkiij
ij
i
k
kiiiii
j
j
>=
<
-
=
Ê<-=
>==
ồ
ồ
-
=
-
=
0
)(
)1(
)1(,
1
1
1
1
2
11
1
11111
 (14) 
Như vậy ta đã thay hệ (12) bằng hai hệ tương đương 
T' y = b, T x = y 15) 
hay 
297 298 
ù
ù
ỵ
ù
ù
ý
ỹ
=+++
=+
=
 nnnnnn bytytyt
bytyt
byt
....
.......................
2211
2222112
1111
 (16) 
ù
ù
ỵ
ù
ù
ý
ỹ
=
=++
=+++
nnnn
nn
nn
yxt
yxtxt
yxtxtxt
.........................
...
....
22222
11212111
 (17) 
Từ đó suy ra các công thức tính: 
)1(,
1
1
11
1
1 >
-
==
ồ
-
= i
t
ytb
y
t
by
ii
i
k
kkii
i (18) )(,
1 ni
t
xty
x
t
yx
ii
n
ik
kiki
i
nn
n
n <
-
==
ồ
+= (19) 
Vậy quá trình thuận gồm tính các phần tử của ma trận T theo các công thức (14). Quá trình nghịch là tính các ma trận cột y và x theo các công thức (18), (19). 
Phụ lục 3: Phương pháp bình phương nhỏ nhất trong phân tích hồi quy 
1. Mô hình tuyến tính 
299 300 
Mô hình hồi quy tuyến tính có dạng: 
baxxfy +== )( . 
Theo phương pháp bình phương nhỏ nhất, các hệ số hồi quy a và b trong phương trình trên được tìm sao cho tổng bình phương sai số bằng 
ồ
=
--=
n
k
kk baxyE
1
2)( 
cực tiểu. Lần lượt lấy đạo hàm biểu thức này theo a , b và cho bằng không, ta được hệ phương trình sau đây để xác định a và b : 
ồồồ
===
=+
n
k
kk
n
k
k
n
k
k yxxbxa
111
2
, 
ồồ
==
=+
n
k
k
n
k
k ynbxa
11
 . 
Vậy các hệ số hồi quy được tính theo các công thức sau: 
ồồ
ồồồ
==
===
-
-
= n
k
k
n
k
k
n
k
kk
n
k
k
n
k
k
xnx
yxnyx
a
1
22
1
111
 )(
 (20) 
ồồ
ồồồồ
==
====
-
-
= n
k
k
n
k
k
n
k
k
n
k
k
n
k
kk
n
k
k
xnx
yxyxx
b
1
22
1
11
2
11
 )(
, (21) 
hay hệ số b còn có thể tính theo công thức: 
n
xay
b
n
k
k
n
k
k ồồ
==
-
= 11
. (22) 
2. Mô hình đa thức 
Phương pháp bình phương nhỏ nhất cũng có thể áp dụng để tính các hệ số hồi quy đa thức dạng 
301 302 
m
n xaxaxaaxf ++++= ...)(
2
210 . 
thí dụ đối với mô hình bậc hai 
2
210)( xaxaaxf ++= . 
Lấy đạo hàm tổng sai số theo các hệ số và cho bằng không ta có hệ sau đây để xác định các hệ số hồi quy bậc hai: 
ù
ù
ù
ù
ợ
ùù
ù
ù
ớ
ỡ
=++
=++
=++
ồồ ồồ
ồồ ồồ
ồồồ
== ==
== ==
===
n
k
kk
n
k
n
k
kk
n
k
k
n
k
kk
n
k
n
k
kk
n
k
k
n
k
k
n
k
k
n
k
k
yxxaxaxa
yxxaxaxa
ynaxaxa
1
2
1 1
2
0
3
1
1
4
2
11 1
0
2
1
1
3
2
11
01
1
2
2
 (23) 
Về nguyên tắc ta có thể sử dụng phương pháp này để tìm phương trình đa thức bậc bất kỳ. Tuy nhiên trong thực tế phương pháp trở thành không ổn định khi bậc đa 
thức lớn hơn vì các sai số làm tròn số trong máy tính. 
3. Mô hình phi tuyến 
Phương pháp bình phương nhỏ nhất có thể áp dụng cho hàm bất kỳ, nhưng hệ các phương trình để tìm các hệ số có thể phi tuyến, và do đó không thể giải được bằng 
cách sử dụng các phương trình tuyến tính. Tuy nhiên, trong một số trường hợp, một hàm phi tuyến có thể chuyển thành một hàm tuyến tính. Thí dụ về một hàm có thể 
tuyến tính hoá là 
axbxf )( = (24) 
Nếu lấy loga hai vế của phương trình này, ta có 
bxaxf lnln)(ln += . (25) 
Nếu ký hiệu 
)(ln)( xfxg = (26) 
bb ln~ = (27) 
xx ln~ = (28) 
303 304 
yy ln~ = (29) 
ta có 
bxaxg ~~)( += (30) 
Với phương trình (30) các hệ số hồi quy a và b~ tính theo các công thức 
ồồ
ồồồ
==
===
-
-
= n
k
k
n
k
k
n
k
kk
n
k
k
n
k
k
xnx
yxnyx
a
1
22
1
111
~ )~(
~~ ~ ~
 (31) 
ồồ
ồồồồ
==
====
-
-
= n
k
k
n
k
k
n
k
k
n
k
k
n
k
kk
n
k
k
xnx
yxyxx
b
1
22
1
11
2
11
~ )~(
~ ~ ~~ ~
~
 (32) 
Vậy công việc tính toán gồm: chuyển đổi các giá trị số liệu kx và ky theo các công thức (28), (29), tính các tổng, kết quả thế vào các phương trình (31), (32) để tìm 
a và b~ . Giải phương trình (27) đối với b và đặt vào phương trình (24). 
305 306 
Phụ lục 4: Sơ đồ ứng dụng phương pháp hồi quy nhiều biến 
Giả sử có n quan trắc đối với biến phụ thuộc y và các biến độc lập mxxx ..., , , 21 . Phương trình hồi quy được thiết lập như sau: 
mm xaxaxaay ++++= ...22110 . 
Các hệ số hồi quy ),...,1( miai = được chọn sao cho thoả mãn 
( )ồ
=
=-----=
n
i
mm xaxaxaay
1
2
22110 min...d 
Lần lượt lấy đạo hàm biểu thức trên theo maaaa ,..., , , 210 và cho các đạo hàm bằng không, ta có hệ 1+m phương trình để xác định các hệ số a 
[ ] [ ] [ ] [ ]
[ ] [ ] [ ] [ ] [ ]
[ ] [ ] [ ] [ ] [ ]
[ ] [ ] [ ] [ ] [ ]mmmmmmm
mm
mm
mm
yxaxxaxxaxxax
yxaxxaxxaxxax
yxaxxaxxaxxax
yaxaxaxna
=++++
=++++
=++++
=++++
...
..................
...
...
...
22110
2222212102
1121211101
22110
 (33) 
307 308 
Hệ phương trình này gọi là hệ phương trình chính tắc để xác định các hệ số hồi quy. Dưới dạng ma trận ta viết hệ này như sau: 
[ ] [ ] [ ]
[ ] [ ] [ ] [ ]
[ ] [ ] [ ] [ ]
[ ] [ ] [ ] [ ] ữ
ữ
ữ
ữ
ữ
ữ
ứ
ử
ỗ
ỗ
ỗ
ỗ
ỗ
ỗ
ố
ổ
=
ữ
ữ
ữ
ữ
ữ
ữ
ứ
ử
ỗ
ỗ
ỗ
ỗ
ỗ
ỗ
ố
ổ
ữ
ữ
ữ
ữ
ữ
ữ
ứ
ử
ỗ
ỗ
ỗ
ỗ
ỗ
ỗ
ố
ổ
mmmmmmm
m
m
m
b
b
b
b
a
a
a
a
xxxxxxx
xxxxxxx
xxxxxxx
xxxn
......
.
...
...............
...
...
...
2
1
0
2
1
0
21
222212
112111
21
(34) 
với dấu [ ] ký hiệu phép lấy tổng ồ
n
1
. 
Để tìm các hệ số hồi quy maaaa ,...,,, 210 ta phải giải hệ phương trình chính tắc theo phương pháp loại biến Gauss hoặc phương pháp căn bậc hai đã mô tả trong 
phụ lục 2 vì ma trận hệ số của các phương trình chính tắc là ma trận đối xứng. Dưới đây dẫn hai thủ tục hỗ trợ cho việc lập hệ phương trình đại số tuyến tính chuẩn tắc 
(34) - SUBROUTINE LHPTCT và giải hệ phương trình đó bằng phương pháp loại biến Gauss - SUBROUTINE GAUSS. 
 SUBROUTINE LHPTCT (Y, X, A, N, M) 
 INTEGER N, M, I, J, K 
 REAL Y(10000), X(10000, 50), A(0 : 50, 0 : 51) 
 A (0, 0) = N 
 DO J = 1, M 
 A (0, J) = 0.0 
 DO K = 1, N 
 A (0, J) = A (0, J) + X (K, J) 
 END DO 
 END DO 
 A (0, M + 1) = 0.0 
 DO K = 1, N 
 A (0, M + 1) = A (0, M + 1) + Y (K) 
 END DO 
 DO I = 1, M 
 A (I, M + 1) = 0.0 
 DO K = 1, N 
 A (I, M + 1) = A (I, M + 1) + Y (K) * X(K, I) 
 END DO 
309 310 
 END DO 
 DO I = 1, M 
 DO J = I, M 
 A (I, J) = 0.0 
 DO K = 1, N 
 A (I, J) = A (I, J) + X (K, I) * X (K, J) 
 END DO 
 ENDDO 
 ENDDO 
 DO I = 1, M 
 DO J = 0, I - 1 
 A (I, J) = A (J, I) 
 END DO 
 END DO 
 RETURN 
 END 
 SUBROUTINE GAUSS (M, A, X) 
 INTEGER M 
 REAL A (0 : 50, 0 : 51), X (0 : 50) 
 REAL MAX 
 DO I = 0, M - 1 
 K = I 
 MAX = ABS (A (K, K)) 
 DO J = I + 1, M 
 IF (MAX .LT. ABS (A (J, I))) THEN 
 MAX = ABS (A (J, I)) 
 K = J 
 END IF 
 END DO 
 IF (K .NE. I) THEN 
 DO J = I, M + 1 
 MAX = A (I, J) 
 A (I, J) =A (K, J) 
 A (K, J) = MAX 
 END DO 
 END IF 
311 312 
 DO J = I + 1, M + 1 
 A (I, J) = A (I, J) / A (I, I) 
 END DO 
 DO J = I + 1, M 
 DO K = I + 1, M + 1 
 A (J, K) = A (J, K) - A (J, I) * A (I, K) 
 END DO 
 END DO 
 END DO 
 X (M) = A (M, M + 1) / A (M, M) 
 DO I = M - 1, 0, -1 
 X (I) = A (I, M + 1) 
 DO J = I + 1, M 
 X (I) = X (I) - A (I, J) * X (J) 
 END DO 
 END DO 
 RETURN 
 END