NGÔN NGỮ LẬP TRÌNH FORTRAN VÀ ỨNG DỤNG TRONG KHÍ TƯỢNG THỦY VĂN part 10

Chia sẻ: Ashfjshd Askfaj | Ngày: | Loại File: PDF | Số trang:4

Thêm vào BST

Báo xấu

143
lượt xem 20
download

Download Vui lòng tải xuống để xem tài liệu đầy đủ

Những thí dụ và hệ thống bài tập tự luyện trong sách này có ý nghĩa minh họa, hướng sinh viên tới vận dụng các lệnh của Fortran để viết ra những chương trình ứng dụng nho nhỏ có tính cụ thể, bước đầu làm quen với những đặc thù xử lý dữ liệu quan trắc trong chuyên môn khí tượng thủy văn. Những đặc điểm khác của nội dung ứng dụng lập trình trong các chuyên ngành này như quản lý cơ sở dữ liệu, các phương pháp thống kê hiện đại, các phương pháp giải số trị những bài toán...

Chủ đề:

Bình luận(0) Đăng nhập để gửi bình luận!

Lưu

Nội dung Text: NGÔN NGỮ LẬP TRÌNH FORTRAN VÀ ỨNG DỤNG TRONG KHÍ TƯỢNG THỦY VĂN part 10

