Thủy lợi - Phương pháp tính

0 Trường Đại học Thủy lợi Phạm Phú Triêm phương pháp tính 1 Carl Friedrich Gauss (1777-1855) vua Toán học 2 MỤC LỤC Lời nói đầu ................................................................................................................................ 4 Chương 1: HỆ PHƯƠNG TRÌNH ĐẠI SỐ TUYẾN TÍNH................................................ 5 1.1 Phương pháp Cholesky .................................................................................................. 5 1.2 Phương pháp lặp Gauss-Seidel ..................................................................................... 8 1.3 Phương pháp nới lỏng.................................................................................................. 13 Chương 2 : PHƯƠNG TRÌNH VÀ HỆ PHƯƠNG TRÌNH PHI TUYẾN....................... 22 2.1 Phương pháp chia đôi ................................................................................................. 22 2.2 Phương pháp dây cung ................................................................................................ 25 2.3 Phương pháp tiếp tuyến............................................................................................... 28 2.4 Phương pháp lặp đơn................................................................................................... 31 2.5 Phương pháp Newton-Raphson cho hệ phương trình .............................................. 33 2.6 Phương pháp lặp Seidel cho hệ phương trình ........................................................... 37 Kiểm tra nhận thức ............................................................................................................ 42 Chương 3 : NỘI SUY GIÁ TRỊ HÀM SỐ ........................................................................... 44 3.1 Công thức nội suy Gregory-Newton tiến................................................................... 44 3.2 Công thức nội suy Gregory-Newton lùi...................................................................... 46 3.3 Công thức nội suy Gauss ............................................................................................. 48 3.4 Công thức nội suy Lagrange ....................................................................................... 50 3.5 Công thức nội suy Newton........................................................................................... 51 3.6 Công thức bình phương nhỏ nhất............................................................................... 54 Bài tập.................................................................................................................................. 57 Kiểm tra nhận thức ............................................................................................................ 58 Chương 4 : XẤP XỈ ĐẠO HÀM VÀ TÍCH PHÂN XÁC ĐỊNH....................................... 59 4.1 Xấp xỉ giá trị đạo hàm theo tỷ sai phân..................................................................... 59 4.2 Xấp xỉ giá trị đạo hàm theo công thức Richardson.................................................. 60 4.3 Xấp xỉ giá trị đạo hàm theo công thức nội suy với các mốc cách đều .................... 62 a- Công thức nội suy Gregory-Newton tiến.................................................................... 62 b- Công thức nội suy Gregory-Newton lùi ..................................................................... 62 c- Công thức nội suy Gauss ............................................................................................ 62 4.4 Xấp xỉ giá trị đạo hàm theo công thức nội suy với các mốc bất kỳ......................... 65 a- Công thức nội suy Lagrange ...................................................................................... 65 b- Công thức nội suy Newton.......................................................................................... 65 c- Công thức bình phương nhỏ nhất .............................................................................. 65 4.5 Xấp xỉ giá trị tích phân xác định A.... 68 a- Công thức hình thang ................................................................................................. 68 b- Công thức Simpson ..................................................................................................... 68 4.6 Dãy quy tắc .................................................................................................................. 71 a- Dãy quy tắc hình thang ............................................................................................... 71 b- Dãy quy tắc Simpson................................................................................................... 72 Bài tập.................................................................................................................................. 75 Kiểm tra nhận thức ............................................................................................................ 78 Chương 5: XẤP XỈ NGHIỆM PHƯƠNG TRÌNH VI PHÂN............................................ 79 5.1 Xấp xỉ nghiệm phương trinh vi phân cấp một................................................................ 79 a- Phương pháp Euler ...................................................................................................... 79 b- Phương pháp Runge-Kutta bậc hai .............................................................................. 80 c- Phương pháp Runge-Kutta bậc bốn ............................................................................ 80 5.2 Xấp xỉ nghiệm hệ phương trình vi phân cấp một........................................................... 82 a- Phương pháp Euler....................................................................................................... 82 3 b- Phương pháp Runge-Kutta bậc bốn ............................................................................. 83 5.3 Xấp xỉ nghiệm phương trình vi phân cấp 2.................................................................... 85 5.4 Xấp xỉ nghiệm phương trình đạo hàm riêng ................................................................ 86 Kiểm tra nhận thức ............................................................................................................ 96 4 Lời nói đầu Nhu cầu nâng cao chất lượng đào tạo sinh viên và các bài toán thực tiễn rất đa dạng, phức tạp là đòi hỏi cấp thiết đưa vào môn học PHƯƠNG PHÁP TÍNH nhằm giúp cho sinh viên khối kỹ thuật- Kỹ sư tương lai, tiếp cận với cách giải gần đúng phương trình, hệ phương trình có đánh giá sai só . . . , kết hợp trên cơ sở làm quen và tự nâng cao khả năng lập trình bằng một ngôn ngữ thường được sử dụng, đó là PASCAL. Bộ môn Toán và tác giả trân trọng giới thiệu Giáo trình này và vô cùng cảm ơn các ý kiến đóng góp quý giá của độc giả. Hà nội 8-2005 5 Chương 1: HỆ PHƯƠNG TRÌNH ĐẠI SỐ TUYẾN TÍNH Chúng ta biết rằng hầu hết các bài toán thực tế đều dẫn đến giải hệ phương trình tuyến tính. Các phương pháp giải hệ được trình bày khá đầy đủ trong các tài liệu về Đại số tuyến tính. Tuy nhiên mọi cách giải về mặt thực tiễn đều phải tính gần đúng gắn liền với sai số. Vì vậy trong chương này chúng ta giải quyết vấn đề nêu ra trên đây. 1.1 Phương pháp Cholesky Giải hệ AX = B (1.1.1) trong đó A = 11 1 1 ... ... ... ... n n nn a a a a é ù ê ú ê ú ê ú ê ú ê úë û , X = 1 ... n x x é ù ê ú ê ú ê ú ê úê úë û , B = 1 ... n b b é ù ê ú ê ú ê ú ê úê úë û (1.1.2) Các bước cơ bản của phương pháp này là ** Tìm 2 ma trận tam giác L = 21 31 32 1 2 3 1 1 0 0 ... 0 0 1 0 ... 0 0 1 ... 0 0 ... ... ... ... ... ... ... 1n n n nn l l l l l l l - é ù ê ú ê ú ê ú ê ú ê ú ê ú ê ú ê ú ê ú ê úë û , U = 11 12 13 1 1 1 22 23 2 1 2 33 3 1 3 ... 0 ... 0 0 ... ... ... ... ... ... ... 0 0 0 ... 0 n n n n n n nn u u u u u u u u u u u u u - - - é ù ê ú ê ú ê ú ê ú ê ú ê ú ê ú ê ú ê ú ê úë û (1.1.3) thoả mãn LU = A (1.1.4) ** Giải hệ sau đây tìm Y LY= B (1.1.5) ** Tìm nghiệm X từ hệ UX= Y (1.1.6) Đó chính là nghiệm của hệ (1.1.1). Thật vậy AX = L(UX) = LY = B. Để tìm L , U ta tìm lần lượt hàng 1 của U , cột 1 của L , hàng 2 của U , cột 2 của L , . . . theo công thức tổng quát sau đây 6 1 1 1 1 , , i i i j i j i k k j k i j i k k j k i j j j u a l u i j a l u l i j u - - = = ìïïï = - £åïïïïïíï - åïïï = >ïïïïî (1.1.7) Để đơn giản ta chỉ xét n = 3 và các bước tìm nghiệm * Tìm L và U L = 21 31 32 1 0 0 1 0 1 l l l é ù ê ú ê ú ê ú ê ú ê ú ê úë û , U = 11 12 13 22 23 33 0 0 0 u u u u u u é ù ê ú ê ú ê ú ê ú ê ú ê úë û (1.1.8) Hàng 1 của U là u11 = a11 , u12 = a12 , u13 = a13 Cột 1 của L là 2121 11 a l u = , 3131 11 a l u = Hàng 2 của U là u22 = a22 – l21u12 u23 = a23 – l21u13 Cột 2 của L là 32 31 1232 22 a l u l u - = Hàng 3 của U là u33 = a33 – (l31u13 + l32u23) * Tìm Y y1 = b1 ∧ y2 = b2 – l21y1 ∧ y3 = b3 – ( l31y1 + l32y2) * Tính nghiệm X 33 33 y x u = 2 23 32 22 y x u x u - Ù = 1 12 2 13 31 11 y ( ) x u x u x u - + Ù = Ví dụ 1.1 Giải hệ sau đây bằng phương pháp Cholesky 1 3 1 2 3 1 2 2 3 3 2 5 10 0 x x x x x x x ìï + =ïïï + + =íïïï - =ïî Giải A = 2 0 1 3 2 5 1 1 0 é ù ê ú ê ú ê ú ê ú-ê úë û , X = 1 2 3 x x x é ù ê ú ê ú ê ú ê úê úë û , B = 3 10 0 é ù ê ú ê ú ê ú ê úê úë û . * Tìm L và U Hàng1 của U là u11 = a11 = 2 , u12 = a12 = 0 , u13 = a13 = 1 Cột 1 của L là 2121 11 3 2 a l u = = , 3131 11 1 2 a l u = = 7 Hàng 2 của U là u22 = a22 – l21u12 = 2 – 32 .0 = 2 ; u23 = a23 – l21u13 = 5 – 3 2 .1 = 7 2 Cột 2 của L là 32 31 1232 22 11 .0 12 2 2 a l u l u - -- = = = - Hàng 3 của U là u33 = a33 – (l31u13 + l32u23) = 0 – 51 1 7.1 . 42 2 2 æ ö÷ç - =÷ç ÷çè ø 1 1 7 5.1 . 2 2 2 4 æ ö÷ç - =÷ç ÷çè ø * Tìm Y y1 = b1 = 3 ∧ y2 = b2 – l21y1 = 10 – 32 .3 = 112 ∧ y3 = b3 – ( l31y1 + l32y2) = 0 – 51 1 11.3 . 42 2 2 æ ö÷ç - =÷ç ÷çè ø * Tính nghiệm X 33 33 y x u = = 5 4 5 4 = 1 2 23 32 22 y x u x u - Ù = = 11 7.1 2 2 2 - = 1 1 12 2 13 31 11 y ( ) x u x u x u - + Ù = = 3 (0.1 1.1) 2 - + = 1 Program CHOLESKY; Var A:Array[1..10,1..10] of real; B,X,Y:Array[1..10] of real; R:Array[1..10] of integer; DK,I,J,K,T,N:Integer; DET,SUM:Real; BEGIN Write('So an N=');Readln(N); Writeln('Nhap ma tran A '); For I:=1 to N do For J:=1 to N do Begin Write('A[',I,J,']=');Readln(A[I,J]); End; Writeln(' Nhap vecto B '); For I:=1 to N do Begin Write('B[',I,']=');Readln(B[I]); End; DK:=1;DET:=1; For I:=1 to N do R[I]:=I; For I:=1 to N-1 do If DK=1 then Begin For K:=I+1 to N do Begin If ABS(A[R[K],I])>ABS(A[R[I],I]) then 8 Begin T:=R[I];R[I]:=R[K];R[K]:=T;DET:==-DET; End; If A[R[I],I]=0 then DK:=2 Else Begin DET:=DET*A[R[I],I]; For K:=I+1 to N do Begin A[R[K],I]:=A[R[K],I]/A[R[I],I]; For J:=I+1 to N do A[R[K],J]:=A[R[K],J]-A[R[K],I]*A[R[I],J]; End; End; End; End; If A[N,N]=0 then DK:=2 Else Begin Y[1]:=B[R[1]]; For K:=2 to N do Begin SUM:=0; For J:=1 to K-1 do SUM:=SUM+A[R[K],J]*Y[J]; Y[K]:=B[R[K]]-SUM; End; X[N]:=Y[N]/A[R[N],N]; For K:=N-1 downto 1 do Begin SUM:=0; For J:= K+1 to N do SUM:=SUM+A[R[K],J]*X[J]; X[K]:=(Y[K]-SUM)/A[R[K],K]; End; End; If DK=1 then For I:=1 to N do Write('X[',I,']=',X[I]); If DK=2 then Writeln('Khong tinh duoc'); Readln; END. 1.2 Phương pháp lặp Gauss-Seidel Giải hệ 9 11 1 1 1 1 1 ... ... ... ... ... ... n n n nn n n a x a x b a x a x b ìï + + =ïïïíïïï + + =ïïî (1.2.1) * Điều kiện lặp: A có đường chéo trội, có nghĩa là 1; n a ai i i j j j i ³ å = ¹ i = 1 , . . . , n (1.2.2) * Điều kiện hội tụ: A xác định dương, có nghĩa là a11 > 0 , 11 12 13 11 12 21 22 23 21 22 31 31 33 0 , 0 , . . . , 0 a a a a a a a a A a a a a a > > > (1.2.3) * Đưa hệ (1.2.1) về hệ tương đương 1 12 2 1 1 2 21 1 23 3 2 2 1 1 2 2 1 1 '... '... . . . . . . . . . . . . . . . . . . . . . . . . . . . '... n n n n n n n n n n n n n x c x c x b x c x c x c x b x c x c x c x b- - - - ìïï = + + +ïïïïï = + + + +ïíïïïïïï = + + + +ïïî (1.2.4) bằng cách: chuyển các xK ( k ≠ i ) của phương trình thứ i sang vế phải và chia 2 vế cho aii , i = 1, . . . , n * Thay xấp xỉ ban đầu X(0) = (x1(0) . . . xn(0)) của nghiệm vào (1.2.4) để tìm xấp xỉ thứ nhất X(1) = (x1(1) . . . xn(1)) của nghiệm (1) (0) (0) 12 1 11 2 (1) (1) (0) (0) 21 23 2 22 1 3 (1) (1) (1) (0) (0) 31 32 34 3 33 1 2 4 (1) (1) (1) (0) 11 1 2 1 11 1 2 '... '... '... . . . . . . . . . . . . . . . . . . . . . . . . . . . ... n n n n n n n n n n n n nn n x c x c x b x c x c x c x b x c x c x c x c x b x c x c x c x b- - - - -- - = + + + = + + + + = + + + + + = + + + + (1) (1) (1)(1) 1 2 11 2 1 ' '...n n n n n n nn nx c x c x c x b- -- - ìïïïïïïïïïïïïïïíïïïïïïïïïïïï = + + + +ïïî (1.2.5) Như vậy ta đã sử dụng các giá trị x1(1), . . . , xi(1) được tính ở các hàng trên để tìm các giá trị xi+1(1), . . . , xn(1) , i =1, . . . , n – 1, ở các hàng dưới. * Nếu (1) (0) 1 n i i i x x = å - ≤ εx (1.2.6) trong đó εx> 0 , cho trước (sai số đối với nghiệm) thì ta dừng ở X(1) . * Ngược lại, có nghĩa là 10 (1) (0) 1 n i i i x x = å - > εx (1.2.7) ta lại thay xấp xỉ thứ nhất X(1) = (x1(1) . . . xn(1)) của nghiệm vào (1.2.4) để tìm xấp xỉ thứ hai X(2) = (x1(2) . . . xn(2)) của nghiệm (2) (1) (1) 12 1 11 2 (2) (2) (1) (1) 21 23 2 22 1 3 (2) (2) (2) (1) (1) 31 32 34 3 33 1 2 4 (2) (2) (2) (1) 11 1 2 1 11 1 2 '... '... '... . . . . . . . . . . . . . . . . . . . . . . . . . . . ... n n n n n n n n n n n n nn n x c x c x b x c x c x c x b x c x c x c x c x b x c x c x c x b- - - - -- - = + + + = + + + + = + + + + + = + + + + (2) (2) (2)(2) 1 2 11 2 1 ' '...n n n n n n nn nx c x c x c x b- -- - ìïïïïïïïïïïïïïïíïïïïïïïïïïïï = + + + +ïïî (1.2.8) * Một cách tổng quát Ta tiếp tục cho đến xấp xỉ thứ k : X(k) = (x1(k) . . . xn(k)) của nghiêm theo công thức ( ) ( 1) ( 1) 12 1 11 2 ( ) ( ) ( 1) ( 1) 21 23 2 22 1 3 ( ) ( ) ( ) ( 1) ( 1) 31 32 34 3 33 1 2 4 ( ) ( ) ( ) 11 1 21 1 2 '... '... '... . . . . . . . . . . . . . . . . . . . . . . . . . . . ... k k k n n k k k k n n k k k k k n n k k k n n n nn n x c x c x b x c x c x c x b x c x c x c x c x b x c x c x c - - - - - - - - - -- - = + + + = + + + + = + + + + + = + + + ( 1)1 1 ( ) ( ) ( )( ) 1 2 11 2 1 ' '... k n n n k k kk n n n n n n nn n x b x c x c x c x b - - - -- - ìïïïïïïïïïïïïïïíïïïïïï +ïïïïïï = + + + +ïïî (1.2.9) và dừng lại ở X(k) khi thoả mãn ( 1)( ) 1 n kk i i i x x - = å - ≤ εx (1.2.10) hoặc ( 1) ( 2) 1 n k k i i i x x x k M e- - = ìïï - >ï åïíïïï =ïî (1.2.11) trong đó M là số bước lặp tối đa cho trước mà ta phải dừng lại mặc dù tại bước lặp thứ M – 1 vẫn chưa đảm bảo sai số. Ví dụ 1.2 Giải hệ sau đây bằng phương pháp Gauss-Seidel 1 2 3 1 2 3 1 2 3 5 2 6 4 2 2 4 1 x x x x x x x x x ìï - + =ïïï - + = -íïïï - + =ïî với X(0) = (0 0 0) , εx = 0,15 và M = 10. Giải 11 * Kiểm tra điều kiện đường chéo trội ⎜a11 ⎜= ⎜5 ⎜> ⎜– 1 ⎜ + ⎜2 ⎜= ⎜a12 ⎜ + ⎜a13 ⎜ ⎜a22 ⎜= ⎜– 5 ⎜> ⎜1 ⎜ + ⎜1 ⎜= ⎜a21 ⎜ + ⎜a23 ⎜ ⎜a33 ⎜= ⎜4 ⎜> ⎜– 2 ⎜ + ⎜– 1 ⎜= ⎜a31 ⎜ + ⎜a32 ⎜ * Kiểm tra điều kiện xác định dương: a11 = 5 > 0. Sau khi đổi dấu 2 vế của phương trình thứ 2 ta tính được 11 12 21 22 5 1 1 4 a a a a - = - = 19 > 0 11 12 13 21 22 23 31 31 33 5 1 2 1 4 1 2 1 4 a a a a a a a a a - = - - - - = 5.4.4 + ( – 1).( – 1).2 + ( – 2).( – 1).( – 1) – ( – 2).4.2 – ( – 1).( – 1).5 – ( – 1).( – 1).4 = 87 > 0 * Đưa hệ đã cho về hệ tương đương theo công thức (1.2.4) 1 2 3 2 1 3 3 1 2 0,2 0,4 1,2 0,25 0,25 0,5 0,5 0,25 0,25 x x x x x x x x x ìï = - +ïïï = + +íïïï = + +ïî * Thay xấp xỉ ban đầu X(0) = (0 0 0) của nghiệm vào vế phải tìm xấp xỉ thứ nhất X(1) = (x1(1) x2(1) x3(1)) của nghiệm theo công thức (1.2.5) với n = 3 (1) 1 (1) 2 (1) 3 0,2.0 0,4.0 1,2 1,2 0,25.1,2 0,25.0 0,5 0,8 0,5.1,2 0,25.0,8 0,25 1,05 x x x ìï = - + =ïïïï = + + =íïïïï = + + =ïî * Tính sai số 3 (1) (0) 1 i i i x x = å - = 1,2 + 0,8 + 1,05 = 3,05 > 0,15 * Tìm xấp xỉ thứ hai X(2) = (x1(2) x2(2) x3(2)) của nghiệm theo công thức (1.2.8) (2) 1 (2) 2 (2) 3 0,2.0,8 0,4.1,05 1,2 0,94 0,25.0,94 0,25.1,05 0,5 0,9975 0,5.0,94 0,25.0,9975 0,25 0,969375 x x x ìï = - + =ïïïï = + + =íïïïï = + + =ïî * Tính sai số 3 (2) (1) 1 i i i x x = å - = ⎜0,94 – 1,2⎜ + ⎜0,9975 – 0,8⎜ + ⎜0,969375 – 1,05⎜ = 0,538125 > 0,15 * Tìm xấp xỉ thứ ba X(3) = (x1(3) x2(3) x3(3)) của nghiệm theo công thức (1.2.9) : (3) 1 (3) 2 (3) 3 0,2.0,9975 0,4.0,969375 1,2 1,01175 0,25.1,01175 0,25.0,969375 0,5 0,995281 0,5.1,01175 0,25.0,995281 0,25 1,004695 x x x ìï = - + =ïïïï = + + »íïïïï = + + »ïî * Tính sai số 12 3 (3) (2) 1 i i i x x = å - = ⎟1,01175 – 0,94⎟ + ⎟0,995281 – 0,9975⎟ + ⎟1,004695 – 0,969375⎟ = 0,109289 < 0,15 Như vậy ta dừng ở xấp xỉ thứ 3 của nghiệm với độ chính xác cho trước X(3) = ( 1,01175 0,995281 1,004695) Program G_SEIDEL; Var AB:Array[1..20,1..20] of real; X0,XK:Array[1..20] of real; I,J,K,M,N,KMAX,DK:Integer; R,EPS,SUM,SUMEPS:Real; BEGIN Write('So an N = ');Readln(N); M:=N+1; Writeln('Nhap ma tran AB : N hang, N+1 cot ’); Writeln('N cot dau la cac cot cua ma tran A’); Writeln('cot thu N+1 la cot cac he so tu do B ’); For I:=1 to N do For J:=1 To M do Begin Write('AB[',I,J,']=');Readln(AB[I,J]); end; Write(' So buoc lap toi da KMAX=');Readln(KMAX); Write(' Sai so EPS =');Readln(EPS); DK:=1;I:=1; { Kiem tra tinh cheo troi } While ((I<=N) and (DK=1)) do Begin R:=0; For J:=1 to N do R:=R + ABS(AB[I,J]); If R>=2*ABS(AB[I,I]) then DK:=2;I:=I+1: end; If DK=1 then { Khi tinh cheo troi thoa man ta nhap xap xi ban dau } Writeln(‘ Nhap xap xi ban dau X’); Begin For I:=1 to N do Begin Write('X0[',I,']=');Readln(X0[I]); end; end; 13 K:=1; While ((K<KMAX) and (DK=1)) do { Tinh cac thanh phan XK[I] o buoc lap thu K khi khong su dung cac gia tri XK[1] , XK[2] , . . . , XK[I] da tinh o tren de tinh XK[I+1], . . . ,XK[N] ; I = 1 , . . . , N– 1 } Begin For I:=1 to N do Begin SUM:=AB[I,M]; For J:=1 to N do If JI then SUM:=SUM – AB[I,J]*X0[J];XK[I]:=SUM/AB[I,I]; end; { So sanh tong sai so cac thanh phan o buoc lap thu K va K-1 voi sai so da cho } SUMEPS:=0; For I:=1 to N do SUMEPS:=SUMEPS + ABS(XK[I]-X0[I]); If SUMEPS<EPS then DK:=3; For I:=1 to N do X0[I]:=XK[I];K:=K+1; end; If DK=3 then Begin Writeln('Nghiem xap xi la'); For I:=1 to N do Writeln('X',K,'[',I,']=',XK[I]); end; If DK=2 then Writeln('Khong cheo troi'); If K=KMAX then Begin Writeln('Chua hoi tu va ta dung o xap xi thu k cua nghiem,k= ', K); For I:=1 to N do Writeln('X',K,'[',I,']=',XK[I]); end; Readln; END. 1.3 Phương pháp nới lỏng Giải hệ 14 11 1 1 1 1 1 ... ... ... ... ... ... n n n nn n n a x a x b a x a x b ìï + + =ïïïíïïï + + =ïïî (1.3.1) * Điều kiện lặp: A có đường chéo trội, có nghĩa là 1; n a ai i i j j j i ³ å = ¹ i = 1 , . . . , n (1.3.2) * Điều kiện hội tụ: A xác định dương, có nghĩa là a11 > 0 , 11 12 13 11 12 21 22 23 21 22 31 31 33 0 , 0 , . . . , 0 a a a a a a a a A a a a a a > > > (1.3.3) * Đưa hệ