Tính toán thiết kế tàu và tự động hóa
Tài liệu giúp các bạn tìm hiểu nội dung các chương sau: chương 1 phương pháp tính và tự động hóa tính toán thiết kế tàu, chương 2 tính toán nối và tính ổn định tàu, chương 3 sức cản vỏ tàu, chương 4 thiết kế chân vịt tàu thủy, chương 5 thiết kế tối ưu tàu thủy.
TRẦN CÔNG NGHỊ
TỰ ĐỘNG HÓA
TÍNH TOÁN
THIẾT KẾ TÀU
ĐẠI HỌC GIAO THÔNG VẬN TẢI TP HỒ CHÍ MINH
7 - 2009
Trang để trống
Trần công nghị
TỰ ĐỘNG HÓA
TÍNH TOÁN
THIẾT KẾ TÀU
THÀNH PHỐ HỒ CHÍ MINH 2009
ĐẠI HỌC GIAO THÔNG VẬN TẢI TP HỒ CHÍ MINH
Mục lục
Mở đầu 5
Chương 1: Phương pháp tính và tự động hóa ính toán thiết kế tàu 6
1.1 Nội suy Lagrange 6
1.2 Tích phân một lớp 8
1.3 Đa thức Legendre 11
1.4 Đa thức Tchebyshev 14
1.5 Tìm nghiệm bằng phương pháp chia đôi đoạn có nghiệm 15
1.6 Phương pháp tổng nhỏ nhất các bình phương 16
1.7 Qui hoạch tuyến tính 21
1.8 Qui hoạch phi tuyến 27
1.8.1 Hàm một biến 28
1.8.2 Hàm nhiều biến 29
1.8.3 Xác định min/max hàm một biến 30
1.8.4 Phương pháp sử dụng gradient 34
1.8.5 Phương pháp tìm trực tiếp (không qua giai đoạn tính gradient). 39
1.8.6 Phương pháp dùng hàm phạt penalty 43
Chương 2: Tính nổi và tính ổn định tàu 51
2.1 Tính nổi tàu thủy 51
2.1.1 Kích thước chính và các hệ số thân tàu 51
2.1.2 Tỷ lệ Bonjean 55
2.1.3 Tính thể tích phân chìm và cá đại lượng liên quan thể tích 55
2.1.4 Tính các đường thủy tĩnh trên máy cá nhân 58
2.1.5 Biểu đồ Firsov 60
2.2 Ổn định tàu. 61
2.2.1 Ổn định ngang ban đầu. 61
2.2.2 Ổn định tại góc nghiêng lớn. 62
2.2.3 Đồ thị ổn định. 65
2.3 Thuật toán xác lập họ đường Cross Curves (pantokaren) 65
2.4 Giới thiệu chương trình tính tính nổi tàu thủy 70
Chương 3: Sức cản vỏ tàu 78
3.1 Sức cản vỏ tàu 77
3.2 Công suất hữu hiệu 81
3.3 Các phương pháp kinh nghiệm tính sức cản vỏ tàu 81
Chương 4: Thiết kế chân vịt tàu thủy 95
4.1 Đặc tính hình học chân vịt 95
4.2 Vẽ chân vịt 97
4.3 Đặc tính thủy động lực 98
4.4 Đồ thị thiết kế chân vịt 102
4.5 Tính hệ số dòng theo, hệ số lực hút 106
4.6 Xâm thực chân vịt 109
4.7 Độ bền cánh chân vịt 114
4.8 Thiết kế chân vịt bước cố định 118
4.9 Lập chương trình thiết kế chân vịt tàu 126
4.10 Vẽ chân vịt trên máy PC 135
3
Chương 5: Thiết kế tối ưu tàu thủy 148
5.1 Đánh giá các chỉ tiêu kinh tế – kỹ thuật của tàu 148
5.2 Sơ đồ tính hiệu quả kinh tế 150
5.3 Tự động thiết kế tàu vận tải 151
Tài liệu tham khảo 172
4
Mở đầu
“Tự động hóa tính toán, thiết kế tàu” trình bày cách tính toán, thuật toán phục vụ việc lập
chương trình tính tính năng tàu thủy, tính di chuyển, thiết bị đẩy tàu và tự động hóa vẽ tàu. Sau mười
năm sử dụng sách cho chuyên đề này những người viết chỉnh, sửa, viết lại phù hợp thực tế. Sửa
chữa và bổ sung lần này nhằm làm cho tài liệu sát đề cương giảng dạy và học tập tại trường Đại học
Giao thông Vận tải Tp Hồ Chí Minh.
Hy vọng rằng sách có ích cho những người đang theo học đóng tàu và công trình nổi cũng
như các đồng nghiệp đang làm việc trong cùng lĩnh vực.
Thành phố Hồ Chí Minh tháng 6 năm 2009.
Người viết
“Mở đầu” lần in thứ nhất
“Tự động hóa tính toán, thiết kế và đóng tàu” bao gồm hướng dẫn tính toán, chương
trình tính phục vụ những môn học tàu thủy tại trường đại học. Những đề tài trong tài liệu này: Thiết
kế tàu, Tính nổi và tính ổn định, Sức cản vỏ tàu và thiết bị đẩy tàu, Qui hoạch tuyến tính, qui
hoạch phi tuyến và ứng dụng của lý thuyết này vào thiết kế tối ưu tàu thủy, Spline và ứng dụng
trong vẽ đường hình, khai triển vỏ tàu.
Tài liệu được bố trí theo cách tiện lợi cho người đọc. Mở đầu mỗi chương bạn đọc có điều
kiện ôn lại những hiểu biết cần thiết về các phương pháp tính liên quan đến nội dung của chương,
có điều kiện làm quen chương trình tính viết bằng ngôn ngữ C áp dụng trong tính toán. Các
chương trình nhỏ này còn được dùng cho những vấn đề liên quan với ngành tàu. Nội dung mỗi
chương chỉ gồm những kiến thức đã được truyền đạt trong trường đại học chuyên ngành. Trên cơ sở
những vấn đề đang được trình bày bạn đọc tìm hiểu thêm giải thuật xử lý những bài toán cụ thể đang
đặt ra và cách hoàn thiện một chương trình máy tính dựa vào giải thuật vừa có.
Tài liệu có thể giúp ích cho sinh viên khoa đóng tàu, kỹ sư làm việc trong lĩnh vực đóng
sửa tàu, thiết kế, nghiên cứu tàu cùng đông đảo bạn đọc quan tâm đến tàu thủy khi tính toán tính
năng, như tính nổi, ổn định, tính sức cản, chọn máy phù hợp, thiết kế mới, lập phương án đóng mới,
lập phương án sửa chữa tàu vv...
Trong quá trình biên soạn tài liệu những người làm công tác chuẩn bị nhận được sự giúp đỡ
chân tình và thiết thực của ban giám hiệu phân hiệu Đại học Hàng hải, phân khoa đóng tàu, bạn
cùng nghề và những bạn bè xa, gần. Những đóng góp quí giá về nội dung, về biên soạn và hiệu chỉnh
tài liệu, hiệu chỉnh các bản in thử vv… đã làm cho tài liệu có nội dung phù hợp hơn, tránh được
nhiều sai sót. Xin chân thành cám ơn về sự đóng góp quí giá trên. Người viết căn cứ sự giúp đỡ, chỉ
dẫn trên đã cố gắng hoàn chỉnh tài liệu kịp ra mắt bạn đọc, tuy nhiên vì khả năng có hạn chắc rằng
trong tài liệu vẫn còn những sai sót khó tránh. Rất mong bạn đọc gần, xa góp thêm ý kiến nhằm làm
cho tài liệu ngày càng hoàn thiện. Thư, bài góp ý, xây dưng xin gửi về phân hiệu Đại học Hàng Hải,
thành phố Hồ Chí Minh.
Thành phố Hồ Chí Minh tháng 12 năm 2000.
5
Chương 1
PHƯƠNG PHÁP TÍNH VÀ TỰ ĐỘNG HÓA TÍNH TOÁN,
THIẾT KẾ TÀU
Trong chương này giới thiệu những phương pháp tính thông dụng dùng xử lý những vấn đề
thường gặp trong tính tính nổi tàu thủy, tính ổn định tàu, thiết kế máy đẩy tàu, thiết kế tối ưu tàu
thủy.
1.1 NỘI SUY LAGRANGE
Đa thức nội suy Lagrange được viết dưới dạng 1 :
f(x) = pn(x) + Rn(x), (1.1)
hoặc dạng đầy đủ:
n
⎡n ⎤ f
( n +1)
(ξ )
f ( x) = ∑ Li ( x) f ( xi ) + ⎢ Π ( x − xi )⎥ ,a < ξ < b (1.2)
i =0 ⎣i = 0 ⎦ (n + 1)!
n ⎛ x− x ⎞
trong đó Li ( x) = ∏ ⎜ ⎟
j
(1.3)
⎜x −x
j =0 ⎝ i
⎟
j ≠i
j ⎠
n
Đa thức p n ( x) = ∑ Li ( x) f ( xi ) mang tên gọi đa thức Lagrange, còn vế sau của phía phải công thức
i =0
gọi hàm sai số.
Đa thức pn(x) mặt khác được hiểu là đa thức bậc n, có dạng:
pn(x) = a0(x – x1) (x – x2)... (x – xn) +
+ a1(x – x0) (x – x2)... (x – xn) +
+ a2(x – x0) (x – x1)... (x – xn) +
...
+ ai(x – x0) (x – x1)... (x – xi-1) (x – xI+1)... (x – xn)
...
an(x – x0) (x – x1)... (x – xn-2)(x – xn-1 ) (1.4)
Các hệ số a0, a1, a2,... tính từ quan hệ:
pn(xi) = f(xi) = yi ; i = 0, 1, 2,... (1.5)
Lần lượt thay x = x0, x = x1,... vào công thức cuối có thể xác định công thức tính các hệ số.
Ví dụ, từ pn(x0) = y0 = a0(x0 – x1) (x0 – x2)... (x0 – xn) sẽ nhận được:
f ( x0 )
a0 =
( x0 − x1 )( x0 − x 2 )...( x0 − x n )
tương tự vậy có thể viết:
1
R.W. Hamming, “Numerical Methods for Scientists and Engineers”, McGraw-Hill, N.Y, 1962,
F.B. Hildebrand, “Introduction to Numerical Analysis”, McGraw-Hill, N.Y., 1956.
6
f ( x1 )
a1 =
( x1 − x0 )( x1 − x 2 )...( x1 − x n )
…
f ( xn )
an =
( x n − x0 )( x n − x1 )...( x n − x n −1 )
Hệ số thứ i mang dạng chung:
f ( xi )
ai = (1.6)
( xi − x0 )( xi − x1 )...( xi − xi −1 )( xi − xi +1 )...( xi − x n )
Thay các biểu thức vừa xác định vào vị trí a0, a1,..., an sẽ nhận được công thức nội suy hay còn gọi
đa thức Lagrange:
( x − x1 )( x − x 2 )...( x − x n ) ( x − x0 )( x − x 2 )...( x − x n )
p n ( x) = y0 + y1 +
( x 0 − x1 )( x0 − x 2 )...( x0 − x n ) ( x1 − x0 )( x1 − x 2 )...( x1 − x n )
(1.7)
( x − x0 )( x − x1 )...( x − x n −1 )
+ ... + yn
( x n − x0 )( x n − x1 )...( x n − x n −1 )
n n ⎛ x− x ⎞
p n ( x ) = ∑ Li ( x) f ( xi ) , với Li ( x) = ∏ ⎜ ⎟.
j
hoặc dưới dạng gọn hơn như đã trình bày
⎜x −x
j =0 ⎝ i
⎟
i =0
j ≠i
j ⎠
Những trường hợp riêng lẻ của hàm nội suy Lagrange như sau.
với n =1:
( x − x1 ) ( x − x0 )
p1 ( x) = y0 + y1 (1.8)
( x0 − x1 ) ( x1 − x0 )
với n = 2:
( x − x1 )( x − x 2 ) ( x − x0 )( x − x 2 ) ( x − x0 )( x − x1 )
p 2 ( x) = y0 + y1 + y2 (1.9)
( x0 − x1 )( x0 − x 2 ) ( x1 − x0 )( x1 − x 2 ) ( x 2 − x0 )( x 2 − x1 )
Hàm p1(x) là đoạn thẳng qua hai điểm (x0, y0) (x1, y1 ), có tên gọi công thức nội suy tuyến tính. Hàm
thứ hai là đường parabol bậc hai qua ba điểm cho trước, gọi là nội suy bậc hai.
Chương trình hóa phương pháp nội suy Lagrange được thể hiện bằng thuật toán Neville và minh họa
tại hàm bằng ngôn ngữ C sau.
#include (math.h>
void Lagrange(xa, ya, n, x, y, dy)
float xn[], ya[], x,*y, *dy;
int n;
{
int i, m, ns=1;
float den, dif, dift, h0, hp, w;
float *c, *d, *vector();
dif = fabs( x-xa[q]);
c = vector(1,n);
d = vector(1,n);
for (i=1; i dif = dift;
}
c[i] = ya[i];
d[i] = ya[i];
}
*y = ya[ns--];
for ( m=1; m n −1 xk +1
b k = n −1 1 n −1
f ( x k ) + f ( x k +1 )
∫ f ( x)dx ≅ ∑ ∫ f ( x k )dx = ∑ h ∫ [ f ( x k ) + pΔf ( x)]du = h∑ =
a k = 0 xk k =0 0 k =0 2
= h[ f ( x 0 ) / 2 + f ( x1 ) + ... + f ( x n −1 ) + f ( x 2 ) / 2]
(3) Công thức Simpson được tính theo luật trên đây khi giữ lại ba thành phần đầu tiên của chuỗi.
b 2 n−2 xk + 2 2
2n−2
p( p − 1) 2
∫ f ( x)dx ≅ ∑ ∫ f ( x )dx = ∑ h∫
a k =0, 2,... xk
k
k =0, 2,... 0
[ f ( xk ) + pΔf ( xk ) +
2!
Δ f ( xk )]du =
2 n−2
h h
= ∑,...f ( xk ) + 4 f (xk +1 ) + f ( xk +2 )] = 3 [ f ( x0 ) + 4 f ( x1 ) + 2 f ( x2 ) + ... + 4 f ( x2n−1 ) + f (x2n )]
3 k =0 , 2
[
Hàm bằng ngôn ngữ C, thực hiện tích phân theo phương pháp hình thang được trình bày tiếp dưới
đây.
#define FUNC(x) ((*func) (x) )
float trapezd( func, a, b, n)
float a, b;
float (*func) ();
int n;
{
float x, tnm, sum, del;
static float s;
static int it;
int j;
if (n == 1 ) {
it=1;
return (s=0.5* (b-a) *(FUNC(a) +FUNC(b) ));
}
else
{
tnm = it;
del = (b-a)/tnm;
x = a + 0.5*del;
for (sum=0.0, j=1; j ⎡ 1 1 x3 − x1 ⎤
a3 = ( x3 − x1 ) ⎢ − ⎥;
⎣ 2 6 x3 − x 2 ⎦
1 ( x3 − x1 ) 3
a2 = ; (1.13)
6 ( x 2 − x1 )( x3 − x 2 )
a1 = x3 − x1 − a 2 − a3
Hàm viết bằng ngôn ngữ C xử lý tích phân giới hạn biến thiên được thể hiện như sau.
x
integrationv = ∫ f ( x)dx
0
void integrationv(int N, float *X, float *f, float *Ans)
{
register i,ii,n2 ;
float d1,d2,d3,a1,a2,a3 ;
n2 = N / 2;
Ans[0]=0.0;
for (i = 1; i + f[k+1]*((X[k]-X[k-1])/(X[k+1]-X[k-1])+2.0)-f[k-1]*
sqr(X[k+1]-X[k])/((X[k]-X[k-1])*(X[k+1]-X[k-1])));
return Sum;
}
1.3 ĐA THỨC LEGENDRE
Đa thức Legendre trực giao trong đoạn [-1, 1], được hiểu như sau:
+1
∫ P ( x) P
−1
n m ( x)dx = 0,..n ≠ m (1.14)
+1
∫ P ( x) P ( x)dx = c(n) ≠ 0
−1
n n
Một số ít đa thức Pn (x) được viết dưới đây:
P0(x) =1,
P1(x) =x,
P2(x) = ½ (3x2 - 1),
P3(x) = ½ (5x3 – 3x),
P(x) = 1/8 (35x4 – 30x2 + 3)
Quan hệ hồi qui của đa thức Legendre có dạng:
2n − 1 n −1
Pn ( x ) = xPn −1 ( x) − Pn − 2 ( x) . (1.15)
n n
Áp dụng tính trực giao của đa thức Pn(x) để xác định các tọa độ xj nhằm giảm thiểu sai số khi tính.
Công thức tính tích phân Gauss-Legendre có dạng:
b b b
∫
a
f ( x)dx = ∫ p n ( x)dx + ∫ Rn ( x)dx
a a
(1.16)
Trong đó hàm f(x) có thể viết: f(x) = pn(x) + Rn(x).
n
⎡n ⎤ f
( n +1)
(ξ )
f ( x) = ∑ Li ( x) f ( xi ) + ⎢ Π ( x − xi )⎥ ,a < ξ < b (1.17)
i =0 ⎣i = 0 ⎦ (n + 1)!
n ⎛ x− x ⎞
trong đó Li ( x) = ∏ ⎜ ⎟
j
⎜x −x
j =0 ⎝ i
⎟
j ≠i
j ⎠
Nếu sử dụng biến z theo cách sau đây:
a+b b−a
z= + x (1.18)
2 2
n ⎛ z− z ⎞
Li ( z ) = ∏ ⎜ ⎟ còn
j
hàm f(x) có thể thay bằng F(z), hàm Lagrangre Li(x) có thể thay bằng
⎜z −z
j =0 ⎝ i
⎟
j ≠i
j ⎠
giới hạn tích phân trở thành –1 và +1, biến ξ nằm trong giới hạn -1 < ξ < +1.
Tích phân (1.19 ) trở thành:
11
+1 +1 +1
n
⎡ n ⎤
∫1 F ( z )dz = ∫ ∑ Li ( z )F ( z i ) dz + ∫ ⎢∏ ( z − z i ) ⎥ q n ( z )dz (1.19)
− −1 i = 0 −1 ⎣ i = 0 ⎦
Nếu bỏ vế sau bên phía phải của công thức có thể viết:
+1 n +1 n
∫ F ( z )dz ≅ ∑ F ( z i ) ∫ Li ( z )dz = ∑ Wi F ( z i )
−1 i =1 −1 i =1
(1.20)
Từ công thức cuối có thể xác định:
+1
+1
∫ ( z − z 0 )...( z − z i −1 )( z − z i +1 )...( z − z n )dx
wi = ∫ L ( z )dz = −1
i (1.21)
−1
( z − z 0 )...( z − z i −1 )( z − z i +1 )...( z − z n )
+1
⎡ n ⎤
Vế thứ hai của phía phải ∫1 ⎢∏ ( z − zi )⎥ q n ( z )dz chính là sai số của phép tích phân số đang được xét.
− ⎣ i =0 ⎦
Trong giai đoạn này cần thiết xác định vị trí của zi trong đa thức Legendre nhằm làm cho vế này
trượt tiêu. Tại đây sử dụng tính trực giao của hàm Legendre để đạt điều momg muốn. Khai triển
n
hai đa thức qn(z) và ∏ ( z − z ) dưới dạng sau đây:
i =0
i
n n +1
∏ ( z − z ) = b P ( z ) + b P ( z ) + ... = ∑ b P ( z
i =0
i 0 0 1 1
i =0
i i ) (1.22)
và
n
q(z) = c0P0(z) +c1P1(z) + … = ∑ c P ( z)
i =0
i i (1.23)
Thay hai đa thức cuối vào biểu thức thuộc vế hai phía phải công thức (1.19) có thể viết:
⎡n n
+1 n ⎤
∫1 ⎢∑∑
− ⎣ i =0 j =0
bi c j Pi ( z ) Pj ( z ) + bn +1 ∑ ci Pi ( z ) Pn +1 ( z )⎥dz
i =0 ⎦
(1.24)
Từ tính trực giao của đa thức Legendre có thể viết:
+1
bi c j ∫ Pi ( z ) Pj ( z )dz = 0; i ≠ j (1.25)
−1
Từ đó, sai số phép tích phân được viết là:
+1 +1 +1
⎡ n ⎤ n n
∏ ( z − z i )⎥ q n ( z )dz = ∫ ∑ bi ci [Pi ( z )] dz = ∑ bi ci ∫ [Pi ( z )] dz
∫1 ⎢ i =0
2 2
(1.26)
− ⎣ ⎦ −1 i = 0 i =0 −1
Từ phương trình cuối có thể xác định các nút tính toán. Các nút tìm bằng cách này gọi là các điểm
zero của đa thức Legendre.
b
Giới hạn a, b trong công thức chung ∫a
f ( z )dz , có quan hệ với giới hạn chuẩn của hàm f(x) như đã
a+b b−a
trình bày z = + x.
2 2
Công thức tính tích phân theo Gauss-Legendre có dạng:
12
b n
∫ f ( x)dx ≅ (b − a)∑ Wi f ( z i ) (1.27)
a i =1
Một số giá trị của nút và hệ số tính cho công thức trên như sau:
Bảng 1.1
n=1 x1= 0,0 w1 = 2,0
n=2 x1,2 = ±0,577350 a1,2 = 1,0
n=3 x3,1 = ±0,774597 a1,3 = 0,555556
x2 = 0,0 a2 = 0,888889
n=4 x3,1 = ±0,861136 a1,4 = 0,347855
x4,2 = 0,339981 a2,3 = 0,652145
Ví dụ tính theo công thức Gauss-Legendre.
2 dx
Tính giá trị tích phân sau ∫
1 x
2 dx
Thay thế biến z = 2x – 3 và hàm f(x) thành F(z) = 2/ (z+3). Tích phân ∫1 x
có giá trị bằng tích
2
+1
phân ∫
−1 z + 3
dz . Kết quả tính, với n = 4 sẽ có dạng:
Bảng 1.2
i zi wi F(zi) wF(zi)
0 0,0 0,5688889
0,333333 0,189629
1 0,5384693 0,4786287
0,282608 0,135264
2 -0,5384693 0,4786287
0,406251 0,1944435
3 0,9061798 0,2369269
0,256004 0,0606544
4 -0,9061798 0,2369269
0,4775959 0,1131553
Cộng: 0,6931471
Hàm viết bằng ngôn ngữ C cho trường hợp n =10 được ghi lại dưới đây.
float qgaus (func, a, b)
/* Gauss-Legendre formula */
float a,b;
float (*func) ();
{
int j;
float xr, xm, dx, s;
static float x[] = {0.0, 0.1488743389, 0.4333953941,
0.679409568, 0.8650633666, 0.9739065285};
static float w[] = { 0.0, 0.2955242247, 0.269266719,
0.2190863625, 0.149451349, 0.0666713443};
xm=0.5* (b+a);
xr=0.5*(b-a);
s=0.0;
for (j=1; j }
return s *= xr;
}
1.4 ĐA THỨC TCHEBYSHEV
Đa thức Tchebyshev, ký hiệu từ ký tự đầu tên nhà bác học Tn(x), trực giao trong [-1, 1] cùng hàm
trọng lượng w(x) = 1/ √ (1 – x2).
+1
1
∫
−1 1− x2
Tn ( x)Tm ( x)dx = 0,..n ≠ m (1.28)
+1
1
−1
∫1− x2
Tn ( x)Tn ( x)dx = c(n) ≠ 0 (1.29)
Những đa thức đầu có dạng sau.
T0(x) =1,
T1(x) =x,
T2(x) = 2x2 - 1,
T3(x) = 4x3 – 3x.
Quan hệ hồi qui của đa thức Tchebyshev có dạng:
T(x) = 2xTn-1(x) – Tn-2(x) (1.30)
Công thức Gauss – Tchebyshev.
Tích phân hàm f(x), giới hạn a, b theo công thức Gauss – Tchebyshev sẽ là:
b +
(b − a ) ⎛ a + b b − a ⎞
∫
a
f ( x)dx = ∫ f⎜
2 −1 ⎝ 2
+
2 ⎠
z ⎟dz (1.31)
Hàm trọng lượng được trình bày trên, cho phép viết tiếp:
b +
(b − a) 1
∫
a
f ( x)dx =
2 − ∫1 1 − z 2 F ( z )dz , trong đó
⎛a+b b−a ⎞
F ( z) = 1 − z 2 f ⎜ + z⎟ (1.32)
⎝ 2 2 ⎠
và
b
π n
∫ f ( x)dx ≅ (b − a)∑ 1 − z 2 f ( xi ) (1.33)
a
2n i =1
a+b b−a
xi = + zi
2 2
⎛ [2(i − 1) + 1]π ⎞
z i = cos⎜ ⎟
⎝ 2n ⎠
Các nút trong công thức Tchebyshev được trình bày tại bảng 1.3.
14
Bảng 1.3
n Wi X
3 2/3 X1 = -X3 = 0,707107
X2 = 0,0
4 ½ X1 = -X4 = 0,794654
X2 = -X3 = 0,187592
5 2/5 X1 = -X5 = 0,832498
X2 = -X4 = 0,374541
X3 = 0,0
1.5 TÌM NGHIỆM PHƯƠNG TRÌNH BẰNG PHƯƠNG PHÁP CHIA ĐÔI ĐOẠN CÓ
NGHIỆM
Phương pháp tìm nghiệm phương trình y = f(x) giản đơn mà nhanh, được dùng tại đây mang
tên gọi bisection ( lưỡng đoạn).
Giả sử phương trình f(x) = 0, trong đó f(x) liên tục trong đoạn [a, b] và f(a).f(b) < 0. Để tìm
nghiệm của phương trình vừa nêu, nằm trong [a, b], cần thiết chia đoạn (a, b) ra làm hai phần bằng
⎛a+b⎞
nhau. Nếu f ⎜ ⎟ = 0 , x = (a + b)/2 sẽ là nghiệm chính xác của phương trình. Ngược lại
⎝ 2 ⎠
⎛a+b⎞
f⎜ ⎟ ≠ 0 , chúng ta phải chọn miền có nghiệm [a, (a+b)/2] hay là [ (a+b)/2, b] theo dấu hiệu,
⎝ 2 ⎠
⎛a+b⎞
giá trị f(x), trong đó x lần lượt mang giá trị a, b, có dấu ngược với giá trị f ⎜ ⎟ , vừa nêu. Trong
⎝ 2 ⎠
phân đoạn mới, ký hiệu [a1, b1], bằng ½ giá trị đoạn [a, b] công việc được lặp lại như vừa trình bày.
Trường hợp chưa xác định được nghiệm tiến hành tiếp các thủ tục cho phân đoạn vừa được chọn
[a2, b2], …, [an, bn ], . . ., vv… thoả mãn điều kiện:
f(an) f(bn) < 0 n = 1, 2, . . . (a)
1
và bn − a n = (b − a ) (b)
2n
Có thể thấy rằng, a1, a1, . . ., an, . . . tạo thành chuỗi không giảm, còn giới hạn bên phải b1, b2,
. . ., bn, . . . tạo thành chuỗi không tăng, phải tồn tại một giới hạn mà tại đó hai chuỗi này tiến đến:
ξ = lim a n = lim bn
n→∞ n →∞
Phép chia được thực hiện với điều kiện (a), nếu n → ∞ điểm chọn x = (an+bn)/2 sẽ đạt đến
ξ, tại đó f(ξ) = 0.
Dưới đây trình bày một hàm viết bằng ngôn ngữ C giúp cho việc tìm nghiệm hàm f(x) được
định nghĩa tại (*func) (float). Số lần chia đôi đoạn chứa nghiệm MAX=40. Hàm nên được hiệu
chỉnh tùy trường hợp sử dụng.
/* Bisection Method */
#include
#define Max 40
15
float rtbis( func, x1, x2, xacc)
float x1, x2, xacc;
float (*func) ();
{
int j;
float dx, f, fmid, xmid, rtb;
void nerror();
f=(*func) (x2);
if (f*fmid >= 0.0) nerror("Bisection Method);
rtb = f < 0.0 ? (dx=x2-x1, x2);
for (j=1; j n
∑a
i =0
ki i c = bk ; k = 0,1,2,..., m (f)
m
a ki = ∑ f i ( x j ) f k ( x j );
j =0
với m
bk = ∑ y j f k ( x j )
j =0
Ví dụ: Xác lập hàm y = c1 ex + c2e2x qua ba điểm (0,1) (1, -2) (2, -40) trên cơ sở phương pháp tổng
nhỏ nhất các bình phương sai số.
c1 e0 + c2 e0 - 1 = δ1
c1 e1 + c2 e2 - 1 = δ2
c1 e2 + c2 e4 - 1 = δ3
Theo cách làm nêu trên, có thể viết:
∂δ 1 ∂δ ∂δ
δ1 +δ2 2 +δ3 3 = 0
∂c1 ∂c1 ∂c1
∂δ 1 ∂δ ∂δ
δ1 +δ2 2 +δ3 3 = 0
∂c 2 ∂c 2 ∂c 2
∂δ 1 ∂δ ∂δ ∂δ 1 ∂δ ∂δ
Thay các giá trị = e 0 ; 2 = e1 ; 3 = e 2 ; = e0 ; 2 = e2 ; 3 = e4 vào hệ
∂c1 ∂c1 ∂c1 ∂c2 ∂c 2 ∂c 2
phương trình cuối có thể xác định : c1 ≈ 2; c2 ≈ -1.
Nếu sử dụng các ký hiệu vector và ma trận vào bài toán này, hệ phương trình nêu trên được
viết lại gọn hơn.
Các hệ số c sắp xếp tại vector {c} gồm n+1 thành phần, giá trị yj, j = 0, 1, 2, …, m xếp vào
{y}, ma trận [F] tập họp m+1 dòng, ứng với các gía trị của fk(xj) .
⎧c 0 ⎫ ⎧ y0 ⎫
⎪c ⎪ ⎪y ⎪
⎪ ⎪ ⎪ ⎪
{c} = ⎨ 1 ⎬ ;{ y} = ⎨ 1 ⎬ ;
⎪•⎪ ⎪•⎪
⎪c n ⎪
⎩ ⎭ ( N +1) x1 ⎪ ym ⎪
⎩ ⎭ ( M +1) x1
⎡ f 0 ( x0 ) f1 ( x0 ) • f n ( x m )⎤
⎢ f (x ) f1 ( x1 ) • f n ( x m )⎥
⎢ 0 1 ⎥
⎢ • • • • ⎥
⎢ ⎥
⎢ • • • • ⎥
⎢ f 0 ( xm )
⎣ f1 ( xm ) • f n ( x m )⎥
⎦ ( M +1) x ( N +1)
Vector sai số xác định bằng biểu thức: [F]{c} – {y} = {δ},
Và công thức (d) trở thành:
17
∂
2[[ F ]{c} − { y}] [ F ]{c} = 0
∂{c}
hoặc : [[ F ]{c} − { y}][ F ] = 0
Vì rằng [F] ≠ {0}, biểu thức trong ngoặc vuông phải bằng 0. Bài toán theo phương pháp tổng
nhỏ nhất các bình phương đưa về dạng:
[F](m1xn1).{c}(n1x1) = {y}(m1x1), trong đó m1 = m + 1; n1 = n + 1. (g)
Ma trận [F] gồm m+1 dòng, n+1 cột, với n < m. Vector {c} chỉ có n+1 dòng, {y} gồm m+1
dòng.
Xử lý phương trình (g) theo nhiều cách khác nhau. Dưới đây trình bày hai cách thông dụng,
dễ dùng.
Cách thứ nhất đưa bài toán (g) về dạng bài toán kinh điển như trình bày tại (f). Nhân hai vế
của phương trình (g) với ma trận chuyển vị của [F] là [F]T :
[F]T[F]{c} =[F]T{y}
Từ đó có thể viết :
[A].{c} = {b}, (h)
T T
với [A] = [F] [F] và {b} = [F] {y}. (i)
Các hệ số {c} được xác định sau khi giải hệ phương trình đại số tuyến tính.
Ví dụ : áp dụng cách làm này xử lý bài toán vừa nêu.
⎡e 0 e0 ⎤ ⎧ 1 ⎫
2 ⎥⎧ 1 ⎫
⎢ 1 c ⎪ ⎪
⎢e e ⎥⎨ ⎬ = ⎨ − 2 ⎬
c
⎢e 2
⎣ e 4 ⎥ ⎩ 2 ⎭ ⎪− 40⎪
⎦ ⎩ ⎭
⎡e 0 e0 ⎤ ⎧ 1 ⎫
⎡e 0 e1
e 2
⎤⎢ 1 2 ⎥⎧ 1 ⎫
c ⎡e 0 e1 e 2 ⎤⎪ ⎪
⎢ 0 ⎥ e e ⎥⎨ ⎬ = ⎢ 0 4 ⎥⎨
−2 ⎬
⎣e e2 e4 ⎦⎢ 2 c e e2 e ⎦⎪
⎢e
⎣ e 4 ⎥⎩ 2 ⎭ ⎣
⎦
⎪
⎩− 40⎭
⎡ 63 424 ⎤ ⎧ c1 ⎫ ⎧ − 296 ⎫
⎢424 3036⎥ ⎨c ⎬ = ⎨− 2198⎬
⎣ ⎦⎩ 2 ⎭ ⎩ ⎭
Từ hệ phương trình cuối có thể xác định : c1 ≈ 2; c2 ≈ -1.
Cách làm thứ hai là tiến hành giải hệ phương trình đại số gồm m+1 phương trình, chứa chỉ n
+1 ẩn, trong đó n < m, bằng phương pháp xử lý ma trận suy biến. Phương pháp mang tên gọi bằng
tiếng Anh: Singular Value Decomposition – SVD.
Một trong những thuật toán hay nhất xử lý hệ phương trình đại số tuyến tính [A]
(MxN){b{(Nx1)={y}(Mx1) được trình bày trong sổ tay toán tính của Wilkinson. Thủ tục tính được nhóm
nhà toán học dưới sự chỉ dẫn của Wilkinson viết ra trong những năm sáu mươi bằng ngôn ngữ Algol
đã dùng có hiệu quả trên những dàn máy tính. Thủ tục trên người viết tài liệu này đã “dịch” sang
ngôn ngữ FORTRAN và sau đó “chuyển “ sang ngôn ngữ C, được giới thiệu tiếp theo, giúp bạn đọc
xử lý những bài toán thực tế. Trong chương trình, người viết đang hạn chế n người dùng cần tăng số hệ số ck đến số lớn hơn 9 ( tính từ 0), cần thay thế JMAX bằng số hệ số cao
nhất.
Hàm LeastSquares được dùng xác định các hệ số thủy động lực chân vịt tàu, xác định các hệ
số đường cong sức cản vỏ tàu trong phần tự động hoá thiết kế chân vịt. Hàm này là phương tiện
chính cho các phép hồi qui dùng trong phần xác định sức cản tàu.
#include
#include
#include
#include
#define JMAX 10
void LeastSquares( int N, int M, double A[][JMAX], double *b,
double *y )
{
register int N1, M1, MM1, L, i, LL, L1, j;
double s, ss, s2, d, pp;
MM1 = M -1; M1 = M+1;
for ( L =0; L < M; L++) {
ss =0.0;
for ( i = L; i < N; i++) ss += A[i][L] * A[i][L];
s2 = ss;
s = sqrt(s2);
if ( A[L][L] < 0.0 ) s = -s;
d = s2 + s*A[L][L];
A[L][L] += s;
if ( L != MM1) {
L1 = L + 1;
for ( j =L1; j < M; j++) {
pp = 0.0;
for ( i =L; i < N; i++) pp += A[i][L] * A[i][j];
A[N][j] = pp/ d;
}
for ( j = L1; j < M; j++)
for ( i =L; i < N; i++)
A[i][j] -= A[i][L] * A[N][j];
}
A[N][L] = -s;
}
for ( i=0; i