n n n ∑x ∑y − n ∑ xk y k k k k =1 k =1 k =1 a= (20) n n (∑ x k ) − n ∑ x k2 2 k =1 k =1 n n n n ∑x ∑x y k − ∑ x k2 ∑y k k k b= k =1 k =1 k =1 k =1 , (21) Phụ lục 3: Phương pháp bình phương nhỏ nhất n n (∑ x k ) 2 − n ∑ x k2 trong phân tích hồi quy k =1 k =1 hay hệ số b còn có thể tính theo công thức: n n ∑y ∑x −a k k 1. Mô hình tuyến tính k= k= b= 1 1 . (22) n Mô hình hồi quy tuyến tính có dạng: 2. Mô hình đa thức y = f ( x) = ax + b . Phương pháp bình phương nhỏ nhất cũng có thể áp dụng để tính các Theo phương pháp bình phương nhỏ nhất, các hệ số hồi quy a và b hệ số hồi quy đa thức dạng trong phương trình trên được tìm sao cho tổng bình phương sai số bằng n f ( x) = a 0 + a1 x + a 2 x 2 + ... + a n x m . E = ∑ ( y k − axk − b) 2 k =1 thí dụ đối với mô hình bậc hai 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 : f ( x) = a 0 + a1 x + a 2 x 2 . n n n n n a ∑ xk2 + b∑ xk = ∑ xk y k , a ∑ xk + b n = ∑ y k . 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: k =1 k =1 k =1 k =1 k =1 Vậy các hệ số hồi quy được tính theo các công thức sau: 108
(29) ⎧ n n n ∑ xk + a1 ∑ xk + a0 n = ∑ yk 2 ⎪ a2 ta có ⎪ k=1 k=1 k=1 ⎪ ~ ⎪ g ( x ) = a~ + b n n n n (23) ∑ xk + a1 ∑ xk + a0 ∑ xk = ∑ xk yk x 3 2 (30) ⎨ a2 ⎪ k=1 k=1 k=1 k=1 ~ ⎪ Với phương trình (30) các hệ số hồi quy a và b tính theo các công n n n n ∑ xk + a1 ∑ xk + a0 ∑ xk = ∑ xk yk 4 3 2 2 ⎪ a2 thức ⎪ ⎩ k=1 k=1 k=1 k=1 n n n Về nguyên tắc ta có thể sử dụng phương pháp này để tìm phương ∑ ~k ∑ ~k − n ∑ ~k ~k x y xy trình đa thức bậc bất kỳ. Tuy nhiên trong thực tế phương pháp trở thành a= k =1 k =1 k =1 (31) n n 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 (∑ ~k ) 2 − n ∑ ~k2 x x tính. k =1 k =1 n n n n ∑~ ∑~ ~ − ∑~ ∑~ 3. Mô hình phi tuyến 2 x xy x y k k k k k ~ b= k =1 k =1 k =1 k =1 (32) Phương pháp bình phương nhỏ nhất có thể áp dụng cho hàm bất kỳ, n n (∑ ~k ) 2 − n ∑ ~k2 x x nhưng hệ các phương trình để tìm các hệ số có thể phi tuyến, và do đó k =1 k =1 không thể giải được bằng cách sử dụng các phương trình tuyến tính. Tuy Vậy công việc tính toán gồm: chuyển đổi các giá trị số liệu x k và y k nhiên, trong một số trường hợp, một hàm phi tuyến có thể chuyển thành 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 một hàm tuyến tính. Thí dụ về một hàm có thể tuyến tính hoá là ~ trình (31), (32) để tìm a và b . Giải phương trình (27) đối với b và đặt vào f ( x) = b x a (24) phương trình (24). Nếu lấy loga hai vế của phương trình này, ta có ln f ( x) = a ln x + ln b . (25) Nếu ký hiệu (26) g ( x) = ln f ( x) ~ b = ln b (27) ~ = ln x (28) ~ = ln y x y 109
[x1 ] [x2 ] [xm ] ⎛n ⎞ ⎛ a0 ⎞ ⎛ b0 ⎞ ... ⎜ ⎟ ⎜⎟⎜⎟ ⎜ [x1 ] [x1 x1 ] [x2 x1 ] [xm x1 ] ⎟ ... ⎜ a1 ⎟ ⎜ b1 ⎟ ⎜ [x ] [x1 x2 ] [x2 x2 ] [xm x2 ] ⎟ . ⎜ a 2 ⎟ = ⎜ b2 ⎟ ... ⎜2 ⎟ ⎜⎟⎜⎟ ⎜ ... ... ⎟ ⎜ ... ⎟ ⎜ ... ⎟ ... ... ... ⎜ [xm xm ]⎟ ⎜⎟⎜⎟ ⎝ [x m ] [x1 xm ] [x2 xm ] ... ⎝ a m ⎠ ⎝ bm ⎠ ⎠ Phụ lục 4: Sơ đồ ứng dụng phương pháp hồi (34) quy nhiều biến n [ ] ký hiệu phép lấy tổng ∑ . với dấu 1 Giả sử có n quan trắc đối với biến phụ thuộc y và các biến độc lập Để tìm các hệ số hồi quy a0 , a1 , a2 ,..., am ta phải giải hệ phương x1 , x 2 , ..., x m . Phương trình hồi quy được thiết lập như sau: 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 y = a 0 + a1 x1 + a 2 x 2 + ... + a m x m . 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 Các hệ số hồi quy ai (i = 1,..., m) được chọn sao cho thoả mãn 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 − n δ = ∑ ( y − a 0 − a1 x1 − a 2 x 2 − ... − a m x m )2 = min SUBROUTINE GAUSS. i =1 Lần lượt lấy đạo hàm biểu thức trên theo a0 , a1 , a 2 ,..., a m và cho các SUBROUTINE LHPTCT (Y, X, A, N, M) INTEGER N, M, I, J, K đạo hàm bằng không, ta có hệ m + 1 phương trình để xác định các hệ số a REAL Y (10000), X (10000, 50), A (0 : 50, 0 : 51) [x1 ]a1 [x 2 ]a 2 [x m ]a m [y ] + + + ... + = na 0 A (0, 0) = N [x1 ]a 0 [x1 x1 ]a1 [x 2 x1 ]a 2 [x m x1 ]a m [yx1 ] + + + ... + = DO J = 1, M [x 2 ]a 0 [x1 x 2 ]a1 [x 2 x 2 ]a 2 [x m x 2 ]a m [yx 2 ] + + + ... + = A (0, J) = 0.0 DO K = 1, N ... ... ... ... ... ... A (0, J) = A (0, J) + X (K, J) [x m ]a 0 + [x1 x m ]a1 + [x 2 x m ]a 2 [x m x m ]a m = [ yx m ] + ... + END DO END DO (33) A (0, M + 1) = 0.0 Hệ phương trình này gọi là hệ phương trình chính tắc để xác định các DO K = 1, N hệ số hồi quy. Dưới dạng ma trận ta viết hệ này như sau: A (0, M + 1) = A (0, M + 1) + Y (K) END DO 110
END IF DO I = 1, M END DO A (I, M + 1) = 0.0 IF (K .NE. I) THEN DO K = 1, N DO J = I, M + 1 A (I, M + 1) = A (I, M + 1) + Y (K) * X(K, I) AMAX = A (I, J) END DO A (I, J) =A (K, J) END DO A (K, J) = AMAX DO I = 1, M END DO DO J = I, M END IF A (I, J) = 0.0 DO J = I + 1, M + 1 DO K = 1, N A (I, J) = A (I, J) / A (I, I) A (I, J) = A (I, J) + X (K, I) * X (K, J) END DO END DO DO J = I + 1, M ENDDO DO K = I + 1, M + 1 ENDDO A (J, K) = A (J, K) - A (J, I) * A (I, K) DO I = 1, M END DO DO J = 0, I - 1 END DO A (I, J) = A (J, I) END DO END DO X (M) = A (M, M + 1) / A (M, M) END DO DO I = M - 1, 0, -1 RETURN X (I) = A (I, M + 1) END DO J = I + 1, M X (I) = X (I) - A (I, J) * X (J) SUBROUTINE GAUSS (M, A, X) END DO INTEGER M END DO REAL A (0 : 50, 0 : 51), X (0 : 50) RETURN DO I = 0, M - 1 END K=I AMAX = ABS (A (K, K)) DO J = I + 1, M R = ABS (A (J, I)) IF (AMAX .LT. R) THEN AMAX = R K=J 111