Giáo trình Phương pháp tính - Phạm Phú Triêm

100 trang cucquyet12 2900

Download

Bạn đang xem 20 trang mẫu của tài liệu "Giáo trình Phương pháp tính - Phạm Phú Triêm", để tải tài liệu gốc về máy bạn click vào nút DOWNLOAD ở trên

Tài liệu đính kèm:

giao_trinh_phuong_phap_tinh_pham_phu_triem.pdf

Nội dung text: Giáo trình Phương pháp tính - Phạm Phú Triêm

Trường Đại học Thủy lợi Phạm Phú Triêm phương pháp tính 0
Carl Friedrich Gauss (1777-1855) vua Toán học 1
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 hKiểm tra n ậ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 hKiểm tra n ậ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 Kih ểm tra n ậ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 2
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 Kih ểm tra n ận thức 96 3
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 4
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 đó éùaa é ù éù êú11 1n êx 1 ú êúb 1 êú ê ú êú A = êú , X = ê ú , B = êú êú aa êx ú êúb ëûêú1nnn ëê n ûú ëûêún (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 é ù éùuu u u u ê 1 0 0 0 0ú êú11 12 13 1nn- 1 1 ê ú êú êl 21 1 0 0 0ú êú0 uu22 23 u 2nn- 1 u 2 ê ú êú L = êll 1 0 0ú , U = êú0 0uuu ê 31 32 ú êú33 3nn- 1 3 ê ú êú ê ú êú êú êll l l 1ú 0 0 0 0 u ëê nn123 n nn- 1ûú ëûêúnn (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 5
ì ï i- 1 ï uaij=- ijå luij ik kj , £ ï k = 1 ï í i- 1 ï aluij- å ik kj ï k = 1 ï lij=>, ï ij u îï jj (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 éù é ù êú100 êuu11 12 u 13 ú êú ê ú L = êúl 10 , U = ê 0 uuú êú21 ê 22 23ú êúll1 ê 00u ú ëûêú31 32 ëê 33ûú (1.1.8) Hàng 1 của U là u11 = a11 , u12 = a12 , u13 = a13 Cột 1 của L là a 21 a 31 l 21 = , l 31 = u 11 u 11 Hàng 2 của U là u22 = a22 – l21u12 u23 = a23 – l21u13 Cột 2 của L là alu32- 31 12 l 32 = u 22 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 y 3 y 2233- ux y(1122133-+ux ux ) x3 = Ù=x2 Ù=x1 u 33 u 22 u 11 Ví dụ 1.1 ïì 23xx+= ï 13 íï 32xxx++ 5 = 10 Giải hệ sau đây bằng phương pháp Cholesky ï 123 ï îï xx12-=0 Giải é ù é ù é ù ê201ú êx 1 ú ê3 ú ê ú ê ú ê ú A = ê325ú , X = êx 2 ú , B = ê10ú. ê ú ê ú ê ú ëê110- ûú ëêx 3 ûú ëê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à a 21 3 a 31 1 l 21 == , l 31 == u 11 2 u 11 2 6
Hàng 2 của U là u = a – l u = 2 – 3 .0 = 2 ; u = a – l u = 5 – 3 .1 = 7 22 22 21 12 2 23 23 21 13 2 2 Cột 2 của L là 1 1.0 alu32- 31 12 2 1 l 32 ===- u 22 22 Hàng 3 của U là æö117÷ 5 æö1175÷ u33 = a33 – (l31u13 + l32u23) = 0 – ç .1-= . ÷ ç .1-= . ÷ èøç222÷ 4 èøç2224÷ * Tìm Y y = b = 3 ∧ y = b – l y = 10 – 3 .3 = 11 1 1 2 2 21 1 2 2 æö1111÷ 5 ∧ y3 = b3 – ( l31y1 + l32y2) = 0 –ç .3-= . ÷ èøç222÷ 4 5 11 7 - .1 y 3 4 y 2233- ux 22 * Tính nghiệm X x3 = = = 1 Ù=x2 = = 1 u 33 5 u 22 2 4 y(1122133-+ux ux ) 3(0.11.1)-+ Ù=x1 = = 1 u 11 2 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 7
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ệ 8
ì ï ax11 1++ a 1nn x = b 1 ï íï ï ï ax++ a x = b îï 11nnnnn (1.2.1) * Điều kiện lặp: A có đường chéo trội, có nghĩa là n aaii³ å i j i = 1 , . . . , n jji=¹1; (1.2.2) * Điều kiện hội tụ: A xác định dương, có nghĩa là aaa11 12 13 aa11 12 a11 > 0 , >>>0, aaa21 22 23 0, ,A 0 aa21 22 aaa31 31 33 (1.2.3) * Đưa hệ (1.2.1) về hệ tương đương ïì ï xcx=+++ cxb' ï 1122 1nn 1 ï ï xcxcx=++++ cxb' íï 2211233 2nn 2 ï ï ï ï x =++cx c x + c x + b' îï nn11 nnnnnn 2 2 1 1 n (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 (0) (0) (0) * Thay xấp xỉ ban đầu X = (x1 . . . xn ) của nghiệm vào (1.2.4) để tìm xấp xỉ thứ nhất (1) (1) (1) X = (x1 . . . xn ) của nghiệm ì ï (1) (0) (0) ' ï xcx12=+++12 cxb1nn 1 ï ï (1) (1) (0) (0) ' ï xcxcx21=+21 23 3 +++ cxb2nn 2 ï ï (1) (1) (1) (0) (0) ' ï xcxcxcx=+ + +++ cxb íï 3131 32 2 34 4 3nn 3 ï ï ï (1) (1) (1) (0) ' ï xcxnn 11=+++nn 11 c1nn 2 x2 cxb1 nnn +1 ï ï xcx(1) =+++(1) c x(1) c x(1) + b' îï nn121121 nn nn nn n (1.2.5) (1) (1) Như vậy ta đã sử dụng các giá trị x1 , . . . , xi được tính ở các hàng trên để tìm các giá trị (1) (1) xi+1 , . . . , xn , i =1, . . . , n – 1, ở các hàng dưới. * Nếu n (1) (0) å xii- x ≤ εx i= 1 (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à 9
n (1) (0) å xii- x > εx i= 1 (1.2.7) (1) (1) (1) ta lại thay xấp xỉ thứ nhất X = (x1 . . . xn ) của nghiệm vào (1.2.4) để tìm xấp xỉ thứ hai (2) (2) (2) X = (x1 . . . xn ) của nghiệm ì ï (2) (1) (1) ' ï xcx12=+++12 cxb 1nn 1 ï ï (2) (2) (1) (1) ' ï xcxcx21=++++21 23 3 cxb2nn 2 ï ï (2) (2) (2) (1) (1) ' ï xcxcxcx=+ ++++ cxb íï 3131 32 2 34 4 3nn 3 ï ï ï (2) (2) (2) (1) ' ï xcxnn 11=+++nn 11 c1nn 2 x2 cxb1 nn +n1 ï ï xcx(2)=+++(2) c x(2) cx(2) + b' îï nn121121 nn nn nn n (1.2.8) * Một cách tổng quát (k) (k) (k) Ta tiếp tục cho đến xấp xỉ thứ k : X = (x1 . . . xn ) của nghiêm theo công thức ì ï ()kk (- 1) (1)k- ' ï xcx12=+++12 cxb1nn 1 ï ï ()kk () ( k- 1) (1)k- ' ï xcxcx21=+21 23 3 +++ cxb2nn 2 ï ï ()kk () () k ( k- 1) (1)k- ' ï xcxcxcx=+ + +++ cxb íï 3131 32 2 34 4 3nn 3 ï ï ï ()kk () () k (1)k- ' ï xcxnn 11=+++nn 11 c1nn 2 x2 c11nnxb+ n- ï ï xcx()k =+++()kkk c x() cx() + b' îï nn121121 nn nn nn n (1.2.9) và dừng lại ở X(k) khi thoả mãn n ()k (1)k- å xxii- ≤ εx i= 1 (1.2.10) hoặc ì n ï (1)(2)kk ï å xii->xxe í i= 1 ï îï kM= (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 ïì 526xx-+= x ï 123 íï xx-+=-42 x Giải hệ sau đây bằng phương pháp Gauss-Seidel ï 123 ï îï 241xx12-+ x 3 = với X(0) = (0 0 0) , εx = 0,15 và M = 10. Giải 10
* 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 aa 51- 11 12 = = 19 > 0 aa21 22 - 14 aaa11 12 13 512- aaa21 22 23 =-14 - 1= 5.4.4 + ( – 1).( – 1).2 aaa31 31 33 214 + ( – 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) ïì xxx=-0,2 0,4 + 1,2 ï 123 íï xxx=+0,25 0,25 + 0,5 ï 213 ï îï xx31=+0,5 0,25 x 2 + 0,25 * 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 (1) (1) (1) (1) X = (x1 x2 x3 ) của nghiệm theo công thức (1.2.5) với n = 3 ïì x(1) =-0,2.0 0,4.0 += 1,2 1,2 ï 1 ï (1) í x2 =++=0,25.1,2 0,25.0 0,5 0,8 ï ï (1) îï x3 =+0,5.1,2 0,25.0,8 += 0,25 1,05 * Tính sai số 3 (1) (0) å xii- x = 1,2 + 0,8 + 1,05 = 3,05 > 0,15 i= 1 (2) (2) (2) (2) * Tìm xấp xỉ thứ hai X = (x1 x2 x3 ) của nghiệm theo công thức (1.2.8) ïì x(2) =-0,2.0,8 0,4.1,05 += 1,2 0,94 ï 1 ï (2) í x2 =++=0,25.0,94 0,25.1,05 0,5 0,9975 ï ï (2) îï x3 =+0,5.0,94 0,25.0,9975 += 0,25 0,969375 * Tính sai số 3 (2) (1) å xii- x = ⎜0,94 – 1,2⎜ + ⎜0,9975 – 0,8⎜ + ⎜0,969375 – 1,05⎜ = i= 1 0,538125 > 0,15 (3) (3) (3) (3) * Tìm xấp xỉ thứ ba X = (x1 x2 x3 ) của nghiệm theo công thức (1.2.9) : ïì x(3) =-0,2.0,9975 0,4.0,969375 += 1,2 1,01175 ï 1 ï (3) í x2 =+0,25.1,01175 0,25.0,969375 +» 0,5 0,995281 ï ï (3) îï x3 =+0,5.1,01175 0,25.0,995281 +» 0,25 1,004695 * Tính sai số 11
3 (3) (2) å xxii- = ⎟1,01175 – 0,94⎟ + ⎟0,995281 – 0,9975⎟ i= 1 + ⎟1,004695 – 0,969375⎟ = 0,109289 =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; 12
K:=1; While ((K I 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ệ 13
ì ï ax11 1++ a 1nn x = b 1 ï íï ï ï ax++ a x = b îï 11nnnnn (1.3.1) * Điều kiện lặp: A có đường chéo trội, có nghĩa là n aaii ³ å i j i = 1 , . . . , n jji=¹1; (1.3.2) * Điều kiện hội tụ: A xác định dương, có nghĩa là aaa11 12 13 aa11 12 a11 > 0 , >>>0, aaa21 22 23 0, ,A 0 aa21 22 aaa31 31 33 (1.3.3) * Đưa hệ (1.3.1) về hệ mới tương đương bằng cách: chuyển tất cả các ẩn sang phải và chia phương trình thứ i cho aii , i = 1 , . . . , n ïì ï -+xcx ++ cxb +=' 0 ï 1122 1nn 1 ï ï cx-+ x cx ++ c x += b' 0 íï 21 1 2 23 3 2nn 2 ï ï ï ï cx++ c x - x + b' =0 îï nnnnnn11 1 1 (1.3.4) (0) (0) (0) * Thay xấp xỉ ban đầu X = (x1 . . . xn ) của nghiệm vào (1.3.4) tìm sai số xấp xỉ ban đầu (0) (0) (0) R = (R1 . . . Rn ) ì ï (0) (0) (0) ' (0) ï -+xcx1212 +++= cxbR1nn 1 1 ï ï (0) (0) (0) (0) ' (0) ï cxxcx21 12-+23 3 ++ cxbR2nn +=2 2 ï ï (0) (0) (0) (0) (0) ' (0) ï cx+-++++= c x x c x c x b R íï 311234 32 34 3nn 3 3 ï ï ï (0) (0) (0) (0) ' (0) ï cxnnn 11 12++ c1- 2 xnn x1++=cxbnnn 11 n Rn- 1 ï ï cx(0) ++ c x(0) + c x(0) - x(0) + b' = R(0) îï nn121121nn nn nnnn (1.3.5) (0) * Xác định Rh thoả mãn (0) (0) (0) ⎜Rh ⎜= max(⎜R1 ⎜, . . . , ⎜Rn ⎜) (1.3.6) * Nếu (0) ⎜Rh ⎜ ≤ εf (1.3.7) trong đó εf > 0 , cho trước (sai số đối với phương trình) thì ta dừng ở X(0). * Ngược lại, có nghĩa là 14
(0) ⎜Rh ⎜ > εf (1.3.8) (1) (1) (1) ta tìm xấp xỉ thứ nhất X = (x1 . . . xn ) của nghiệm theo công thức ì ï xx(1)=+ (0) R (0) íï hh h ï (1) (0) îï x ii==xi; 1, , nih ; ¹ (1.3.9) (1) (1) (1) và tính xấp xỉ thứ nhất R = (R1 . . . Rn ) của sai số theo công thức ïì (0) ï R h = 0 íï ï R (0)=+=¹cR (0) R (0) ; i 1, , ni ; h îï iiih h (1.3.10) (1) * Xác định Rh thoả mãn (1) (1) (1) ⎜Rh ⎜= max(⎜R1 ⎜, . . . , ⎜Rn ⎜) (1.3.11) * Nếu (1) ⎜Rh ⎜ ≤ εf (1.3.12) thì ta dừng ở X(1) . * Ngược lại, có nghĩa là (1) ⎜Rh ⎜ > εf (1.3.13) tương tự ta tìm xấp xỉ thứ 2 của nghiệm và xấp xỉ thứ 2 của sai số . . . * Một cách tổng quát (k) (k) (k) Ta tiếp tục cho đến xấp xỉ thứ k : X = (x1 . . . xn ) của nghiệm được tính theo công thức ì ï xx()k =+(1)kk R (1) íï hh h ï ()k (1)k- îï x ii==¹xinih; 1, , ; (1.3.14) (k) (k) (k) và tính xấp xỉ thứ k : R = (R1 . . . Rn ) của sai số theo công thức ïì (0) ï R h = 0 íï ï R ()k =+=¹cR(1)(1)kk R; i 1, , ni ; h îï iiih h (1.3.15) và dừng lại ở X(k) khi (k) (k) (k) ⎜Rh ⎜= max(⎜R1 ⎜, . . . , ⎜Rn ⎜) (1.3.16) thoả mãn (k) ⎜Rh ⎜≤ εf (1.3.17) hoặc 15
ì ï R(1)k- > e f íï h ï îï kM= (1.3.18) trong đó M là số bước lặp tối đa cho trước mà ta phải dừng mặc dù tại bước lặp thứ M – 1 vẫn chưa đảm bảo sai số. Nội dung tư tưởng của phương pháp nới lỏng: Biến đổi sao cho trong xấp xỉ tiếp theo, sai số tương ứng với phương trình thứ h (sai số xấp xỉ có trị tuyệt đối lớn nhất) sẽ bằng không. Điều này thể hiện bởi các công thức (1.3.14), (1.3.15). Ví dụ 1.3 ì ï 526xx123-+= x ï (0) ïí xxx-+=-42 Giải hệ ï 123 bằng phương pháp nới lỏng với X = (0 0 0) , εf = ï îï 241xx123 + x = 0,15 và M = 3. Giải * Kiểm tra điều kiện đường chéo trội : ⎜a11 ⎜= ⎜5 ⎜> ⎜– 1 ⎜ + ⎜2 ⎜= ⎜a12 ⎜ + ⎜a13 ⎜ ⎜a22 ⎜= ⎜– 4 ⎜> ⎜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 aa 51- 11 12 = = 19 > 0 aa21 22 - 14 aaa11 12 13 512- aaa21 22 23 =-14 - 1= 5.4.4 + ( – 1).( – 1).2 aaa31 31 33 214 + ( – 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.3.4) ïì -+xxx0,2 - 0,4 + 1,2 = 0 ï 123 íï 0,25xxx-++= 0,25 0,5 0 ï 123 ï îï 0,5xxx123 .+-+= 0,25 0,25 0 * Thay xấp xỉ ban đầu X(0) = (0 0 0) của nghiệm vào vế phải tìm sai số ban đầu (0) (0) (0) (0) R = (R1 R2 R3 ) theo công thức (1.3.5) với n = 3 ïì -+0 0,2.0 - 0,4.0 + 1,2 = 1,2 =R (0) ï 1 ï (0) í 0,25.0-+ 0 0,25.0 + 0,5 = 0,5 =R 2 ï ï (0) îï 0,5.0+-+== 0,25.0 0 0,25 0,25 R 3 * Xác định theo công thức (1.3.6) (0) ⎜R1 ⎜= max(1,2 ; 0,5 ; 0,25) = 1,2 > 0,15 16
cho nên ta chuyển sang bước sau với h = 1. (1) (1) (1) (1) * Tìm X = (x1 x2 x3 ) theo công thức (1.3.9) ïì xxR(1)=+ (0) (0) =+=01,21,2 ï 11 1 ï (1) (0) í xx22==0 ï ï (1) (0) îï xx33==0 * Tìm xấp xỉ thứ nhất của sai số theo (1.3.10) ïì (1) ï R 1 = 0 ï íï RcRR(1)=+=+= (0) (0) 0,25.1,2 0,5 0,8 ï 21221 ï ï RcRR(1)=+=+= (0) (0) 0,5.1,2 0,25 0,85 îï 31331 * Xác định theo công thức (1.3.11) (1) ⎜R3 ⎜= max(0 ; 0,8 ; 0,85) = 0,85 > 0,15 cho nên ta chuyển sang bước sau với h = 3. (2) (2) (2) (2) * Tìm X = (x1 x2 x3 ) theo công thức (1.3.14) ïì ï xxR(2)=+ (1) (1) =+=0 0,85 0,85 ï 333 ï (2) (1) í xx11= = 1,2 ï ï (2) (1) ï xx= = 0 îï 2 2 * Tìm xấp xỉ thứ hai của sai số theo (1.3.15) ïì (2) ï R 3 = 0 ï íï RcRR(2) =+=-+=-(1) (1) 0,4.0,85 0 0,36 ï 1 13 31 ï ï RcRR(2)=+=+= (1) (1) 0,25.0,85 0,8 1,0125 îï 23223 * Xác định theo công thức (1.3.16) (2) ⎜R2 ⎜= max(0 ; 0,36 ; 1,0125) = 1,0125 > 0,15 cho nên ta chuyển sang bước sau với h = 2. (3) (3) (3) (3) * Tìm X = (x1 x2 x3 ) theo công thức (1.3.14) ì (2) ï xxR(3) =+=+=(2) 0 1,0125 1,0125 ï 22 2 ï (3) (2) í xx11==1,2 ï ï xx(3)== (2) 0,85 îï 33 * Tìm xấp xỉ thứ ba của sai số theo (1.3.15) ïì (3) ï R 2 = 0 ï íï RcRR(3)=+=-=- (2) (2) 0,2.1,0125 0,36 0,1575 ï 12112 ï ï RcRR(3)=+=+= (2) (2) 0,25.1,0125 0 0,2531 îï 32332 * Xác định theo công thức (1.3.16) (3) ⎜R3 ⎜= max(0 ; 0,1575 ; 0,2531) = 0,2531 > 0,15 nhưng vì số bước lặp tối đa M = 3 nên ta dừng ở X(3) = ( 1,2 1,0125 0,85) Program PP_NOILONG; Var C: Array[1 20,1 20] of real; 17
X,R:Array[1 20] of real; I,J,K,K1,L,LMAX,M,N,DK:Integer; AB,ABMAX,DELX,EPS: Real; BEGIN Write('So an N = ');Readln(N); M:=N+1; Writeln('Nhap ma tran C '); {C là ma trận N hàng, M cột- các hệ số vế trái của hệ (1.3.4): ïì ï -+xcx ++ cxb +=' 0 ï 1122 1nn 1 ï ï cx-+ x cx ++ c x += b' 0 íï 21 1 2 23 3 2nn 2 ï ï ï ï cx++ c x - x + b' =0 îï nn11nnnn 1 1 For I:=1 to N do For J:=1 to M do Begin Write('C[',I,J,']=');Readln(C[I,J]); End; Writeln('Xap xi ban dau'); For I:=1 to N do Begin Write('X[',I,']=');Readln(X[I]); end; Write('So buoc lap toi da la LMAX=');Readln(LMAX); Write('Sai so EPS =');Readln(EPS); DK:=1;L:=0; While ((L<LMAX) and (DK=1)) do Begin For I:=1 to N do Begin R[I]:=C[I,M]; For J:=1 to N do R[I]:=R[I]+C[I,J]*X[J]; end; AB:=0; For K:=1 to N do Begin ABMAX:=ABS(R[K]); If AB<ABMAX then Begin AB:=ABMAX;K1:=K; End; end; 18
If ABS(R[K1}) K1 then R[I]:=R[I]+C[I,K1]*DELX; R[K1]:=0; end; L:=L+1; end; If DK=2 then Begin Writeln('Nghiem xap xi o buoc lap thu l =',l); For I:=1 to N do Writeln('X[',I,']=',X[I]); end; If L=LMAX then Begin Writeln('Chua hoi tu va ta tam dung o xap xi thu l cua nghiem, l=',l); For I:=1 to N do Writeln('X[',I,']=',X[I]); end; Readln; END. Bài tập 1- Giải các hệ sau đây bằng phương pháp Cholesky ⎧ x2x4x123−+=− 1 ⎧2x++= y 3z 13 ⎧ x2y3z++= 6 ⎪ ⎪ ⎪ a- ⎨−−+= 2x123 3x x 0 b-⎨ x2y3z14++= c-⎨2x+ 3y+=− z 1 d- ⎪ ⎪ ⎪ ⎩ 3x123+−= x 5x 7 ⎩3x++= 2y z 10 ⎩3x++= y 2z 7 ⎧ 2x123−+= x 3x 7 ⎪ ⎨−−3x12 2x + x 3 =− 3 ⎪ ⎩ x3x4x12+− 3 =− 4 Đs: a- ( 1 – 1 – 1 ) b- ( 1 2 3 ) c- ( 1 – 2 3 ) d- ( 1 1 2 ) {Vô số nghiệm} 2- Giải các hệ sau đây bằng phương pháp Gauss-Seidel và nới lỏng 19
⎧ 5x123 −−= x 3x 1 ⎪ a- ⎨ x123 +−= 4x 2x 3 b- ⎪ ⎩−++= 2x12 x 4x 3 3 ⎧ 4x12 −+ x 2x 3 = 5 ⎪ ⎨−+ 3x12 8x − 2x 3 = 3 ⎪ ⎩ x12 −+ 2x 10x 3 = 9 a1- Với X(0) = ( 0,8 0,9 1,1 ) , ε = 0,1 , M = 4 a2- Với X(0) = ( 0,8 0,9 1,1 ) , ε = 0,1 , M = 10 a3- Với X(0) = ( 0,8 0,9 1,1 ) , ε = 0,01 , M = 10 a4- Với X(0) = ( 1,1 0,8 0,9 ) , ε = 0,01 , M = 8 b1- Với X(0) = ( 1,1 0,8 0,9 ) , ε = 0,1 , M = 4 b2- Với X(0) = ( 1,1 0,8 0,9 ) , ε = 0,1 , M = 8 b3- Với X(0) = ( 1,1 0,8 0,9 ) , ε = 0,01 , M = 15 b4- Với X(0) = ( 0,9 1,1 0,8 ) , ε = 0,01 , M = 15 ⎧10x123−+ x 5x = 16 ⎪ ⎪ x12−−+= 10x 2x 34 6x 3 (0) c- ⎨ với X = (0,9 – 0,8 1,1 – 1,2) ; ε = 0,01 ; M = 20 ⎪−+2x12 4x + 20x 34 − 10x = 24 ⎪ ⎩−+ x12 2x − 3x3 + 25x 4 =− 31 ⎧ 10x123−−= 2x x 3 ⎪ d- ⎨−+ x123 10x − 2x = 13 e- ⎪ ⎩−− 2x123 x + 10x = 26 ⎧10x123++= 2x x 9 ⎪ ⎨ x123+−=− 10x x 22 ⎪ ⎩−+2x12 3x + 10x 3 = 22 với X(0) = (0 0 0) ; ε = 0,1 ; M = 8 với X(0) = (1,2 – 1,9 3,1) ; ε = 0,1 ; M = 6 Đs: a- Phương pháp Gauss-Seidel a1- Dừng ở xấp xỉ thứ 2 của nghiệm X(2) = (1,014 1,0015 1,006625) a2- Nghiệm xấp xỉ là X(2) = (1,014 1,0015 1,006625) a3- Nghiệm xấp xỉ là X(7) = ( 0,998953 1,000680 0,000047 ) a4- Dừng ở xấp xỉ thứ 8 X(8) = ( 1,000160 1,000646 1,000141 ) Phương pháp nới lỏng a1- Dừng ở xấp xỉ thứ 2 X(4) = ( 1,04 1,04 1,1 ). a2- Dừng ở xấp xỉ thứ 2 X(4) = ( 1,04 1,04 1,1 ). a3- Dừng ở xấp xỉ thứ 7 X(7) = ( 1,05 0,995 1,01 ). a4- Dừng ở xấp xỉ thứ 8 X(8) = ( 0,995 1,003125 0,990625 ). b- Phương pháp Gauss-Seidel b1-Dừng ở xấp xỉ thứ 4 X(4) = ( 0,995625 1,011172 0,994688 ). b2-Dừng ở xấp xỉ thứ 4 X(4) = ( 0,995625 1,011172 0,994688 ). b3- Dừng ở xấp xỉ thứ 7 X(7) = ( 1,001247 0,998936 1,000750 ). b4- Dừng ở xấp xỉ thứ 9 X(9) = ( 0,998971 1,001005 0,999428 ). Phương pháp nới lỏng 20
b1-Dừng ở xấp xỉ thứ 4 X(4) = ( 1,1 1,7475 1,1395 ). b2-Dừng ở xấp xỉ thứ 4 X(4) = ( 1,1 1,7475 1,1395 ). b3-Dừng ở xấp xỉ thứ 8 X(8) = ( 1,132094 1,830643 1,1152919 ). b4-Dừng ở xấp xỉ thứ 11 X(11) = ( 1,130200 1,884884 1,162042 ). c - Phương pháp Gauss-Seidel Dừng ở xấp xỉ X(6) = ( 0,999477 –1,000426 0,999579 – 0,999930 ). Phương pháp nới lỏng Dừng ở xấp xỉ X(18) = ( 1,41845 0,989316 0,543982 – 1,2 ) Kiểm tra nhận thức 1* Trình bày cụ thể các bước tìm nghiệm bằng phương pháp Cholesky với n = 4 và giải hệ ⎧ x2x3x4x1234+++=− 2 ⎪ ⎪2x1234+++= 3x 4x x 2 ⎨ ⎪3x1234+++=− 4x x 2x 2 ⎪ ⎩4x1234+++= x 2x 3x 2 2* Nêu những điểm giống nhau, khác nhau của hai phương pháp Gauss-Seidel và nới lỏng. 3* Hãy chạy chương trình PASCAL đã có để giải các ví dụ. Nếu có điều gì chưa thoả mãn thì chủ động viết lại chương trình (có thể bằng ngôn ngữ lập trình bậc cao khác) và giải lại các ví dụ, các bài tập. 21
Chương 2 : PHƯƠNG TRÌNH VÀ HỆ PHƯƠNG TRÌNH PHI TUYẾN Trong thực tế ta gặp không ít các bài toán gắn liền với phương trình đại số cấp cao hoặc siêu việt và hệ phương trình phi tuyến. Ta biết rằng cách xấp xỉ hệ bởi hệ tuyến tính dẫn đến sai số rất lớn. Vì vậy trong chương này ta sẽ trình bày một số phương pháp trực tiếp xấp xỉ nghiệm. 2.1 Phương pháp chia đôi Giải gần đúng phương trình F(x) = 0 (2.1.1) thoả mãn điều kiện * F '(x) không đổi dấu trên [a ; b] (2.1.2) * F(a)F(b) 0, cho trước (sai số đối với hàm) thì ta dừng ở x(1). * Ngược lại, có nghĩa là ⎜F(x(1)) ⎜ > εf (2.1.6) thì ta tìm xấp xỉ thứ hai của nghiệm là x(2) bởi công thức (1) (1) (2) ab+ x = 2 (2.1.7) trong đó a(1) = a(0) , b(1) = x(1) nếu F(a(0))F(x(1)) < 0 hoặc a(1) = x(1) , b(1) = b(0) nếu F(b(0))F(x(1)) < 0 * Nếu 22
⎜x(2) – x(1) ⎜ ≤ εx (2.1.8) trong đó εx > 0 cho trước (sai số đối với nghiệm) thì ta dừng ở x(2). * Ngược lại, có nghĩa là ⎜x(2) – x(1) ⎜ > εx (2.1.9) thì ta tính F(x(2)). * Nếu ⎜F(x(2)) ⎜ ≤ εf (2.1.10) thì ta dừng ở x(2). * Ngược lại, có nghĩa là ⎜F(x(2)) ⎜ > εf (2.1.11) ta tính tương tự x(3) . . . * Một cách tổng quát Ta dừng ở xấp xỉ thứ k của nghiệm là (1)(1)kk ()k ab+ x = 2 (2.1.12) trong đó a(k-1) = a(k-2) , b(k-1) = x(k-1) nếu F(a(k-2))F(x(k-1)) xxe ï ï í Fx(1)k- > e f ï () ï ï kM= îï (2.1.15) trong đó M là số bước lặp tố đa cho trước mà ta phải dừng lại, mặc dù đến bước lặp thứ M – 1 vẫn chưa đảm bảo sai số đối với nghiệm và sai số đối với hàm. Chú ý Điều kiện (2.1.13) được thoả mãn khi 23
ba- ln k ³ ex ln 2 (2.1.16) Ví dụ 2.1 Xấp xỉ bằng phương pháp chia đôi nghiệm của phương trình sinx – 3x + 2 = 0 với [a ; b] = [0 ; 1] , εx = 0,01 , εf = 0,01 và M = 10 Giải * F(x) = sinx – 3x + 2 nên ta có F '(x) = cosx – 3 0 , F(1) = – 0,158529 0 * ⎜F(x(1)) ⎜ = 0,979426 > 0,01 * Vì F(1) F(x(1)) 0,01 * F(x(2)) = sin0,75 – 3.0,75 + 2 » 0,431639 > 0 . * ⎜F(x(2) ⎜ » 0,431639 > 0,01 . * Vì F(1) F(x(2)) 0,01 * F(x(3)) = sin0,875 – 3.0,875 + 2 » 0,142544 > 0 * ⎜F(x(3)) ⎜» 0,142544 > 0,01 * Vì F(1) F(x(3)) 0,01 * F(x(4)) = sin0,9375 – 3.0,9375 + 2 » – 0,006419 < 0 * ⎜F(x(4)) ⎜ » 0,006419 < 0,01 Vậy ta dừng ở xấp xỉ thứ tư của nghiệm x(4) = 0,9375 Program CHIADOI; Var 24
K,DK,M:Integer; A,B,C,YA,YB,YC,EPSX:Real; Function F(X:Real):Real; Begin F:=sin(X) – 3*X + 2; End; {Giải phương trình sinx – 3x + 2 = 0 trên [A;B] } BEGIN Write('A= ');Readln(A); Write('B = ');Readln(B); Write(' Sai so theo X : EPSX = ');Readln(EPSX); M:=1+Trunc((ln(B- A)- ln(EPSX))/ln(2)); Writeln('So buoc lap toi thieu la',M); YA:=F(A);YB:=F(B);DK:=1; If YA*YB > 0 Then DK:=2;K:=1; While ((K 0 Then Begin B:=C;YB:=YC; End Else Begin A:=C;YA:=YC; End; K:=K+1; End; If DK=3 Then Writeln('Nghiem chinh xac la ',C); If DK=1 Then Writeln('Nghiem xap xi la X(',K,')=',C); If DK=2 Then Writeln('Chon A va B sai.Phai chon lai A,B'); Readln; END. 2.2 Phương pháp dây cung Giải gần đúng phương trình F(x) = 0 (2.2.1) thoả mãn điều kiện 25
* F '(x) không đổi dấu trên [a,b] (2.2.2) * F(a)F(b) 0 , cho trước (sai số đối với hàm) thì ta dừng ở x(1). * Ngược lại, có nghĩa là ⎜F(x(1)) ⎜ > εf (2.2.6) thì tìm xấp xỉ thứ hai của nghiệm là x(2) bởi công thức Fb(1) b (1)- a (1) (2) (1) ()( ) xb=- Fb()(1)- Fa () (1) (2.2.7) trong đó a(1) = a(0) , b(1) = x(1) nếu F(a(0))F(x(1)) 0 , cho trước (sai số đối với nghiệm) thì ta dừng ở x(2). * Ngược lại, có nghĩa là ⎜x(2) – x(1) ⎜> εx (2.2.9) ta tính F(x(2)) . . . * Một cách tổng quát :Ta lặp cho đến xấp xỉ thứ k của nghiệm là Fb(1)(1)(1)kk b- a k - ()kk (- 1) ()( ) xb=- Fb()()(1)kk Fa (1) (2.2.10) trong đó 26
a(k - 1) = a(k - 2) , b(k - 1) = x(k - 1) nếu F(a(k - 2))F(x(k - 1)) xxe ï ï í Fx(1)k- > e f ï () ï ï kM= îï (2.2.13) 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ù đến bước lặp thứ M – 1 vẫn chưa đảm bảo sai số đối với nghiệm và sai số đối với hàm. Ví dụ 2.2 Xấp xỉ nghiệm của phương trình sau đây bằng phương pháp dây cung : sinx – 3x + 2 = 0 với [a ; b] = [0 ; 1] , εx = 0,01 , εf = 0,01 và M = 10 Giải * F(x) = sinx – 3x + 2 nên ta có F '(x) = cosx – 3 0 , F(1) » – 0,158529 0 * ⎜F(x(1)) ⎜ = 0,019886 > 0,01 . * Vì F(1) F(x(1)) < 0 nên ta tính (2) 0,158529() 1 0,926557 x »-1 » 0, 934774 0,158529 0,019886 Dấu của F(x) + + + – * * * * 0 x(1) x(2) 1 x * ⎜x(2) – x(1) ⎜ = ⎜ 0,934774 – 0,926557 ⎜ = 0,008217 < 0,01 Vậy ta dừng ở xấp xỉ thứ hai của nghiệm x(2) » 0,934774 27
Program DAYCUNG; Var K,M,DK:Integer; A,B,C,YA,YB,YC,EPSF:Real; { Giải phương trình sinx – 3x + 2 = 0 trên [A;B] } Function F(X:Real):Real; Begin F:=sin(X) – 3*X + 2; end; BEGIN Write('A=');Readln(A); Write('B=');Readln(B); Write(' Sai so theo ham : EPSF =');Readln(EPSF); Write(' So buoc lap toi da M =');Readln(M); YA:=F(A);YB:=F(B);K:=1;DK:=1; While ((K 0 then Begin B:=C;YB:=YC; end else Begin A:=C;YA:=YC; end; end; If ABS(YC) < EPSF then DK:=3;K:=K+1; end; If DK = 3 then Writeln('Nghiem xap xi la X(',K–1,') = ',C); If DK = 2 then Writeln('Nghiem chinh xác la',C); If DK = 1 then Writeln('Chua dam bao sai so va ta dung o X',K,')=',C); Readln; END . 2.3 Phương pháp tiếp tuyến Giải gần đúng phương trình F(x) = 0 (2.3.1) 28
thoả mãn điều kiện * F '(x), F ”(x) không đổi dấu trên [a,b] (2.3.2) * F(a)F(b) 0) . * Ta tính xỉ thứ nhất của nghiệm theo công thức y Fx(0) (1) (0) () xx=- (2.3.4) y = F(x) Fx'()(0) x trong đó 0 a x x* b x(0) = a nếu F(a)F ”(a) > 0 hoặc x(0) = b nếu F(b)F ”(b) > 0 * Nếu ⎜x(1) – x(0) ⎜ ≤ εx (2.3.5) 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à ⎜x(1) – x(0) ⎜ > εx (2.3.6) ta tính F(x(1)). * Nếu ⎜F(x(1)) ⎜ ≤ εf (2.3.7) trong đó εf > 0 cho trước (sai số đối với hàm) thì ta dừng ở x(1). * Ngược lại, có nghĩa là ⎜F(x(1)) ⎜ > εf (2.3.8) ta tính xấp xỉ thứ 2 theo công thức Fx(1) (2) (1) () xx=- Fx'()(1) (2.3.9) . . . * Một cách tổng quát Ta tính cho đến xấp xỉ thứ k của nghiệm là 29
Fx(1)k- ()kk (- 1) () xx=- Fx'()(1)k- (2.3.10) và dừng lại ở x(k) khi thoả mãn một trong các điều kiện ⎜x(k) – x(k-1) ⎜ ≤ εx (2.3.11) hoặc ⎜ F(x(k)) ⎜ ≤ εf (2.3.12) hoặc ì (1)(2)kk ï x ->xxe ï ï í Fx(1)k- > e f ï () ï ï kM= îï (2.3.13) 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ù đến bước lặp thứ M – 1 vẫn chưa đảm bảo sai số đối với nghiệm và sai số đối với hàm. Ví dụ 2.3 Xấp xỉ nghiệm của phương trình sau đây bằng phương pháp tiếp tuyến sinx – 3x + 2 = 0 với [a ; b] = [0 ; 1] , εx = 0,01 , εf = 0,001 và M = 10 Giải * F(x) = sinx – 3x + 2 nên ta có F '(x) = cosx – 3 0 . F(1) = sin1 – 3.1 + 2 » – 0,158529 0 (0) ()10 () F()x sin1-+ 3.1 2 * xx=- = 1-»0,935549 Fx'()(0) cos1- 3 * ⎜x(1) – x(0) ⎜ = ⎜0,935549 – 1 ⎜ = 0,064451 > 0,01 * F(x(1)) = sin0,935549 – 3. 0,935549 + 2 » – 0,001722 * ⎜F(x(1)) ⎜= 0,001722 > 0,001 Dấu của F(x) + – * * * * Dấu của F”(x) 0 – – x(1) – x(2) 1 x 30
(1) ()21 () F()x sin0,935549-+ 3.0,935549 2 * xx=- »-0,935549 » 0,936265 Fx'()(1) cos0,935549- 3 ()21 () * xx-=0,936265 - 0,935549 = 0,000716 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à ⎜x(1) – x(0) ⎜ > εx (2.4.8) ta tính F(x(1)). * Nếu ⎜F(x(1)) ⎜ ≤ εf (2.4.9) trong đó εf > 0, cho trước (sai số đối với hàm) thì ta dừng ở x(1). * Ngược lại, có nghĩa là ⎜F(x(1) ⎜ > εf (2.4.10) 31
ta tính xấp xỉ thứ hai của nghiệm là x(2) = ϕ (x(1)) . . . * Một cách tổng quát Ta tính cho đến xấp xỉ thứ k của nghiệm là x(k) = ϕ (x(k - 1)) và dừng ở x(k) khi thoả mãn một trong các điều kiện ⎜x(k) – x(k-1) ⎜ ≤ εx (2.4.11) hoặc ⎜F(x(1) ⎜ ≤ εf (2.4.12) hoặc ì (1)(2)kk ï x ->xxe ï ï í Fx(1)k- > e f ï () ï ï kM= îï (2.4.13) 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ù đến bước lặp thứ M – 1 vẫn chưa đảm bảo sai số đối với nghiệm và sai số đối với hàm. Chú ý Lúc này phương pháp lặp hội tụ và hơn nữa với ký hiệu nghiệm chính xác là x* thì xx(1)- (0) k xxq()k -£* 1 - q (2.4.14) Ví dụ 2.4 Xấp xỉ nghiệm của phương trình sau đây bằng phương pháp lặp đơn x3 – x – 1 = 0 với x ( 0) = 1 , εx = 0,05 , εf = 0,05 , [a ; b] = [1 ; 2] và M = 10 Giải * F(x) = x3 – x – 1 nên ta có F '(x) = 3x2 – 1 > 0 , ∀x ∈ [1 ; 2]. * F(1) = – 1 0 cho nên F(1)F(2) < 0 . * Cách biến đổi x = x3 – 1 = ϕ(x) không thoả mãn (2.4.6) vì ϕ (2) = 7 ∉ [1 ; 2]. * Cách biến đổi 11 x =+= ϕ(x) x 2 x 11 3 không thoả mãn (2.4.6) vì ϕ(2) = + = ∉ [1 ; 2]. 2 2 2 4 * Cách biến đổi xx=+3 1 = ϕ(x) 32
thoả mãn (2.4.5) và (2.4.6) vì + ϕ’(x) = 1 là hàm giảm 313 ()x+ 2 Giá trị lớn nhất của ϕ’(x) là ϕ’(1) = 1 » 0,209987 3113 ()+ 2 Giá trị nhỏ nhất của ϕ’(x) là ϕ’(2) = 1 » 0,160250 3213 ()+ 2 cho nên ⎜ϕ’(x) ⎜ 0 nên ϕ(x) là hàm tăng Giá trị lớn nhất của ϕ(x) là ϕ (2) = 3 21+ » 1,442250 Giá trị nhỏ nhất của ϕ(x) là ϕ (1) = 3 11+ » 1,259921 cho nên ϕ(x) ∈ [1 ; 2] , ∀x ∈ [1 ; 2] (1) (0) * x = ϕ (x ) = 3 11+ » 1,259921 . * ⎜x(1) – x(0) ⎜ = ⎜1,259921 – 1 ⎜ = 0,259921 > 0,05 . * F(x(1)) = 1,2599213 – 1,259921 – 1 » – 0,259921 * ⎜F(x(1)) ⎜ » 0,259921 > 0,05 . (2) (1) * x = ϕ (x ) = 3 1,259521+ 1 » 1,312294 * ⎜x(2) – x(01 ⎜ = ⎜1,312294 – 1,259921 ⎜ = 0,052373 > 0,05 * F(x(2)) = 1,3122943 – 1,312294 – 1 » – 0,052372 * ⎜F(x(2)) ⎜ » 0,052372 > 0,05 (3) (2) * x = ϕ (x ) = 3 1,312294+ 1 » 1,322354 * ⎜x(3) – x(2) ⎜ = ⎜1,322354 – 1,312294 ⎜ = 0,01006 < 0,05 Vậy ta dừng ở xấp xỉ thứ ba của nghiệm x(3) » 1,322354 2.5 Phương pháp Newton-Raphson cho hệ phương trình Giải hệ phương trình phi tuyến ïì fx( , , x )= 0 ï 11 n íï ï îï fxnn(1 , , x )= 0 (2.5.1) trong đó các hàm fi(x1, . . . , xn), i =1, . . . , n có các đạo hàm riêng cấp 2 liên tục và bị chặn. Nội dung của phương pháp 33
Xấp xỉ hệ đã cho bởi hệ dừng ở đạo hàm riêng cấp 1 theo khai triển Taylor của hàm nhiều biến n (0) (0) (0) (0) ¶ fxin(1 , , x ) (0) f i (x1, . . . , xn) » f i (x1 , . . . , xn ) + å ()xxj - j ; i = 1 , . . . , n j = 1 ¶ x j (2.5.2) (0) (0) (0) trong đó X = (x . . . xn ) là xấp xỉ ban đầu của nghiệm được cho trước. * Ký hiệu với X = (x1 . . . xn) é ù ê f xx, , ú ê 11()n ú ê ú F(X) = ê ú ê ú êf xx, , ú ëê nn()1 ûú (2.5.3) é ù ê¶¶ff11()xx11, , nn() xx, , ú ê ú ê ¶¶xx1 n ú ê ú J(X) = ê ú ê ú ê¶¶ffnx(), , x ()x, , x ú ê 11nn n ú êë ¶¶xx1 n ûú (2.5.4) (0) (0) (0) * Tính tại xấp xỉ ban đầu X = (x1 . . . xn ) của nghiệm éù(0) êúfx, , x(0) êú1()1 n (0) êú F(X ) = êú êú êúfx(0) , , x(0) ëûêúnn()1 (2.5.5) é ()00()00() () ù ê¶¶ff11()xx11, , nn() xx, , ú ê ú ê ¶¶xx1 n ú (0) ê ú J(X ) = ê ú ê ú ê nx()00, , x()00x(), , x() ú ê¶¶ff()11nnn()ú ê ú êë ¶¶xx1 n ûú (2.5.6) * Tính ma trận nghịch đảo (J(X(0))) - 1 (1) (1) (1) * Tìm xấp xỉ thứ nhất X = (x1 . . . xn ) của nghiệm theo công thc X(1) = X(0) – (J(X(0)))- 1F(X(0)) (2.5.7) * Nếu n (1) (0) å xxjj- ≤ εx j= 1 (2.5.8) 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à n (1) (0) å xxjj- > εx j= 1 (2.5.9) ta tính 34
X(2) = X(1) – (J(X(1)))- 1F(X(1)) (2.5.10) . . . * Một cách tổng quát Ta dừng lại ở xấp xỉ thứ k của nghiệm X (k) = X (k - 1) – (J(X (k - 1))) - 1F(X (k - 1)) (2.5.11) thoả mãn một trong các điều kiện sau đây n ()kk (- 1) å xxjj- ≤ εx j= 1 (2.5.12) hoặc ì n ï (1)(2)kk ï å xii->xxe íï ï i= 1 ï îï kM= (2.5.13) 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ù đến bước lặp thứ M – 1 vẫn chưa đảm bảo sai số đối với nghiệm. Ví dụ 2.5 Giải hệ phương trình sau đây bằng phương pháp Newton-Raphson ì 32 2 ï xx+-10 x - 2 x += 1 0 (0) íï 11 1 2 với X = (0 0) , εx ï 22 îï 232040xx12+- x 2 -= = 0,1, M=15 Giải é ù ê fxx, ú é 32 2ù ê 112()ú êxx11+-10 x 1 - 2 x 2 + 1ú * F(X) = ê ú = ê ú êf xx, ú ê 23204xx22+- x -ú ëê 21() 2ûú ë 12 2û é ù ê¶¶ff11()xx12,,() x1 x 2ú é 2 ù ê ¶¶xx2 ú ê32104xx+- - xú J(X) = ê 1 ú = ê 11 2ú ê ú ê ú ê¶¶ff2,()x12xxx2()12, ú ê ú ëê 462xx12- 0ûú ëê ¶¶xx1 2 ûú é (0) (0) ù êfx1 , x ú é 32 2ù é ù (0) ê ()12ú ê0+- 0 10.0 - 2.0 + 1ú ê1 ú * F(X ) = ê ú = ê ú = ê ú êfx(0), x (0) ú ê 2.022+ 3.0 20.0 4 ú ë- 4û ëê 2 ()12ûú ë û - 1 éù2 éù- 1 (0) - 1 êú3.0+- 2.0 10 - 4.0 êú- 10 0 * (J(X )) = êú= êú ëûêú4.0 6.0- 20 ëû020- éùéù 1 êú- 20 0 êú- 0,1 0 = êú = êú ( 10).( 20) 0.0 ëû010- ëû00,05- Chú ý éù- 1 é ù êúab 1 êd- bú êú= ê ú ëûcd ad - bc ë- c aû 35
éù é ù é ù é ù (1) (0) (0) - 1 (0) êú0 ê- 0,1 0 ú ê1 ú ê 0,1 ú * X = X – (J(X )) .F(X ) = êú – ê ú. ê ú= ê ú ëû0 ë 00,05- û ë- 4û ë- 0,2û (1) (0) (1) (0) * ⎜x1 – x1 ⎜ + ⎜x2 – x2 ⎜ = ⎜0,1 – 0 ⎜ + ⎜– 0,2 – 0 ⎜ = 0,3 > 0,1 é (1) (1) ù êfx1 , x ú é 32 2ù é ù (1) ê ()12ú ê0,1+- 0,1 10.0,1 + 2.( 0,2) 1ú ê- 0,069ú * F(X ) = ê ú = ê ú = ê ú êfx(1), x (1) ú ê2.0,122+- 3.( 0,2) - 20.( - 0,2) - 4ú ë 0,14 û ëê n ()12ûú ë û - 1 éù2 (1) - 1 êú3.0,1+- 2.0,1 10 4.( 0,2) * (J(X )) = êú ëûêú4.0,1 6.( 0,2) 20 éù- 1 éù êú- 9,77 0,8 1 êú 21,2 0,8 = êú = êú ëû0,4- 21,2 ( 9,77).( 21,2) 0,8.0,4 ëû 0,4 9,77 éù - 1 êú21,2 0,8 = êú 206,804 ëû0,4 9,77 é ù éùé ù é ù (2) (1) (1) - 1 (1) ê 0,1 ú - 1 êú21,2 0,8 ê- 0,069ú ê 0,093468 ú * X = X – (J(X )) .F(X ) = ê ú – êúê ú » ê ú ë- 0,2û 206,804 ëû0,4 9,77 ë 0,14 û ë- 0,193519û (2) (1) (2) (1) * ⎜x1 – x1 ⎜ + ⎜x2 – x2 ⎜ = ⎜0,093468 – 0,1⎜+ ⎜– 0,193519 – (– 0,2) ⎜ = 0,013013 < 0,1 Vì vậy ta dừng ở xấp xỉ thứ hai của nghiệm X(2) = ( 0,093468 – 0,193519 ) Program NEWRAPHSON_HPT; Var DK,K,K1,M: Integer; X10,X20,X1,X2,D11,D12,D21,D22,DX1,DX2, F10,F20,F11,F22,EPS,DET: Real; ⎧ 32 2 ⎪x11 +−x 10x 1 − 2x 2 += 1 0 { Giải hệ ⎨ } 22 ⎩⎪2x12 +− 3x 20x 2 −= 4 0 Function F1(X1,X2:Real):Real; Begin F1:= X1*X1*X1 + X1*X1 – 10 *X1 – 2*X2*X2 + 1; end; Function F2(X1,X2:Real):Real; Begin F2:= 2*X1*X1 + 3*X2*X2 – 20*X2 – 4 ; end; BEGIN Write('Xap xi ban dau X(0) '); (0) (0) Write('X1 =');Readln(X10);Write('X2 =');Readln(X20); Write(' ε =');Readln(EPS); Write('So buoc lap toi da la');Readln(M); DK:=1;K:=1; F10:=F1(X10,X20); F20:=F2(X10,X20); While (K<M) and (DK=1) do Begin D11:= 3*X10*X10 + 2*X10 – 10 ; D12:= – 4*X20 ; D21:= 4*X10*X10 ; D22:= 6*X20 – 20 ; DET:=D11*D22–D12*D21; 36
{ D11 là đạo hàm riêng của F1 theo X1 tại (X10,X20), D12 là đạo hàm riêng của F1 theo X2 tại (X10,X20), D21 là đạo hàm riêng của F2 theo X1 tại (X10,X20), D22 là đạo hàm riêng của F2 theo X2 tại (X10,X20). } If DET=0 then Begin DX1:=0;DX2:=0; End ELSE Begin DX1:=(F10*D22–F20*D12)/DET; DX2:=(F20*D11–F10*D12)/DET; End; X1:=X10–DX1;X2:=X20-DX2; F11:=F1(x1,X2);F22:=F2(x1,X2); If (ABS(F11)+ABS(F22) < EPS) then DK:=3 ELSE Begin X10:=X1;X20:=X2;F10:=F11;F20:=F22; End; K1:=K;K:=K+1; End; If K1 < M then Begin Writeln('Xap xi cua nghiem la '): Writeln('X1[',K1,']=',X1,'X2[',K1,']=',X2); End ELSE If DK=1 then Begin Writeln('Chua hoi tu va ta dung o xap xỉ'); Writeln('X1[',K1,']=',X1,'X2[',K1,']=',X2); End; Readln; END. 2.6 Phương pháp lặp Seidel cho hệ phương trình Giải hệ phương trình phi tuyến ïì fx11( , , xn )= 0 ï íï ï îï fxnn(1 , , x )= 0 (2.6.1) trong đó các hàm fi(x1, . . . , xn), i =1, . . . , n có các đạo hàm riêng cấp 2 liên tục và bị chặn. Nội dung của phương pháp 37
Giải hệ sau đây tương đương vơi hệ (2.6.1) : ïì x = gx( , , x ) ï 111 n íï ï îï x nn= gx(1 , , x n ) (2.6.2) thoả mãn điều kiện n ¶ g å i 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à n (1) (0) å xxjj- > εx j= 1 (2.6.6) (1) (1) (1) thay xấp xỉ thứ nhất X = (x1 . . . xn ) của nghiệm vào vế phải của (2.6.2) (2) (2) (2) để tìm xấp xỉ thứ hai của nghiệm là X = (x1 . . . xn ) ïì (2) (1) (1) ï xgxx11= 1(), , n ï ï (2) (2) (1) (1) ï xgxx= , , , x íï 2122()n ï ï ï ï xgxx(2) = (2), (2) , , xx (2) , (1) îï nn()12n- 1n (2.6.7) . . . * Một cách tổng quát (2) (2) (2) Ta dừng lại ở xấp xỉ thứ k của nghiệm là X = (x1 . . . xn ) theo công thức 38
ïì ()kk (- 1) (1)k- ï xgxx11= 1(), , n ï ï ()kkk () (- 1) (1)k- ï xgxx= , , , x íï 2122()n ï ï ï ï xgxx()kk= ()kk, () , , xx () k , (- 1) îï nn()12n- 1n (2.6.8) thoả mãn một trong các điều kiện sau đây n ()kk (- 1) å xxjj- ≤ εx j= 1 (2.6.9) hoặc ì n ï (1)(2)kk ï å xii->xxe íï ï i= 1 ï îï kM= (2.6.10) 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ù đến bước lặp thứ M – 1 vẫn chưa đảm bảo sai số đối với nghiệm. Ví dụ 2.6 Giải hệ phương trình sau đây bằng phương pháp lặp Seidel ì 32 2 ï xx+-10 x - 2 x += 1 0 (0) íï 11 1 2 với X = (0 0) , εx = 0,3, M=5 trong miền 0 ≤ x ≤ 1 ï 22 1 îï 232040xx12+- x 2 -= , – 1 ≤ x2 ≤ 1 Giải * Đưa hệ đã cho về hệ tương đương ïì 32 2 ï x 111=+-0,1xx 0,1 0,2 x 2 += 0,1 gxx 112() , íï ï 22 ï x = 0,1xx+- 0,15 0,2 = gxx() , îï 21 2 212 * Kiểm tra điều kiện (2.6.3) với 0 ≤ x1 ≤ 1 , – 1 ≤ x2 ≤ 1 ¶¶gg11 2 +=0,3x11 + 0,1xx +- 0,4 2 ¶¶xx12 2 2 ≤ 0,3⎜x1⎜ + 0,1⎜x1⎜ + 0,4⎜x2⎜ ≤ 0,3.1 + 0,1.1 + 0,4.1 = 0,8 < 1 ¶¶gg22 +=0,2x 12 + 0,3x = 0,2⎜x1⎜ + 0,3⎜x2⎜ ≤ 0,2.1 + 0,3.1 = ¶¶xx12 0,5 < 1 (0) (1) (1) (1) * Thay xấp xỉ ban đầu của nghiệm là X để tìm xấp xỉ thứ nhất X = (x1 x2 ) theo công thức (2.6.4) ïì (1) 32 2 ï x 1 =+-+=0,1.0 0,1.0 0,2.0 0,1 0,1 íï ï x (1) =+-=-0,1.0,122 0,15.0 0,2 0,199 îï 2 * Tính sai số theo công thức (2.6.5) 39
(1) (0) (1) (0) ⎜x1 – x1 ⎜ + ⎜x2 – x2 ⎜ = ⎜0,1 – 0 ⎜ + ⎜ – 0,199 – 0 ⎜ = 0,299 > 0,03 (1) (2) (2) (2) * Thay xấp xỉ thứ nhất của nghiệm là X để tìm xấp xỉ thứ hai X = (x1 x2 ) theo công thức (2.6.7) ïì (2) 322 ï x 1 =+ +»0,1.0,1 0,1.0,1 0,2.( 0,199) 0,1 0,09318 íï ï x (2) =+ »-0,1.0,0931822 0,15.( 0,199) 0,2 0,193192 îï 2 * Tính sai số theo công thức (2.6.9) (2) (1) (2) (1) ⎜x1 – x1 ⎜ + ⎜x2 – x2 ⎜ = ⎜0,09318 – 0,1⎜ + ⎜– 0,193192 – (– 0,199) ⎜ = 0,012628 EPS ,EPSX=');Readln(EPSX); Write('So buoc lap toi da la M = ');Readln(M); K:=1; While ((K EPS)) do Begin X1[k]:=G1(X1[K-1],X2[K–1]); X2[K]:=G2(X1[K],X2[K–1]); EPSX:=ABS(X1[K] –X1[K–1])+ABS(X2[K] –X2[K–1]); K1:=K;K:=K+1; end; If K1<M Then Begin Writeln('Xấp xi cua nghiem la X1[',K1,']=',X1[K1],'X2[ ',K1,']=',X2[K1]); 40
end ELSE Begin Writeln(' Chua hoi tu va ta dung o xap xi '); Writeln(' X1[',K1,']=',X1[K1],'X2[',K1,']=',X2[K1]); end; Readln; END. Bài tập 1- Giải phương trình ex – x – 2 = 0 trên [a;b] = [1;2] với εx = εf = 0,01 , M = 10 và với εx = εf = 0,001 , M = 20 . a- Phương pháp chia đôi b- Phương pháp dây cung c- Phương pháp tiếp tuyến d- Phương pháp lặp đơn khi x(0) = 1,5 và khi x(0) = 1,3 Đs: a- x(7) = 1,140625 x(10) = 1,146844 b- x(5) = 1,143013 x(8) = 1,145885 c- x(4) = 1,147913 x(5) = 1,146152 d- x = ln(x+2): Nghiệm xấp xỉ khi x(0) = 1,5 là x(4) = 1,149535 Nghiệm xấp xỉ khi x(0) = 1,3 là x(5) = 1,146346 2- Giải phương trình xcosx – 3x +1 = 0 trên [a;b] = [0;1] với εx = εf = 0,01 , M = 12 và với εx = εf = 0,001 , M = 20. a- Phương pháp chia đôi b- Phương pháp dây cung c- Phương pháp tiếp tuyến d- Phương pháp lặp đơn khi x(0) = 0,5 và khi x(0) = 0,2 Đs: a- x(7) = 0,484375 x(10) = 0,474609 b- x(3) = 0,472150 x(4) = 0,473608 c- x(4) = 0,474091 x(4) = 0,474091 41
x cosx+ 1 d- x = : 3 Nghiệm xấp xỉ khi x(0) = 0,5 là x(2) = 0,475163 Nghiệm xấp xỉ khi x(0) = 0,2 là x(5) = 0,473679 3- Giải hệ phương trình sau bằng phương pháp Seidel, Newton-Raphson ⎧ 3 ⎪ −+ x12 x = 0 ⎨ , 0 ≤ x ≤ 2 , 0,5 ≤ x ≤ 3 22 1 2 ⎩⎪−−+= 4x12 9x 36 0 a- X(0) = ( 1,2 0,6 ) , εx = 0,0001 , M = 10 b- X(0) = ( 1,3 0,7 ) , εx = 0,001 , M = 12 Đs: Phương pháp Seidel (7) x1 = 1,229450 (7) x2 = 1,858395 (5) x1 = 1,229315 (5) x2 = 1,858391 Phương pháp Newton-Raphson (5) x1 = 1,222370 (5) x2 = 1,826449 (4) x1 = 1,222372 (4) x2 = 1,826451 4- Giải hệ phương trình sau bằng phương pháp Seidel, Newton-Raphson ⎧ x ⎪ sin(x121 +− 2x ) x e2 = 0 ⎨ 1 ≤ x ≤ 2 , 5 ≤ x ≤ 7 32 1 2 ⎪⎩(x12 + x ) −+= 9cosx 2 3 0 a- X(0) = ( 1,2 6 ) , ε = 0,01 , M = 5 b- X(0) = ( 1,3 6,5 ) , ε = 0,01 , M = 6 5- Giải hệ phương trình sau bằng phương pháp Seidel, Newton-Raphson (Viết lại chương trình tương ứng cho n = 3) 22 ⎧x −++−= x x x 5 0 ⎪ 1123 ⎪ 22 2 ⎨x1223 +−+−= x x x 4 0 ⎪ 222 ⎩⎪x12 ++−−= x x3 x 3 6 0 a- X(0) = ( – 0,7 2,1 1,9 ) , ε = 0,01 , M = 4 b- X(0) = (– 0,8 2,2 1,8 ) , ε = 0,001 , M = 5 Kiểm tra nhận thức 1* Nêu ưu điểm, nhược điểm của từng phương pháp chia đôi, dây cung, tiếp tuyến, lặp đơn. 2* Viết chương trình PASCAL cho phương pháp tiếp tuyến, lặp đơn. 3* Nêu ưu điểm, nhược điểm của từng phương pháp Seidel, Newton-Raphson. 42
4* Bạn hãy chạy chương trình PASCAL đã có để giải các ví dụ. Nếu có điều gì chưa thoả mãn thì chủ động viết lại chương trình (có thể bằng ngôn ngữ lập trình bậc cao khác) và giải lại các ví dụ, các bài tập. 43
Chương 3 : NỘI SUY GIÁ TRỊ HÀM SỐ Ta thường gặp bài toán sau đây: Với các mốc quan sát xk , k = – M, – M + 1, . . . , 0, 1, . . . , N – 1, N và các giá trị quan sát tương ứng f k = f (x k), cần nội suy giá trị của f (x) , x- M < x < xN. Trong chương này ta sẽ xét một số phương pháp giải quyết bài toán nêu ra. Nội dung tư tưởng của bài toán nội suy: Tìm đa thức P(x) sao cho tại các mốc quan sát xk : f(xk) = P(xk). Như vậy ta tính được giá trị gần đúng của f(x): f(x) ≈ P(x). 3.1 Công thức nội suy Gregory-Newton tiến * Các mốc quan sát xk cách đều với h là khoảng cách giữa 2 mốc quan sát liên tiếp. * Các giá trị quan sát fk = f(xk) . * Tính f(x) khi x ở khoảng đầu các mốc quan sát. * Chọn x0 là một trong các mốc đầu quan sát và gần x nhất. * Tính u = x- x0 . h * Lập bảng số gia hữu hạn tiến 2 3 4 xk = x0 + kh fk Δfk Δ fk Δ fk Δ fk . . . x - 1 f - 1 Δf - 1 2 x 0 f 0 Δ f - 1 3 Δf0 Δ f - 1 2 4 x1 f 1 Δ f0 Δ f - 1 3 Δf1 Δ f0 2 x2 f2 Δ f1 Δf2 x3 f3 . . . . . . . . . . . . trong đó : Số gia hữu hạn tiến bậc một là Δfk = fk + 1 – fk 2 Số gia hữu hạn tiến bậc hai là Δ fk = Δfk + 1 – Δfk 3 2 2 Số gia hữu hạn tiến bậc ba là Δ fk = Δ fk + 1 – Δ fk . . . . . . Chú ý * Để thuận tiện cho cho việc tra cứu các số liệu, ta đóng khung các giá trị và lấy đó làm chuẩn. * Giá trị cột thứ 3 trở đi được tính bằng cách lấy giá trị dưới trừ giá trị trên tương ứng cùng cột trước. Nguyên tắc này còn áp dụng cho các bảng số gia hữu hạn về sau. 44
DDff23 D f DN f * f() x»+ f00 u + uu( -+ 1) 0 uu ( - 1)( u - 2) ++ 0uu ( - 1)( u - 2) ( u - N + 1) (3.1.1) 0 1! 2! 3! N! Vế phải là đa thức bậc N. * Sai số là f()(N + 1) x N + 1 E =- h u(u 1) . . . (u - N + 1) ; x < ξ < x N (N+ 1)! 0 N (3.1.2) Ví dụ 3.1 Tính gần đúng giá trị f( – 3) từ bảng số liệu xk – 4 – 2 0 2 4 6 fk – 64 – 8 0 8 64 216 Giải * Ở đây ta thấy rằng xk là các mốc cách đều với h = 2 và x = – 3 ở khoảng đầu các mốc quan sát xk nên ta chọn x0 = – 4. * u = x- x0 = 3(4)= 0,5. h 2 * Bảng số gia hữu hạn tiến tương ứng 2 3 4 5 xk = – 4 + 2k fk Δfk Δ fk Δ fk Δ fk Δ fk x 0 = – 4 – 64 56 x 1 = – 2 – 8 – 48 8 48 x 2 = 0 0 0 0 8 48 0 x 3 = 2 8 48 0 56 48 x 4 = 4 64 96 152 x 5 = 6 216 * Theo công thức (3.1.1) ta tính được 56- 48 48 f(–3) »-64 + 0, 5 + 0, 5(0, 5- 1) + 0, 5(0, 5 - 1)(0, 5 - 2) 1! 2! 3! 0 0 + 0, 5(0, 5 1)(0, 5 2)(0, 5 3) + =0, 5(0, 5 1)(0, 5 2)(0, 5 3)(0, 5 4) – 4! 5! 27 Program NEWTONTIEN; Var 45
A:Array[0 20] of Integer; D: Array[0 20,0 20] of real; X,Y,B: Array[0 20] of real; XNS,H,U,SUM: Real; I,J,N:Integer; BEGIN Write('So moc quan sat la N+1,N=');Readln(N); Writeln('Cac moc quan sat cach deu la'); For I:=0 to N do Begin Write('X[',I,']=');Readln(X[I]); end; Writeln('Cac gia tri quan sat tuong ung'); For I:=0 to N do Begin Write('Y[',I,']=');Readln(Y[I]); end; Write('Can noi suy gia tri ham tai XNS=');Readln(XNS); H:=X[1]-X[0];SUM:=0; For I:=0 to N do D[I,0]:=Y[I]; For J:=1 to N do For I:=0 to N-J do D[I,J]:=D[I+1,J-1]-D[I,J-1]; B[0]:=Y[0];U:=(XN–X[0])/H;A[1]:=1; For I:=2 to N do A[I]:=A[I–1]*I; For I:=1 to N do Begin B[I]:=1; For J:=1 to I do B[I]:=B[I]*(U–J+1);SUM:=SUM+B[I]*D[0,I]/A[I]; end; SUM:=SUM+B[0]; Writeln('Gia tri noi suy Y(‘,XNS.’) =',SUM); Readln; END. 3.2 Công thức nội suy Gregory-Newton lùi * Các mốc quan sát xk cách đều với h là khoảng cách giữa 2 mốc quan sát liên tiếp. * Các giá trị quan sát fk = f(xk). * Tính f(x) khi x ở khoảng cuối các mốc quan sát. 46
* Chọn x0 là một trong các mốc cuối quan sát và gần x nhất. * Tính u = x- x0 . h * Lập bảng số gia hữu hạn lùi 2 3 4 4 xk = x0 + k.h fk ∇fk ∇ fk ∇ fk ∇ fk ∇ fk . . . x - 4 f - 4 ∇f - 3 2 x - 3 f - 3 ∇ f - 2 3 ∇f - 2 ∇ f - 1 2 4 x - 2 f - 2 ∇ f - 1 ∇ f 0 3 5 ∇f - 1 ∇ f0 ∇ f 1 2 4 x - 1 f - 1 ∇ f0 ∇ f1 3 ∇f0 ∇ f1 2 x0 f 0 ∇ f1 ∇f1 x 1 f 1 . . . . . . . . . . . . trong đó Số gia hữu hạn lùi bậc một là ∇fk = fk – fk - 1 2 Số gia hữu hạn lùi bậc hai là ∇ fk = ∇fk – ∇fk - 1 3 2 2 Số gia hữu hạn lùi bậc ba là ∇ fk = ∇ fk – ∇ fk - 1 . . . . . . ÑÑff23 Ñ f ÑN f f() x»+ f00 u + uu( ++ 1) 0uu ( + 1)( u +++ 2) 0uu ( + 1)( u + 2) ( u +- N 1) (3.2.1) 0 1! 2! 3! N! Vế phải là đa thức bậc N. * Sai số là f()(N + 1) x N + 1 E =+ h u(u 1) . . . (u + N - 1) ; x < ξ < x N (N+ 1)! - N 0 (3.2.2) Ví dụ 3.2 Từ bảng số liệu sau đây hãy tính gần đúng giá trị f(2,2) xk – 2 – 1 0 1 2 3 fk 16 1 0 1 16 81 Giải * Ở đây ta thấy rằng xk là các mốc cách đều với h = 1 và x = 2,2 ở vào khoảng cuối bảng giá trị xk nên ta chọn x0 = 2. * u = x- x0 = 2,2- 2 = 0,2. h 1 47
* Bảng số gia hữu hạn lùi tương ứng 2 3 4 5 xk = 2 + k fk ∇fk ∇ fk ∇ fk ∇ fk ∇ fk x - 4 = – 2 16 – 15 x - 3 = – 1 1 14 – 1 – 12 x - 2 = 0 0 2 24 1 12 0 x - 1 = 1 1 14 24 15 36 x 0 = 2 16 50 65 x 1 = 3 81 * Theo công thức (3.2.1) ta tính được 15 14 12 24 f(2,2) »++++++++++»16 0, 2 0, 2(0, 2 1) 0, 2(0, 2 1)(0, 2 2) 0, 2(0, 2 1)(0, 2 2)(0, 2 3) 23,4256 1! 2! 3! 4! 3.3 Công thức nội suy Gauss * Các mốc quan sát xk cách đều với h là khoảng cách giữa 2 mốc quan sát liên tiếp. * Các giá trị quan sát fk = f(xk). * Tính f(x) khi x ở khoảng giữa các mốc quan sát. * Chọn x0 là một trong các mốc giữa quan sát và gần x nhất. * Tính u = x- x0 . h * Lập bảng số gia hữu hạn trung tâm 2 3 4 5 xk = x0 + k.h fk δfk + 0,5 δ fk δ fk + 0,5 δ fk δ fk + 0,5 . . . x - 2 f - 2 δf - 1,5 2 x - 1 f - 1 δ f - 1 3 δf - 0,5 δ f - 0,5 2 4 x 0 f 0 δ f 0 δ f 0 3 5 δf 0,5 δ f0,5 δ f 0,5 2 4 x 1 f 1 δ f 1 δ f1 3 δf 1,5 δ f1,5 2 x 2 f 2 δ f 2 δf 2,5 x 3 f 3 . . . . . . . . . . 48
trong đó Số gia hữu hạn trung tâm bậc một là δfk + 0,5 = fk + 1 – fk 2 Số gia hữu hạn trung tâm bậc hai là δ fk = δfk + 0,5 – δfk - 0,5 3 2 2 Số gia hữu hạn trung tâm bậc ba là δ fk + 0,5 = δ fk + 1 – δ fk 4 3 3 Số gia hữu hạn trung tâm bậc bốn là δ fk = δ fk + 0,5 – δ fk - 0,5 . . . . . . 2 3 ddff0,5 d f 0,5 f() x»+ f u +0 uu( -+ 1)uu ( - 1)( u + 1) 0 1! 2! 3! 4 5 d f d f0,5 N +-+-0 uu(1)(1)(2) u u ++ uu ( - 1)( u + 1)( u - 2)( u + 2) + +d 4! 5! (3.3.1) Vế phải là đa thức bậc N, trong đó ïì ï d 2M f ï 0 uu(-+ 1)( u 1) ( u -++ M 1)( u M 1)( u M ) , N = 2 M N ï (2M )! d = íï ï 21M + ï d f0,5 ï uu(-+ 1)( u 1) ( u - M )( u + M ) , N =+ 2 M 1 îï (2M + 1)! (3.3.2) * Sai số là ïì 21M + ï h (2M + 1) ï fuuuuMuMuMNM(x ) (-+ 1)( 1) ( -++ 1)( 1)( ) , = 2 ï (2M + 1)! E = íï x < ξ < x N ï 22M + - N N ï h (2M + 2) ï f(x ) uu (-+ 1)( u 1) ( uMuMuM - )( + )( 1) , N =+ 2 M 1 ïî (2M + 2)! (3.3.3) Ví dụ 3.3 Tính gần đúng giá trị f(2,6) từ bảng số liệu xk – 4 – 1 2 5 8 11 fk – 1 3 5 8 16 25 Giải * Ở đây ta thấy rằng xk là các mốc cách đều với h = 3 và x = 2,6 ở vào khoảng giữa các mốc quán sát nên ta chọn x0 = 2. * u = x- x0 = 2,6- 2 = 0,2. . h 3 * Bảng số gia hữu hạn trung tâm tương ứng 2 3 4 5 xk = 2 + 3k fk δfk + 0,5 δ fk δ fk + 0,5 δ fk δ fk + 0,5 x - 2 = – 4 – 1 4 x - 1 = – 1 3 – 2 2 3 x 0 = 2 5 1 1 49
3 4 – 9 x 1 = 5 8 5 – 8 8 – 4 x 2 = 8 16 1 9 x 3 = 11 25 * Theo công thức (3.3.1) ta tính được 31 4 f (2,6) »+5 0,2 + 0,2(0,2 - 1) + 0,2(0,2 - 1)(0,2 + 1) 1! 2! 3! 1 - 9 +-+-0, 2(0, 2 1)(0, 2 1)(0, 2 2). +-+-+»0, 2(0, 2 1)(0, 2 1)(0, 2 2)(0, 2 2) 5,349376 4! 5! 3.4 Công thức nội suy Lagrange * Các mốc quan sát bất kỳ xk , k = 0 , . . . , N * Các giá trị quan sát tương ứng fk= f(xk). * Tìm giá trị gần đúng của f(x) , x0 < x < xN . * Công thức nội suy Lagrange được cho như sau ()() ()xx 12 xx xx -NN ()() () xx -02 xx - xx - f ()xff»+01 (xxxx0102 )( ) ( xx 0 -N ) ( xxxx1 012 )( ) ( xx 1 -N ) (xx 013 )( xx )( xx ) ( xx -N ) ++f2 (xxxxxx202123 )( )( ) ( xx 2 -N ) (xx 02 ) ( xxNN- )( xx - ) (xx 01 )( xx ) ( xxN ) ++fNN- 1 f (xxxxxxNNNNNNNNN 10 ) ( 1 2)( 1 - ) ( xxxxxx 0 )( 1 ) (- 1 ) (3.4.1) Ta có thể rút ra cách đơn giản như sau để triển khai (3.4.1): Vế phải có N+1 số hạng Số hạng thứ k, k = 0, . . . , N có 2 thừa số: Thừa số thứ hai là fK còn thừa số thứ nhất là phân số mà tử số có N thừa số không có (x – xK), mẫu số cũng có N thừa số như tử số khi thay x bởi xK. * Sai số là f()(N + 1) x E = (x xxx )( ) ( xx ) ; x < ξ < x N (N+ 1)! 01 N 0 N (3.4.2) Ví dụ 3.4 Tìm giá trị gần đúng của f(2,6) từ bảng số liệu: xk – 4 –1 2 5 8 11 fk – 1 3 5 8 16 25 Giải Áp dụng công thức (3.4.1) với N = 5 , x0 = – 4 , . . . , x5 = 11 , x = 2,6 50
()()()()()2,6+- 1 2,6 2 2,6 - 5 2,6 - 8 2,6 - 11 ()()()()()2,6+- 4 2,6 2 2,6 - 5 2,6 - 8 2,6 - 11 f (2,6) »-()1 + .3 ()()()()()-+41 42 45 48 411 ()()()()()-+14 12 15 18 111 ()()()()()2,6++ 4 2,6 1 2,6 5 2,6 8 2,6 11 ()()()()()2,6+++ 4 2,6 1 2,6 2 2,6 8 2,6 11 + .5 + .8 ()()()()()24212528211++ ()()()()()54515258511++ ()()()()()2,6++ 4 2,6 1 2,6 2 2,6 5 2,6 11 ()()()()()2,6++ 4 2,6 1 2,6 2 2,6 5 2,6 8 + .16 + .25 » ()()()()()84818285811++ ()()()()()11++ 4 11 1 11 2 11 5 11 8 5,356123 Program LAGRANGE; Var X,F,Array[0 15] of real; J,K,N:Integer; XNS,SH,SUM:Real; BEGIN Write('So moc quan sat la N+1, N=');Readln(N); Writeln('Cac moc quan sat bat ky la'); For J:=0 to N do Begin Write('X[',J,']=');Readln(X[J]); end; Writeln('Cac gia tri quan sat tuong ung la'); For J:=0 to N do Begin Write('F[',J,']=');Readln(F[J]); end; Write('Can noi suy gia tri ham tai XNS=');Readln(XNS); SUM:=0; For K:=0 to N do Begin SH:=F[K]; For J:=0 to N do If J<>K then SH:=SH*(XNS–X[J])/(X[K] –X[J]);SUM:=SUM+SH; end; Writeln('Gia tri noi suy F(‘,XNS,’)=',SUM); Readln; END. 3.5 Công thức nội suy Newton * Các mốc quan sát bất kỳ xk , k = 0 , . . . , N. * Các giá trị quan sát tương ứng fk = f(xk). * Tìm giá trị gần đúng của f(x), x0 < x < xN. * Công thức nội suy Newton được cho như sau 51
f(x) ≈ D00 + D11(x – x0) + D22(x – x0)(x – x1) + . . . + DNN(x – x0)(x – x1) . . . (x – xN-1) (3.5.1) trong đó Dk 0 = fk , k = 0 , . . . , N (3.5.2) Dkj 111- D k j D = ; j =1 , . . . , N ; k = j , . . . , N kj xx- k kj- (3.5.3) * Sai số là f()(N + 1) x E = (x xxx )( ) ( xx ) ; x < ξ< x N (N+ 1)! 01 N 0 N (3.5.4) Ví dụ 3.5 Tính giá trị gần đúng của f(2,6) từ bảng số liệu xk – 4 –1 2 5 8 11 fk – 1 3 5 8 16 25 Giải Theo công thức (3.5.2) D00 = –1 , D10 = 3 , D20 = 5 , D30 = 8 , D40 = 16 , D50 = 25 Theo công thức (3.5.3) D10- D 00 31+ 4 j = 1 k = 1 D11 = = = k = 2 x10- x -+14 3 D20- D 10 53- 2 D21 = = = x21- x 21+ 3 D30- D 20 85- k = 3 D31 = = = 1 k = 4 x32- x 52- D40- D 30 16- 8 8 D41 = = = x43- x 85- 3 D50- D 40 25- 16 k = 5 D51 = = = 3 x54- x 11- 8 24 - D21- D 11 33 1 D31- D 21 j = 2 k = 2 D22 = = = - k = 3 D32 = x20- x 24+ 9 x31- x 2 1- 1 = 3 = 51+ 18 8 - 1 D41- D 31 3 5 D51- D 41 k = 4 D42 = = = k = 5 D52 = x42- x 82- 18 x53- x 8 3- 1 = 3 = 11- 5 18 52
11 51 + - D32- D 22 18 9 1 D42- D 32 18 18 j = 3 k = 3 D33 = = = k = 4 D43 = = x30- x 54+ 54 x41- x 81+ 2 = 81 15 - D52- D 42 18 18 2 k = 5 D53 = = = - x52- x 11- 2 81 21 22 - D43- D 33 81 54 1 D53- D 43 81 81 j = 4 k = 4 D44 = = = k = 5 D54 = = = x40- x 84+ 1944 x51- x 11+ 1 1 - 243 11 D54- D 44 243 1944 1 j = 5 k = 5 D55 = = = - x50- x 11+ 4 3240 Theo công thức (3.5.1) 41 1 fx()»-1 + (x + 4) - (x + 4)(x + 1) + (x + 4)(x + 1)(x - 2) 39 54 1 1 +++ (x 4)(x 1)(x 2)(x 5) -++ (x 4)(x 1)(x 2)(x 5)(x 8) ⇒ f(2,6) 1944 3240 » 5,349376 Program NEWTON; Var A,X,Y:Array[0 15] of real; I,J,N:Integer; XA,SUM:Real; BEGIN Write('So diem moc N+1,N=');Readln(N); Writeln('Cac moc quan sat bat ky va cac gia tri quan sat tuong ung'); For I:= 0 to N do Begin Write('X[',I,']=');Readln(X[I]); Write('Y[',I'']=');Readln(Y[I]); end; Write('XA=');Readln(XA); For I:= 0 to N do A[I]:=Y[I]; For J:= 1 to N do For I:=N downto J do A[I]:=(A[I]– A[I–1])/(X[I] –X[I–J]);SUM:=A[N]; For I:=N–1 downto 0 do SUM:=SUM*(XA–X[I])+A[I]; Writeln('Gia tri noi suy Y(‘,XNS,’)=',SUM); Readln; 53
END 3.6 Công thức bình phương nhỏ nhất * Các mốc quan sát bất kỳ xk , k = 0 , . . . , N. * Các giá trị quan sát tương ứng fk = f(xk). * Tìm giá trị gần đúng của f(x), x0 < x < xN * Ta tìm đa thức bậc M 2 PM(x) = a0 + a1x + a2x + . . . + aMxM (3.6.1) V v tnh f(x) » PM(x) , x0 < x < xN (3.6.2) Các hệ số của PM(x) được tìm từ hệ M + 1 phương trình , M + 1 ẩn số sau đây: ì ï ï aN(++ 1) aåå x + a x2 + + a åå xM = f ï 012kkMkk ï kk kk ï ï 23M + 1 ï axaxax012åååkkk++ ++ a Mkkk å x = å xf í kkk k k ï ï ï ï ï MM++12 M 2 MM ï ax01ååkk++ ax ax 2 å k ++= a Mkkk åå x xf îï kk k kk (3.6.3) Khi M = 1 ta có ì ï ï aN(1)++ aåå x = f ï 01kk íï kk ï ï axax+=2 xf ï 01åååkkkk îï kkk (3.6.4) Khi M = 2 ta có hệ : ïì ï 2 ï aN012(1)++ aååå xkkk + a x = f ï kkk ï ï íï axaxax++23 = xf ï 012ååååkkkkk ï kkkk ï ï 23 42 ï axaxax012ååååkk++ kkk = xf îï kk kk (3.6.5) Ví dụ 3.6 Tính giá trị gần đúng của f(2,6) với M =1 và M = 2 từ bảng số liệu : xk – 4 – 1 2 5 8 11 fk – 1 3 5 8 16 25 Giải * Với M = 1 theo công thức (3.6.4) ta cần giải hệ 54
ïì aa(51)++ (4125811) ++++ =-++++ 13 58 16 + 25 ï 01 í ï éù222222 2 ï aa01(- 4 - 1 + 2 + 5 + 8 + 11) +êú()() - 4 +- 1 + 1 + 2 + 5 + 8 + 11 =()()() - 4 - 1 +- 1 .3 + 2.5 + 5.8 + 8.16 + 11.25 îï ëû ïì 6a += 21a 56 621 ÛÞíï 01 D==6.231- 21.21= 945 ï îï 21a01 += 231a 454 21 231 56 21 6 56 D= =56.231 - 454.21 = 3402 ; D= =6.454 - 21.56 = 1548 aa01454 231 21 454 Δa0 Δa1 3402 1548 Þ P1(x) = +=+x x Þ P1(x) = 3,6 +1,638095x ΔΔ 945 945 Vậy f(2,6) » P1(2,6) » 7,859048 * Với M = 2 theo công thức (3.6.5) giải hệ ì ï 6a01 ++ 21a 231a 2 = 56 ï í 21a ++ 231a 1911a =Þ 454 ï 01 2 ï îï 231a01 ++ 1911a 19635a 2 = 4256 6 21 231 6.231.19635++ 21.1911.231 21.1911.231 D= 21 231 1911 = 231.231.231 21.21.19635 6.1911.1911 231 1911 19635 =-27214110+9270261+9270261 12326391- 8659035- 21911526 = 2857680 56 21 231 56.231.19635++ 454.1911.231 21.1911.4256 D=454 231 1911 = a0 4256.231.231 21.454.19635 56.1911.1911 4256 1911 19635 =-253998360+200414214+170797536 227104416- 187200090- 204507576 = 6398028 6 56 231 6.454.19635++ 21.4256.231 56.1911.231 D=21 454 1911 = a1 231.454.231 21.56.19635 6.4256.1911 231 4256 19635 =-53485740+20645856+24720696 24225894- 23090760- 48799296 = 2736342 62156 6.231.4256++ 21.1911.56 21.454.231 D=21 231 454 = a2 231.231.56 21.21.4256 6.1911.454 231 1911 4256 =-5898816+2247336+2202354 2988216- 1876896- 5205564= 277830 1122 P2(x) = ()Δ+Δ+Δaaxax01.2. = ()6398028+ 2736342x + 277830 x Δ 2857680 » 2,238889 + 0,957540x + 0,097222x2. Vậy f(2,6) » P2(2,6) » 5,385714 Program BPNN; Var A:Array[0 15,0 15] of real; X,Y,Z,COEF:Array[0 15] of real; I,J,K,K1,M,M1,M2,N,DK:Integer; X0,Y0,P0,AB,ABMAX,SUM:Real; BEGIN 55
Write('So diem quan sat la N+1,N=');Readln(N); Write('Bac cua da thuc noi suy la M=');Readln(M); M1:=2*M;M2:=M+1; Writeln('Cac moc quan sat bat ky va cac gia tri quan sat tuong ung'); For K:= 0 to N do Begin Write('X[',K,']=');Readln(X[K]); Write('Y[',K,']=');Readln(Y[K]); end; For I:= 0 to M do A[I,M2]:=0; For K:= 0 to N do Begin X0:=X[K];Y0:=Y[K];P0:=1; For I:= 0 to M do Begin A[I,M2]:=A[I,M2]+Y0*P0; P0:=P0*X0; end; end; COEF[0]:=N+1; For J:=1 to M1 do COEF[J]:=0; For K:= 0 to N do Begin X0:=X[K];P0:=X[K]; For J:=1 to M1 do Begin COEF[J]:=COEF[J]+P0;P0:=P0*X0; end; end; For I:= 0 to M do For J:= 0 to M do A[I,J]:=COEF[I+J]; DK:=1; For I:= 0 to M-1 do If DK=1 then Begin AB:=0; For K:=1 to M do Begin ABMAX:=ABS(A[K,I]); If AB<ABMAX then Begin 56
AB:=ABMAX;K1:=K; end; end; If AB = 0 then DK:=2 Else Begin If K1<>I then Begin For J:=I to M2 do Begin P0:=A[I,J];A[I,J]:=A[K1,J];A[K1,J]:=P0; end; end; P0:=A[I,I]; For K:=I+1 to M do For J:=I+1 to M2 do A[K,J]:=A[K,J]-A[K,I]*A[I,J]/P0; end; end; If A[M,M]=0 then DK:=2 Else Z[M]:=A[M,M2]/A[M,M]; For I:=M–1 downto 0 do Begin SUM:=0;P0:=A[I,I]; For J:=I+1 to M do SUM:=SUM-A[I,J]*Z[J];Z[I]:=(SUM+A[I,M2])/P0; end; end; Writeln('Cac he so theo luy thua tang dan cua da thuc noi suy la') For I:= 0 to M do Write('Z[',I,']=',Z[I]); Readln; END. Bài tập 1- Cho bảng số liệu xk – 2 0 2 4 6 8 10 fk 1 5 8 12 18 20 25 Tính gần đúng f(– 1,5), f(– 0,4), f(3,5) và f(8,2). −−−1,5 ( 2) Đs: x = – 2 , h = 2 , u ==0, 25 ⇒ f(– 1,5) ≈ 1,49112 . 0 2 −−0,4 0 x = 0 , h = 2 , u ==−0, 2 ⇒ f(– 0,4) ≈ 2,75648 . 0 2 3,5− 4 . x = 4 , h = 2 , u ==−0, 25 ⇒ f(3,5) ≈ 10,686432 . 0 2 57
8,2− 8 . x = 8 , h = 2 , u ==0,1⇒ f(8,2) ≈ 19,393337 . 0 2 2- Tìm xấp xỉ của f(0,5) , f(1,5) , f(– 0,5) từ bảng số liệu xk – 4 – 2 – 1 0 1 3 fk – 1 0 2 5 8 10 a- Bằng phương pháp nội suy Newton . (f(0,5) ≈ 6,558594) b- Bằng phương pháp bình phương nhỏ nhất với M = 1 và M = 2. * f(x) ≈ 4,881356 + 1,762712x. f(0,5) ≈ 5,762712 , f(1,5) ≈ 7,525424 , f(– 0,5) ≈ 4. * f(x) ≈ 4,542646 + 1,835293x + 0,072581x2 , f(– 0,5) ≈ 3,643144. c- Bằng phương pháp Lagrange. f(0,5) ≈ 6,75366 3- Cho bảng số liệu xk(độ cao : m) 0 400 800 1200 1600 2000 fk(áp suất: mm thuỷ ngân) 750 675 610 580 540 500 Tính áp suất ở độ cao 125m, 900m, 1950m . 4- Cho bảng số liệu xk(khoảng cách đến nhà máy:m) 50 100 200 250 350 500 3 fk(gam bụi than/10m không khí) 150 120 70 50 40 30 Tính lượng bụi than trong 10m3 không khí ở vị trí cách nhà máy 300m ,150m , 400m . a- Bằng phương pháp nội suy Newton. b- Bằng phương pháp bình phương nhỏ nhất với M = 1 và M = 2. c- Bằng phương pháp Lagrange. 5- Cho bảng số liệu xk(độ sâu: m) 10 20 30 40 50 60 3 fk(gam phù sa/1m nước) 1 3 7 10 12 14 Tính lượng phù sa trong 1m3 nước ở độ sâu 32m . (7,697856 g) 6- Cho bảng số liệu về độ sâu của hồ (đơn vị là m) xk (khoảng cách từ bờ) 5 15 25 35 45 55 fk (độ sâu) 2,1 9,3 11,2 15,4 18,6 22,5 Tìm độ sâu của hồ tại điểm cách bờ 50m . (19,726171 m) Kiểm tra nhận thức 1* Nêu những điểm giống nhau, khác nhau của ba phương pháp nội suy Gregory- Newton tiến, lùi và Gauss. 2* Viết chương trình PASCAL cho các phương pháp Gregory-Newton lùi và Gauss. 3* Nêu ưu điểm, nhược điểm của từng phương pháp Lagrange, Newton và bình phương nhỏ nhất. 4* Hãy chạy chương trình PASCAL đã có để giải các ví dụ. Nếu có điều gì chưa thoả mãn thì chủ động viết lại chương trình (có thể bằng ngôn ngữ lập trình bậc cao khác) và giải lại các ví dụ, các bài tập. 58
Chương 4 : XẤP XỈ ĐẠO HÀM VÀ TÍCH PHÂN XÁC ĐỊNH Mặc dù ta có thể tính đạo hàm một cách dể dàng, nhưng khi tìm giá trị của nó tại một điểm cụ thể lại bắt buộc ta phải tính gần đúng, còn việc tính tích phân xác định không phải lúc nào cũng đơn giản vì không phải lúc nào cũng tìm được nguyên hàm dưới dạng một biểu thức giải tích. Vì vậy ta sẽ đề cập trong chương này một số phương pháp tính gần đúng đạo hàm tại một điểm và tích phân xác định. 4.1 Xấp xỉ giá trị đạo hàm theo tỷ sai phân Df Dựa vào định nghĩa fx¢()= lim D®x 0 Dx ta có giá trị gần đúng khi Δf = f(x + h) – f(x – h) , Δx = 2h f ()()xh+- fxh - fx¢()» 2h (4.1.1) với sai số là h2 Ef= (3)()x , x – h < ξ < x + h 1 6 (4.1.2) Suy ra từ (4.1.1) -++ f(x 2h) 8f(x + +- h) 8f(x h) f(x 2h) f(x)¢ » 12h (4.1.3) với sai số là h4 Ef=-6(5) (x ) 648 f (5) (x ) , x – 2h < ξ , ξ < x + 2h 2 1440 12 1 2 (4.1.4) Ví dụ 4.1 Cho bảng số liệu xk fk 1,01 0,004321 1,02 0,008600 1,03 0,012837 1,04 0,017033 1,05 0,021189 1,06 0,025306 1,07 0,029384 1,08 0,033424 1,09 0,037426 1,10 0,041393 Tính gần đúng f '(1,06) Giải * Theo công thức (4.1.1) Với h = 0,01 59
ff(1,06+- 0,01) (1,06 - 0,01) 0,029384- 0,021189 f ¢(1, 06) » = » 0,40495 2.0,01 0,02 Với h = 0,02 ff(1,06+− 0,02) (1,06 − 0,02) 0,033424− 0,017033 f ′(1, 06) ≈ = = 0,409775 2.0,02 0,04 Với h = 0,03 ff(1,06+− 0,03) (1,06 − 0,03) 0,037426− 0,012837 f ′(1, 06) ≈ = = 0,409817 2.0,03 0,06 Với h = 0,04 ff(1,06+− 0,04) (1,06 − 0,04) 0,041393− 0,008600 f ′(1, 06) ≈ = = 0,409913 2.0,04 0,08 * Theo công thức (4.1.3) Với h = 0,01 -+ f(1,08) 8f(1,07) - 8f(1,05) + f(1,04) f '(1,06) » 12.0,01 = -+ 0,033424 8.0,029384 - 8.0,021189 + 0,017033 » 0,409741 0,12 Với h = 0,02 -+ f(1,10) 8f(1,08) - 8f(1,04) + f(1,02) f '(1,06) » 12.0,02 = -+ 0,041393 8.0,033434 - 8.0,017033 + 0,008600 = 0,409729 0,24 4.2 Xấp xỉ giá trị đạo hàm theo công thức Richardson * Cho các mốc quan sát cách đều xk = x0 + kh , k = – N , . . . , N , trong đó h là khoảng cách giữa 2 mốc quan sát liên tiếp. * Các giá trị quan sát tương ứng fk = f(xk). Tính f '(x) với x là một trong các mốc quan sát. m * h0 = h.2 trong đó m ≥ 1 , nguyên và chọn m lớn nhất sao cho x + h0 và x – h0 trùng với các mốc quan sát và không vượt ra ngoài các mốc quan sát. Lúc này đảm bảo sai số nhỏ nhất. * Tính h0 Hk = , k = 1,2, . . . , m + 1 2k- 1 (4.2.1) f(x+ Hk )fx ( Hk ) Dk1 = 2Hk (4.2.2) j- 1 4 DDkj 111- k j D = , j = 2 , . . . , k k j 41j- 1- (4.2.3) f ’(x) » Dkk (4.2.4) 60
với sai số 2k 2k Ek = Hằng số.Hk º OH()k (4.2.5) Ví dụ 4.2 Cho bảng số liệu xk fk 1,01 0,004321 1,02 0,008600 1,03 0,012837 1,04 0,017033 1,05 0,021189 1,06 0,025306 1,07 0,029384 1,08 0,033424 1,09 0,037426 1,10 0,041393 Tính gần đúng giá trị f'(1,06) . Giải 2 Ở đây x = 1,06, h = 0,01 nên ta chọn h 0 = 0,01.2 = 0,04 (m = 2) là tốt nhất h0 k = 1 H 1 = = 0,04 211- f(1,06+- 0,04)f (1,06 - 0,04) 0,041393- 0,008600 D = = » 0,409913 11 2.0,04 0,08 h0 k = 2 H 2 = = 0,02 221- f(1,06+- 0,02)f (1,06 - 0,02) 0,033424- 0,017033 D = = » 0,409775 21 2.0,02 0,04 421- DD- 4.0,409775- 0,409913 j = 2 D = 21 11 = » 0,409729 22 4121- - 3 h0 k = 3 H 3 = = 0,01 231- f(1,06+- 0,01)f (1,06 - 0,01) 0,029384- 0,021189 D = = » 0,409750 31 2.0,01 0,02 421- DD- 4.0,409750- 0,409775 j = 2 D = 31 21 = » 0,409742 32 4121- - 3 431- DD- 16.0,409742- 0,409729 j = 3 D = 32 22 = » 0,409743 33 4131- - 15 f ’(1,06) » D33 » 0,409743 Chú ý * Nếu chọn h 0 =0,01.2 (m = 1 thì chỉ tính được D22 nên kém chính xác hơn. 3 * Nếu chọn h 0 = 0,01.2 (m = 3) thì x + H1 = 1,06 + 0,08 = 1,12 vượt ra ngoài mốc quan sát (và cũng như vậy với m ≥ 4 ) nên không tính được f’(x) . 61
Program RICHARDSON; Var D:Array[1 20,1 20] of real; J,K,N:Integer; DK,H,X0,EPS:Real; Function F(X:Real):Real; Begin F:=ln(X)/ln(10); end; {Tìm đạo hàm của F(X)= log10X = ln(X)/ln(10) tại X0 tương ứng} BEGIN Write(' Tinh dao ham cua F(X) tai X0 =');Readln(X0); Write(' Sai so EPS =');Readln(EPS); H:=1;DK:=1; D[1,1]:=0.5*(F(X0+H)-F(X0-H))/H;J:=2; While ((J EPS)) do Begin H:=H/2;D{J,1]:=0.5*(F(X0+H)-F(X0-H))/H; For K:=2 to J do D[J,K]:=D[J,K–1]+(D[J,K–1] –D[J–1,K–1])/(EXP(K*LN(4)) –1); DK:=ABS(D[J,J] –D[J–1,J–1]);N:=J;J:=J+1; end; If DK<EPS then Writeln('Dao ham cua F(',X0,')=',D[N,N]) ELSE Writeln('Chua hoi tu'); Readln; END. 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 a- Công thức nội suy Gregory-Newton tiến Khi x ở khoảng đầu các mốc quan sát ' é DDff23 D f D N f ù 1 fx'( ) »+ê f00 u + uu() -+1 0 uu()() 1 u 2 ++ 0 uu()()( -1 u - 2 u - N + 1 )ú . (4.3.1) ê 0 ú h ë 1! 2! 3! N! ûu với sai số là Oh( N ). b- Công thức nội suy Gregory-Newton lùi Khi x ở khoảng cuói các mốc quan sát ' é ÑÑff23 Ñ f Ñ N f ù 1 f'( x ) »+ê f00 u + uu() ++1 0uu()() ++1 u 2 ++ 0 uu()()( +1 u + 2 u +- N 1 )ú . (4.3.2) ê 0 ú h ë 1! 2! 3! N! ûu với sai số là Oh( N ). c- Công thức nội suy Gauss Khi x ở khoảng giữa các mốc quan sát ' é 2 3 N ù ê ddff0,5 d f 0 0,5 d f 0(0,5) ú 1 fx'( ) »+ f0 u + uu() -+1 uu()() -+1 u 1 ++ uu()() -1 u + 1 ú . (4.3.3) ê 1! 2! 3! N! h ëê ûúu với sai số là Oh( N ). 62
Ví dụ 4.3 Cho bảng số liệu : xk fk 1,06 0,025306 1,07 0,029384 1,08 0,033424 1,09 0,037426 1,10 0,041393 1* Tính gần đúng giá trị của f '(1,065). Giải * h = 0,01, x = 1,065 ở vào khoảng đầu các mốc quan sát nên ta chọn x0 = 1,06. * Bảng số gia hữu hạn tiến 2 3 xk = 1,06 + 0,01k fk Δfk Δ fk Δ fk 4 Δ fk x0 = 1,06 0,025306 0,004078 x1 = 1,07 0,029384 – 0,000038 0,004040 0 x2 = 1,08 0,033424 – 0,000038 0,000003 0,004002 0,000003 x3 = 1,09 0,037426 – 0,000035 0,003967 x4 = 1,10 0,041393 ' é DDff23 D f D 4 f ù 1 * f ’(x) » êf++00 u uu(1) -+ 0 uu (1)(2) u + 0 uu(1)(2)(3). u u ú ê 0 ú h ë 1! 2! 3! 4! ûu ' é DDff23 D f D 4 f ù 1 = êf ++00uuuuuu()23 -+ 0 (32) -+2+-+-0 (6116).uu43 uu 2 ú ê 0 ú h ë 1! 2! 3! 4! ûu é 23 4 ù êDDff00 D f 02 D f0 32 ú 1 = ê +-+-+(2uuu 1) (3 6 2) +-+-(4uu 18 22 u 6)ú . ë 1! 2! 3! 4! û h xx- 1,065- 1,06 Thay h = 0,01 và u = 0 = = 0,5 ta tính được h 0,01 é 0,000038 0 0,000003 ù 1 f '(1,065) »-ê0,004078() 2.0,5 -+-+ 1() 3.0,52 6.0,5 2 + ()4.0,532 18.0,5 +22.0,5 6 ú ëê 2! 3! 4! ûú0,01 » 0,4078 2* Tính gần đúng giá trị của f '(1,092). Giải * h = 0,01, x = 1,092 ở vào khoảng cuối các mốc quan sát nên ta chọn x0 = 1,09. * Bảng số gia hữu hạn lùi 63
2 3 4 xk = 1,06 + 0,01k fk ∇fk ∇ fk ∇ fk ∇ fk x- 3 = 1,06 0,025306 0,004078 x- 2 = 1,07 0,029384 – 0,000038 0,004040 0 x- 1 = 1,08 0,033424 – 0,000038 0,000003 0,004002 0,000003 x0 = 1,09 0,037426 – 0,000035 0,003967 x1 =1,10 0,041393 ' éùÑÑff23 Ñ f 1 * f ’(x) » êúf++00 u uu(1) ++++ 0 uu (1)(2). u êú0 h ëû1! 2! 3! u ' éùÑÑff23 Ñ f 1 = êúfu++00() uuuuu23 ++++ 0 (32).2 êú0 h ëû1! 2! 3! u éù23 1 êúÑÑff00 Ñ f 02 = êú+++++(2uuu 1) (3 6 2 . ëû1! 2! 3! h x x 1,092 1,09 Thay h = 0,01 và u == 0 = 0,2 ta tính được h0,01 é ê 0,000038 0 2 ù 1 f '(1,092) »-ê0,004002() 2.0,2 ++++ 1() 3.0,2 6.0,2 2 ú » 0,39754 ë 2! 3! û0,01 3* Tính gần đúng giá trị của f '(1,077). Giải * h = 0,01 , x = 1,077 ở vào khoảng giữa các mốc quan sát nên ta chọn x0 = 1,08. * Bảng số gia hữu hạn trung tâm 2 3 xk = 1,08 + 0,01k fk δfk + 0,5 δ fk δ fk + 0,5 4 δ fk x- 2 =1,06 0,025306 0,004078 x- 1 =1,07 0,029384 – 0,000038 0,004040 0 x 0 =1,08 0,033424 – 0,000038 0,000003 0,004002 0,000003 x 1 =1,09 0,037426 – 0,000035 0,003967 x 2 =1,10 0,041393 64
3 ' é ddffd2 f d4 f ù 1 * f ’(x) » êf++0,5 u0 uu(1) -+-+0,5 uu (1)(1) u +-+-0 uu(1)(1)(2). u u ú ê 0 ú h ëê 1! 2! 3! 4! ûu 3 ' éùddffdd24ff1 = êúfuuuuuuuuu+0,5 +00()234 -+0,5 () -+ (2 -3 -+2 2). êú0 1! 2! 3! 4! h ëûêúu 3 ù éd f d 2 f d f d f 1 = ê 0,5 +-+-0 ()21uu0,5 32 1+ +0 46uu3 2 22 uú ê () ()úh ëê 1! 2! 3! 4! ûú x x 1,077 1,08 Thay h = 0,01 và u == 0 = – 0,3 ta tính được h0,01 é 0,000038 0,000003 f '(1,077)» ê0,004002 + ()2.() 0,3 1 3.() 0,32 1 ëê 2! 3! () 0,000003 ù 1 + +()4.( 0,3)32 6.( 0,3) 2.( 0,3) 2 ú » 0,403226 4! ûú0,01 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ỳ a- Công thức nội suy Lagrange (xx 12 )( xx20 ) ( xx -NN ) ( xx - )( xx - ) ( xx - ) fx( )»++f01f (xxxx0102 )( ) ( xx 0 -NN ) ( xxxx1012 )( ) ( xx1 - ) (xx 02 ) ( xxNN- )( xx - ) (xx 01 )( xx ) ( xxN ) ++fNN- 1 f (xNNNN 10121 x ) ( x x)( x - xN ) ( xxxx NNNN 0 )( 1 ) ( xx - 1 ) (4.4.1) Sau đó ta tính f ’(x). b- Công thức nội suy Newton f(x) ≈ D00 + D11(x – x0) + D22(x – x0)(x – x1) + . . . + DNN(x – x0)(x – x1) . . . (x – xN–1) (4.4.2) trong đó DK0 = fK , k = 0 , . . . , N DDkj 111- k j Dkj = ; j =1 , . . . , N ; k = j , . . . , N xxk - kj- Sau đó ta tính f ’(x). c- Công thức bình phương nhỏ nhất Tính f(x) theo công thức bình phương nhỏ nhất. Sau đó tính f ’(x). Ví dụ 4.4 Cho bảng số liệu : xk – 2,1 – 1,5 – 0,5 0,4 1 fk 12 9,8 6,2 5 4,3 a- Bằng công thức Lagrange tính f ’(0,3) b- Bằng công thức Newton.tính f ’(0) c- Bằng công thức bình phương nhỏ nhất với M = 1 tính f ’(x) Giải (xx++ 1,5)( 0,5)( x - 0,4)( x - 1) a- f(x) » .12 (-+ 2,1 1,5)( -+ 2,1 0,5)( 2,1 0,4)( 2,1 1) (xxxx++- 2,1)( 0,5)( 0,4)( - 1) + .9, 8 (-+ 1,5 2,1)( -+ 1,5 0,5)( 1,5 0,4)( 1,5 1) (xxxx++- 2,1)( 1,5)( 0,4)( - 1) + .6, 2 (-+ 0,5 2,1)( -+ 0,5 1,5)( 0,5 0,4)( 0,5 1) 65
(xxxx+++- 2,1)( 1,5)( 0,5)( 1) + .5 (0,4+++- 2,1)(0,4 1,5)(0,4 0,5)(0,4 1) (xxxx+++- 2,1)( 1,5)( 0,5)( 0,4) + .4, 3 (1+++- 2,1)(1 1,5)(1 0,5)(1 0,4) 12 =+ +()x43 0,6x 1,65x 2 0,25x 0,3 7,44 9,8 -+ +()x43 1,2x 2,19x 2 0,43x 0,42 2,85 6,2 ++ +()x43 2,2x 1,49x 2 2,97x 1,26 2,16 5 -++-+()x43 3,1x 0,85x 2 3,375x 1,1575 2,565 4,3 +++ ()x43 3,7x 3,31x 2 0,405x 0,63 6,975 12 ⇒ fx¢( )»+ () 4x32 1,8x 3,3x 0,25 7,44 9,8 -+ ()4x32 3,6x 4,38x 0,43 2,85 6,2 ++ ()4x3 6,6x2 2,98x 2,97 2,16 5 -++-()4x32 9,3x 1,7x 3,375 2,565 4,3 +++-()4x32 11,1x 6,62x 0,405 6,975 12 ⇒ f ’(0,3) »+ ()4.0,332 1,8.0,3 3,3.0,3 0,25 7,44 9,8 -+ ()4.0,332 3,6.0,3 4,38.0,3 0,43 2,85 6,2 ++ ()4.0,332 6,6.0,3 2,98.0,3 2,97 2,16 5 -++-()4.0,332 9,3.0,3 1,7.0,3 3,375 2,565 4,3 +++-()4.0,332 11,1.0,3 6,62.0,3 0,405 » – 0.72938 6,975 b- D00 = 1,2 ; D10 = 9,8 ; D20 = 6,2 ; D30 = 5 ; D40 = 4,3 ; D10- D 00 9,8 - 12 j = 1 k = 1 D11 = = » – 3,666667 x10- x -+1,5 2,1 D20- D 10 6,2- 9,8 k = 2 D21 = = = – 3,6 x21- x -+0,5 1,5 D30- D 20 56,2- k = 3 D31 = = » – 1,333333 x32- x 0,4+ 0,5 D40- D 30 4,3- 5 k = 4 D41 = = » – 1,166667 x43- x 10,4- D21- D 11 -+3,6 3,666667 j = 2 k = 2 D22 = = » 0,041667 x20- x -+0,5 2,1 D31- D 21 -+1,333333 3,6 k = 3 D32 = = » 1,192983 x31- x 0,4+ 1,5 D41- D 31 -+1,166667 1,333333 k = 4 D42 = = = 0,083333 x42- x 10,5+ 66
D32- D 22 1,192983- 0,041667 j = 3 k = 3 D33 = = » 0,460526 x30- x 0,4+ 2,1 D42- D 32 0,083333- 1,192983 k = 4 D43 = = = – 0,44386 x41- x 11,5+ D43- D 33 0,44386 0,460526 j = 4 k = 4 D44 = = = – 0,291737 x40- x 12,1+ f(x)» D00+ D11(x–x0) + D22(x–x0)(x–x1) + D33(x–x0)(x–x1)(x–x2) + D44(x– x0)(x–x1)(x–x2)(x–x3) = 1,2 – 3,666667(x + 2,1) + 0,041667(x + 2,1)(x + 1,5) + 0,460526(x + 2,1)(x + 1,5)(x + 0,5) – 0,291737(x + 2,1)(x + 1,5)(x + 0,5)(x – 0,4) = 1,2 – 3,666667(x + 2,1) + 0,041667(x2 + 3,6x + 3,15) + 0,460526(x3 + 4,1x2 + 4,95x + 1,575) – 0,291737(x4 + 3,7x3 + 6,361x2 – 0,405x – 0,63) f '(x) » – 3,666667 + 0,041667(2x + 3,6) + 0,460526(3x2 + 8,2x + 4,95) – 0,291737(4x3 + 11,1x2 + 12,722x –0,405) f '(0) » – 3,666667 + 0,041667(2.0 + 3,6) + 0,460526(3.02 + 8,2.0 + 4,95) – 0,291737(4.03 + 11,1.02 + 12,722.0 –0,405) » – 1,118909 c- Với M = 1 theo công thức (3.6.4) ta cần giải hệ ïì ï aa01(4++ 1) ( - 2,1 - 1,5 - 0,5 + 0,4 + 1) = 12 + 9,8 + 6,2 ++ 5 4,3 ï ï éù22 222 í aa01( 2,1 1,5 0,5 + 0,4 ++- 1)êú()()() 2,1 +-+-+ 1,5 0,5 0,4 + 1 ï ëû ï îï =-()2,1 .12 +- () 1,5 .9,8 +- () 0,5 .6,2 + 0,4.5 + 1.4,3 ïì 5a -= 2,7a 37,3 Û íï 01 ï îï -+= 2,7a01 8,07a 36,7 2,7 PT12+® PT PT 2 5 ì ï 5a01 -= 2,7a 37,3 Û í ï îï 6,612a1 = 56,842 ïì a» 12,102269 Þ íï 0 ï îï a1 » 8,596794 Þ f(x) » 12,102269 + 8,596794x Þ f ’(x) » 8,596794 Program DHNEWTON Var A,X,F:Array[0 15] of real; J,K,N:Integer; P,X0,DF:Real; BEGIN Write('So diem moc quan sat la N+1, N=');Readln(N); Writeln('Cac moc quan sat bat ky la'); For J:=0 to N do Begin 67
Write('X[',J,']=');Readln(X[J]); end; Writeln('Cac gia tri quan sat tuong ung’ ); For J:=0 to N do Begin Write('F[',J,']=');Readln(F[J]); end; For K:=0 to N do A[K]:=F[K]; For J:=1 to N do For K:=N downto J do A[K]:=(A[K] –A[K–1])/(X[K] –X[K-J]); X0:=X[0];DF:=A[1];P:=1; For K:=2 to N do Begin P:=P*(X0–X[K–1]);DF:=DF+P*A[K]; end; Writeln('Gia tri gan dung cua dao ham tai X[0] la',DF); Readln; END. y B 4.5 Xấp xỉ giá trị tích phân xác định A Xuất phát từ ý nghĩa hình học của tích phân xác định : y = f(x) b ∫ f ()xdx= Diện tích hình thang cong aABb 0 a a b x a- Công thức hình thang ba- * Chia [a ; b] thành N đoạn nhỏ bằng nhau bởi các điểm chia x = a + kh : h = k N , k = 0 , . . . , N . b h * ò f ()xdx» [f(x0) + f(xN) + 2(f(x1) + . . . + f(xN - 1))] a 2 (4.5.1) * Sai số là hf2 ¢¢()x E = ()ba- , a < ξ < b N 12 (4.5.2) b- Công thức Simpson ba- * Chia [a ; b] thành 2N đoạn nhỏ bằng nhau bởi các điểm chia x = a + kh : h = , k 2N k = 0 , . . . , 2N . 68
b h * ò f ()xdx » [f(x0) + f(x2N) + 4(f(x1) + . . . + f(x2N - 1)) + 2(f(x2) + . . . + f(x2N - 2))] a 3 (4.5.3) * Sai số là hf4(4)()x E = ()ba- , a < ξ < b 2N 180 (4.5.4) Ví dụ 4.5 2 sin x 1* Tính gần đúng ò dx theo công thức hình thang 1 x với các điểm chia xk =1 + 0,1k ; k = 0 , . . . , 10 ( N = 10 Û h = 0,1 ) và tính sai số. Giải 2 sin x 0,1 ésin1 sin2æö sin1,1 sin1,2 sin1,9 ù dx ê ++2ç + ++ ÷ú * ò » ê ç ÷ú 1 x 2 êë 1 2èø 1,1 1,2 1,9 ûú » 0,05[0,841471 + 0,454649 + 2(0,810189 + 0,776699 + 0,741199 + 0,703893 + 0,664997 + 0,624734 + 0,583332 + 0,541026 + 0,498053)] » 0,659218 22 0,1 (-+xxxx 2)sin - 2 cos 1 * Sai số là E = (2-£ 1) ; 1 < ξ < 2 . 10 12 1200 3,8 2 2* Tính gần đúng ò edxx theo công thức Simpson 2 với các điểm chia xk = 2 + 0,15k , k = 0, . . . ,12 (2N = 12 Û h = 0,15) . Giải 3,8 2 0,15 23,822 22222 2 ò edxx » ee++4() e2,15 + e 2,45 + e 2,75 + e 3,05 + e 3,35 + e 3,65 2 3 22222 +++++2()eeeee2,32,62,93,23,5 » 0,05 [54,598150 + 1867292,353035 + 4(101,748085 + 404,438627 + 1924,651132 +10965,398466 + 74794,527535 + 610784,820388) + 2(198,343425 + 862,642196 + 4491,760512 + 28001,125926 + 208981,288870)] » 257414,180497 3* Tính gần đúng diện tích hình thang cong aABb sau đây : k 0 1 2 3 4 5 6 y xk 1 1,2 1,4 1,6 1,8 2 2,2 A fk 2 1,2 1,6 1 1,1 1 1,6 2 1,6 B 1,2 1 a b 69
0 1 1,4 1,8 2,2 x Giải Theo công thức hình thang S » 0,2 [2 + 1,6 + 2(1,2 + 1,6 + 1 + 1,1 + 1)] = 1,54 (đvdt) aABb 2 Theo công thức Simpson S » 0,2 [2 + 1,6 + 2(1,6 + 1,1) + 4(1,2 + 1 + 1)] » 1,4533 (đvdt) aABb 3 Program HINHTHANG; Var N,K: Integer; A,B,H,X,SUM: Real; Function F(X:Real):Real; Begin F:= Sin(X)/X; End; { Tính tích phân xác định của hàm f(x) = sin x trên [a;b]. } x BEGIN Write('Can duoi A=');Readln(A); Write('Can tren B=');Readln(B); Write('So diem chia N+1,N=');Readln(N); H:=(B– A)/N;SUM:=0; For K:=1 to N–1 do Begin X:=A+K*H;SUM:=SUM+F(X); End; SUM:=(2*SUM+F(A)+F(B))*H/2; Write('Gia tri xap xi cua tich phan la',SUM); Readln; END. Program SIMPSON; Var N,K:Integer; A,B,H,X,TC,TL,SUM:Real; Function F(X:Real):Real; Begin F:=EXP(X*X); End; b x2 { Tính tích phân xác định ∫ edx. } a BEGIN Write('Can duoi A=');Readln(A); Write('Can tren B=');Readln(B); 70
Write('So diem chia 2N+1,N=);Readln(N); H:=(B– A)/(2*N);SUM:=0; TC:=0; For K:=1 to N–1 do Begin X:=A+2*K*H;TC:=TC+F(X); End; TL:=0; For K:=1 to N do Begin X:=A+(2*K–1)*H;TL:=TL+F(X); End; SUM:=(F(A)+F(B)+2*TC+4*TL)*H/3; Writeln('Gia trị gan dung cua tich phan la',SUM); Readln; END. 4.6 Dãy quy tắc Với công thức hình thang độ chính xác tương ứng là 0(h2), nhưng còn phụ thuộc vào hàm đã cho, nên việc đánh giá sai số rất phức tạp. Vì vậy vấn đề là chia nhỏ [a ; b] đến mức nào ta có thể tin tưởng rằng đã đủ đảm bảo độ chính xác cho trước mà không cần chia nhỏ hơn nữa. Muốn vậy ta dựa vào dãy quy tắc như sau. b- a * Với các bước chia h j = , j = 1 , 2 , . . . ta có các điểm chia tương ứng là 2j j xk = a + kh j , k = 0 , . . . , 2 = n , x0 = a , xn = b và ký hiệu b * ST j » ò f(x)dx tính theo công thức hình thang. a b * SS j » ò f(x)dx tính theo công thức Simpson. a a- Dãy quy tắc hình thang * Tính h ST » 1 [f(x ) + f(x ) +2f(x )] 1 2 0 2 1 (4.6.1) ST » h2 {f(x ) + f(x ) +2[f(x ) + f(x ) + f(x )]} 2 2 0 4 1 2 3 (4.6.2) * Nếu ⎜ST2 – ST1 ⎜ ≤ ε (4.6.3) trong đó ε > 0, cho trước (sai số đối với tích phân) thì ta dừng ở ST2. * Ngược lại, nếu 71
⎜ST2 – ST1 ⎜ > ε (4.6.4) ta tính h ST ≈ 3 {f(x ) + f(x ) + 2[f(x ) + f(x ) + . . . + f(x )]} 3 2 0 8 1 2 7 (4.6.5) * Nếu ⎜ST3 – ST2 ⎜ ≤ ε (4.6.6) thì ta dừng ở ST3 . * Ngược lại, nếu ⎜ST3 – ST2 ⎜ > ε (4.6.7) ta tính h ST ≈ 4 {f(x ) + f(x ) + 2[f(x ) + f(x ) + . . . + f(x )]} 4 2 0 16 1 2 15 (4.6.8) . . . * Một cách tổng quát Ta tính với n = 2 j h ST ≈ j {f(x ) + f(x ) + 2[f(x ) + f(x ) + . . . + f(x )]} j 2 0 n 1 2 n - 1 (4.6.9) và dừng lại nếu thoả mãn một trong 2 điều kiện sau đây ⎜STj – STj -1 ⎜ ≤ ε (4.6.10) hoặc ì ï ST-> ST e íï jj 12 ï îï jM= (4.6.11) trong đó M là số bước lặp tối đa mà ta phải dừng, mặc dù đến bước lặp thứ M – 1 vẫn chưa đảm bảo sai số. b- Dãy quy tắc Simpson * Tính 4ST - ST SS =2 1 2 3 (4.6.12) SS = 4ST32 - ST 3 3 (4.6.13) * Nếu 72
⎜SS3 – SS2 ⎜ ≤ ε (4.6.14) thì ta dừng ở SS3 . * Ngược lại, nếu ⎜SS3 – SS2 ⎜ > ε (4.6.15) ta tính SS =4ST43 - ST 4 3 (4.6.16) . . . * Một cách tổng quát Với n = 2j ta tính SS = 4STjj - ST - 1 j 3 (4.6.17) và dừng lại nếu thoả mãn một trong 2 điều kiện sau đây ⎜SSj – SSj - 1 ⎜ ≤ ε (4.6.18) hoặc ì ï SS-> SS e íï jj 12 ï îï jM= (4.6.19) trong đó M là số bước lặp tối đa mà ta phải dừng, mặc dù đến bước lặp thứ M – 1 vẫn chưa đảm bảo sai số. Ví dụ 4.6 2 Với ε = 0,01 và M = 10 hãy tính gần đúng giá trị ò 3 5x2 + 1dx - 1 a- Theo dãy quy tắc hình thang. b- Theo dãy quy tắc Simpson. Giải a- Theo dãy quy tắc hình thang * h = 2(1) = 1,5 1 2 h ST » 1 [f(x ) + f(x ) + 2f(x )] 1 2 0 2 1 1,5 é 22 2ù = ê 335.(-++ 1) 13 5.(2) ++ 1 2 5.(0,5) + 1 ú » 5,397590 2 ë û 2(1) * h 2 = = 0,75 22 ST » h2 {f(x ) + f(x ) +2[f(x ) + f(x ) + f(x )]} 2 2 0 4 1 2 3 73
0,75 é 22 22 2 = ê335.(-++ 1) 1 5.(2) + 1 +2333 5.( - 0,25) ++ 1 5.(0,5) ++ 1 5.(1,25) + 1 » 5,069108 2 ë () ⎜ST2 – ST1 ⎜ = 0,328482 > 0,01 2(1) * h 3 = = 0,375 23 h ST » 3 {f(x ) + f(x ) + 2[f(x ) + f(x ) + . . . + f(x )]} 3 2 0 8 1 2 7 0,375 é 22 2 = ê335.(-++ 1)22 1 5.(2) + 1 +-233 5.( 0,625) ++-++++ 1 5.( 0,250) 1 3 5.(1,625) 1 » 4,999049 2 ë () ⎜ST3 – ST2 ⎜ = 0,070059 > 0,01 2(1) * h 4 = = 0,1875 . 24 h ST » 4 {f(x ) + f(x ) + 2[f(x ) + f(x ) + . . . + f(x )]} 4 2 0 16 1 2 17 0,1875 é 22 = ê335.(-++ 1)22 1 5.(2) + 1 +-233 5.( 0,8125) ++-+ 1 5.( 0,625) 1 2 ë ( ù ++. . .3 5.(1,8125)2 + 1 ú » )ûú 4,982439 ⎜ST4 – ST3 ⎜ = 0,016610 > 0,01 2(1) * h 5 = = 0,09375 25 h ST » 5 {f(x ) + f(x ) + 2[f(x ) + f(x ) + . . . + f(x )]} 5 2 0 32 1 2 31 0,09375 é = ê335.(-++ 1)22 1 5.(2) + 1 2 ë ù +-233 5.( 0,90625)22 ++-+++ 1 5.( 0,8125) 1 . . . 3 5.(1,90625) 2 + 1 ú » 4,978296 ()ûú ⎜ST5 – ST4 ⎜ = 0,004143 0,01 . SS = 4ST43 - ST ≈ 4.4,982439 - 4,999049 ≈ 4,976903 4 3 3 ⎜SS4 – SS3 ⎜= 0,001207 < 0,01 Ta dừng ở SS 4 » 4,976903 Program DHT_SS; Var ST,SS: Array[0 15] of real. I,J,K,M,N,S:Integer; 74
H,A B,SUM:Real; Function F(X:Real):Real; Begin F:=Exp(ln(5*X*X+1)/3); End; b 3 {Tính gần đúng gía trị ∫ 5x2 + 1dx .} a BEGIN Write('Can duoi A=');Readln(A); Write('Can tren B=');Readln(B); Write('So lan lap,M=');Readln(M); H:=B-A;N:=1;ST[0]:=H*(F(A)+F(B))/2; For J:=1 to M do Begin S:=N;N:=2*N;SUM:=0;H:=H/2; For K:=1 to S do SUM:=SUM+F(A+H*(2*K-1));ST[J]:=ST[J– 1]/2+H*SUM; End; Writeln('Gia tri gan dung cua tich phan'); Writeln(' Day hinh thang Day Simpson'); For J:=1 to M do Writeln('ST[',J,']=',ST[J],' SS[',J,']=',(4*ST[J] –ST[J–1])/3); Readln; END. Bài tập x−1 1- Cho f(x) = sin . x+1 a- Tính f(0,05) , f(0,06) , . . . , f(0,09) , f(0,10) , . . . , f(0,14) , f(0,15) . b- Bằng phương pháp Richardson hãy tính gần đúng f'(0,1) , ε = 0,001 , M = 10. Đs: 1,129579 f'(0,1) , ε = 0,001 , M = 20. Đs: 1,129579 f'(0,1) , ε = 0,0001 , M = 10. Đs: 1,129818 f'(0,1) , ε = 0,0001 , M = 5. Đs: Chưa hội tụ. f'(– 1,3) , ε = 0,01 , M = 5. Đs: 4,136889 f'(– 1,3) , ε = 0,001 , M = 10. Đs: 4,138170 2- Tính f’(0,06) theo công thức Richardson theo bảng số liệu xk 0,01 0,02 0,03 0,04 0,05 0,06 0,07 0,08 0,09 0,10 fk 1,12 1,18 1,30 1,50 1,75 1,90 2,10 2,15 2,25 2,30 Đs : 17,977778 3- Dùng công thức Gregory-Newton tiến a- Tính f'(1,01) với bảng số liệu xk 1,01 1,02 1,03 1,04 1,05 fk 4,321 8,600 12,837 17,033 21,189 Đs: 430,03333 b- Tính f'(10,6) với bảng số liệu 75
xk 10,6 10,7 10,8 10,9 fk 25,306 29,384 33,424 37,426 Đs: 40,97. 4- Cho bảng số liệu x 1,01 1,02 1,03 1,04 1,05 fk 4,321 8,600 12,837 17,033 21,189 Tính f'(1,028) bằng công thức Gauss, Lagrange,Newton. 5- Cho bảng số liệu về độ cao của ngọn núi (đơn vị là m) : xk ( khoảng cách từ chân núi ) 10 20 30 40 50 60 fk ( độ cao của sườn núi ) 5,2 8,1 15,3 30,4 60,5 82,6 Tìm độ dốc của sườn núi tại điểm cách chân núi 11m ( f’(11)?). Đs : – 0,568763 2,2 cos(2x− 1) 6- Tính gần đúng bằng công thức Simpson với 2N = 10 : ∫ dx . 1 x Đs : – 0,270676 3,2 cos2x 7- Tính gần đúng bằng công thức Simpson với 2N = 10 : ∫ dx 2 21x+ Đs : – 0,15416 8- Tính gần đúng bằng qui tắc Hình thang, Simpson. Qui tắc Hình thang 2 2 a- ∫ exdxx sin , N = 8 Đs: 14,598274 N = 10 1 Đs: 14,519171 3 2 ∫ exdxx sin , N = 15 Đs: 457,497100 0 1 b- ∫ exdx− x cos 2 , N = 5 Đs: 0,590651 N = 10 0 Đs: 0,590202 Qui tắc Simpson 2 2 a- ∫ exdxx sin , 2N = 8 Đs: 14,381322 2N = 10 1 Đs: 14,379375 1 b- ∫ exdx− x cos 2 , 2N =10 Đs: 0,590053 2N = 20 0 Đs: 0,590051 9- Tính gần đúng diện tích hình phẳng sau đây bằng công thức hình thang và Simpson y 0,4 0,2 0 x 76
– 0,2 – 0,4 xk 3,5 2 4 1,5 4,1 1,2 4,2 1 4 fk – 0,4 – 0,3 – 0,3 – 0,2 – 0,2 – 0,1 – 0,1 0 0 1,5 3,7 1,7 3,5 1,9 3 2 0,1 0,1 0,2 0,2 0,3 0,3 0,4 Đs: S Hìnhthang = S1Hìnhthang + S2 Hìnhthang = 1,57 S Simpson = S1Simpson +S2 Simpson = 1,466667 10- Khi đo một khu đất (đơn vị là 100m) ta có bảng số liệu sau đây : xk fk xk fk xk fk f – 0,15 – 0,6 – 0,45 – 0,2 – 0,35 0,3 0,35 – 0,6 0,4 – 0,2 0,35 0,3 A 0,6 – 0,2 – 0,5 – 0,47 – 0,1 – 0,25 0,4 0,3 – 0,5 0,45 – 0,1 0,25 0,4 – 0,25 – 0,4 – 0,5 0 – 0,15 0,5 – 0,5 E 0 B 0,5 x 0,25 – 0,4 0,5 0 0,15 0,5 – 0,35 – 0,3 – 0,47 0,1 0 0,6 0,3 – 0,3 0,47 0,1 – 0,45 0,2 0,45 0,2 Tính gần đúng diện tích S ABCDEA của khu đất trên bằng công thức Simpson . Đs : 0,784667(10000m2) D – 0,6 C 11- Tính gần đúng bằng dãy qui tắc Hình thang, Simpson với ε = 0,01 và M = 15. 2 2 a- ∫ exdxx sin 1 1 b- ∫ exdx− x cos 2 0 Đs: a- Dãy qui tắc Hình thang Dãy qui tắc Simson N = 8 TS 1 = 18,215313 SS 1 = 14,593737 TS 2 = 15,499131 SS 2 = 14,464655 TS 3 = 14,723274 SS 3 = 14,353214 TS 4 = 14,445729 SS 4 = 14,415516 77
TS 5 = 14,423069 SS 5 = 14,388420 TS 6 = 14,397082 SS 6 = 14,383209 TS 7 = 14,386679 SS 7 = 14,380607 TS 8 = 14,382125 b- Dãy qui tắc Hình thang Dãy qui tắc Simpson N = 6 TS 1 = 1,093529 SS 1 = 0,756799 TS 2 = 0,840981 SS 2 = 0,673389 TS 3 = 0,715287 SS 3 = 0,631718 TS 4 = 0,652610 SS 4 = 0,610884 TS 5 = 0,621316 SS 5 = 0,600468 TS 6 = 0,605 Kiểm tra nhận thức 1* Nêu ưu điểm, nhược điểm của từng phương pháp tỷ sai phân và Richardson. 2* Viết chương trình PASCAL cho các phương pháp tỷ sai phân. 3* Viết chương trình PASCAL cho các phương pháp nội suy với các mốc cách đều và bất kỳ. 4* Nêu những điểm giống nhau, khác nhau của hai phương pháp hình thang, Simpson. 5* Hãy chạy chương trình PASCAL đã có để giải các ví dụ. Nếu có điều gì chưa thoả mãn thì chủ động viết lại chương trình (có thể bằng ngôn ngữ lập trình bậc cao khác) và giải lại các ví dụ, các bài tập. 78
Chương 5: XẤP XỈ NGHIỆM PHƯƠNG TRÌNH VI PHÂN VÀ PHƯƠNG TRÌNH ĐẠO HÀM RIÊNG Chúng ta biết rằng lớp các phương trình vi phân và phương trình đạo hàm riêng giải được nghiệm dưới dạng giải tích là rất hạn chế, mà ngay trong trường hợp này công thức nghiệm cũng rất phức tạp. Vì vậy trong chương này ta nêu ra một số phương pháp tìm giá trị gần đúng của nghiệm. 5.1 Xấp xỉ nghiệm phương trinh vi phân cấp một Tại x* cho trước, tính gần đúng nghiệm y(x*) của phương trình y' = f(x,y) (5.1.1) thoả mãn điều kiện y(x0) = y0 (5.1.2) a- Phương pháp Euler yx()()+- h yx Trên cơ sở định nghĩa đạo hàm y' (x) = lim ta có y(x+h) » y(x) + h® 0 h h.y’(x). * Đặt h = x*- x0 ⇔ n = x*- x0 n h (5.1.3) trong đó n là số khoảng chia nhỏ, h là độ dài các đoạn chia nhỏ. * Các điểm chia tương ứng xk = x0 + kh , k = 0 , 1 , 2 , . . . , n ta có yk = y(xk). * Ta tính yk+1 » yk + hf(xk , yk) (5.1.4) * Lần lượt ta tìm được y1 , y2 , . . . yn » y(x*) (5.1.5) * Sai số là hQ 2 E = yx()-£ y eRxk - x0 -1 º O(h ) k kk2R () (5.1.6) ¶ f trong đó ⎜y” ⎜≤ Q , £ R ¶ y Bảng Euler xk = x0 + kh yk Ak = hf(xk , yk) yk + 1 = yk + Ak x0 y0 A0 = hf(x0 , y0) y1 = y0 + A0 x1 = x0 + h y1 A1 = hf(x1 , y1) y2 = y1 + A1 x2 = x0 + 2h y2 A2 = hf(x2 , y2) y3 = y2 + A2 79
. . . . . . . . . . . . . . . . . . . . xn - 1 = x0 + (n – 1)h yn - 1 An - 1 = hf(xn – 1 , yn - 1) yn = yn - 1 + An - 1 b- Phương pháp Runge-Kutta bậc hai * Đặt h = x*- x0 ⇔ n = x*- x0 n h trong đó n là số khoảng chia nhỏ, h là độ dài các đoạn chia nhỏ. * Các điểm chia tương ứng xk = x0 + kh , k = 0 , 1 , 2 , . . . , n Û yk = y(xk) * Tính y » y + A kk+ C k + 1 k 2 (5.1.7) trong đó Ak = hf(xk , yk) , Bk = yk + Ak , Ck = hf(xk + 1 , Bk) (5.1.8) * Lần lượt ta tìm được y1 , y2 , . . . yn » y(x*) (5.1.9) với sai số là O(h3). Bảng Runge-Kutta bậc hai xk = x0 + kh yk Ak = hf(xk , yk) Bk = yk + Ak Ck = hf(xk + 1 , Bk) yk + 1 = yk + Akk+ C 2 x0 y0 A0 = hf(x0 , y0) B0 = y0 + A0 C0 = hf(x1 , B0) y1 = y0 + A00+ C 2 x1 y1 A1 = hf(x1 , y1) B1 = y1 + A1 C1 = hf(x2 , B1) y2 = y1 + A11+ C 2 . . . . . . . . . . . . . . . . . . . . xn - 1 yn - 1 A n - 1 = hf(x n - 1 , yn - 1) B n – 1 = y n – 1 + A n – 1 C n - 1 = hf(x n , B n - 1) yn = y n - 1 AC+ + nn 11 2 c- Phương pháp Runge-Kutta bậc bốn A2B2CD+++ y ≈ y + kkkk k+1 k 6 (5.1.10) trong đó ⎛⎞hAk ⎛⎞hBk Ak = hf(xk,yk) Bk = h f,⎜⎟xyk ++k Ck = h f,⎜⎟xykk++ Dk = ⎝⎠22 ⎝⎠22 hf(xk+1,yk+Ck) với sai số là 0(h5). Tương tự ta lập Bảng Runge-Kutta bậc bốn có 7 cột. Ví dụ 5.1 Tìm giá trị gần đúng của y(0,4) từ phương trình y' = cos(x + 2y) – xy2 80
thoả mãn điều kiện y(0) = 1 với h = 0,1. a- Bằng phương pháp Euler b- Bằng phương pháp Runge-Kutta bậc hai c- Bằng phương pháp Runge-Kutta bậc bốn Giải n = x*- x0 = 0,4- 0 = 4 h 0,1 a- Bảng Euler 2 xk = 0 + 0,1k yk Ak = 0.1[cos(xk + 2yk) – xky k] yk + 1 = yk + Ak x0 = 0 1 – 0,041615 0,958385 x1 = 0,1 0,958385 – 0,052319 0,906067 x2 = 0,2 0,906067 – 0,059134 0,846933 x3 = 0,3 0,846933 – 0,062575 0,784358 Vậy y(0,4) » 0,784358 b- Bảng Runge-Kutta bậc hai 2 xk = 0+0,1k yk Ak = 0,1[cos(xk + 2yk) – xky k] Bk = yk + Ak x0 = 0 1 – 0,041615 0,958385 x1 = 0,1 0,953033 – 0,051248 0,901785 x2 = 0,2 0,898307 – 0,057446 0,840861 x3 = 0,3 0,839006 – 0,060723 0,778283 C = 0,1[cos(x + 2B ) – x B 2] y = y + Akk+ C k k + 1 k k + 1 k k + 1 k 2 – 0,052319 0,953033 – 0,058203 0,898307 – 0,061157 0,839006 – 0,061856 0,777716 Vậy y(0,4) » 0,777716 c Giá trị gần đúng của nghiệm theo Runge-Kutta bậc bốn y(0) = 1 y(0,1) ≈ 0,953242 y(0,2) ≈ 0,898575 y(0,3) ≈ 0,839216 y(0,4) ≈ 0,777798 Program RKUTTA4; Var X,Y:Array[0 10]; K,M:Integer; A,B,C,D,E:Real; Function F(X,Y:Real):Real; Begin F:=cos(X+2*Y) – X*Y*Y; End; 2 {Giải phương trình y' = cos(x+2y) – xy thỏa mãn điều kiện y(x0) = y0 .} BEGIN Writeln('Dieu kien ban dau '); Write('Diem dau X[0]=');Readln(X[0]\); Write('Y(‘,X[0],’)=');Readln(Y[0]); Write('Diem cuoi X[M]=');Readln(E); 81
Write('So khoang M=');Readln(M); H:=(E-X[0])/M; For K:=0 to M-1 do Begin A:=H*F(X[K],Y[K]);B:=H*F(X[K]+H/2,Y[K]+A/2); C:=H*F(X[K]+H/2,Y[K]+B/2);D:=H*F(X[K]+H,Y[K]+C); Y[K+1}:=Y[K]+(A+2*B+2*C+D)/6;X[K+1]:=X[K]+H; End; Writeln('Gia tri gan dung’); For K:=1 to M do Write('Y[',X[K],’]=',Y[K]’); Readln; END. 5.2 Xấp xỉ nghiệm hệ phương trình vi phân cấp một Tính gần đúng nghiệm x(t*), y(t*) tại t* cho trước của hệ phương trình ì ï x¢= ftxy(, , ) í ï îï ygtxy¢= (, , ) (5.2.1) thoả mãn điều kiện x(t0) = x0 , y(t0) = y0 (5.2.2) a- Phương pháp Euler * Đặt h = tt*- 0 ⇔ n = tt*- 0 n h (5.2.3) trong đó n là số khoảng chia nhỏ, h là độ dài đoạn chia nhỏ. * Các điểm chia tương ứng tk = t0 + kh, k = 0,1,2, . . . , n . Û xk = x(tk) , yk = y(tk) . * Tính ì ï xk1+ »+x khf (t,x,y) k k k íï ï yy(t,x,y) »+hg îï k1+ k k k k (5.2.4) Như vậy lần lượt ta tìm được x1 , y1 , x2 ,y2 , . . . xn » x(t*) , yn » y(t*) (5.2.5) * Sai số theo x là hQ 2 x Rtx k - t0 Exk = xt()kk-£ x() e -1 º O(h ) 2Rx (5.2.6) ¶ f trong đó ⎜x” ⎜≤ Qx , £ Rx . ¶ x Sai số theo y là hQy 2 Rtyk- t0 Eyk = yt()kk-£ y() e -1 º O(h ) 2Ry (5.2.7) 82
¶ g trong đó ⎜y” ⎜≤ Qy , £ Ry . ¶ y b- Phương pháp Runge-Kutta bậc bốn * Đặt h = tt*- 0 ⇔ n = tt*- 0 n h trong đó n là số khoảng chia nhỏ, h là độ dài đoạn chia nhỏ. * Các điểm chia tương ứng tk = t0 + kh, k = 0,1,2, . . . , n Û xk = x(tk) , yk = y(tk) * Với mỗi k ta tính A1 = h.f(tk , xk , yk) B1 = h.g(tk , xk , yk) æö ç hA11 B÷ A2 = h. ftç kk++,, x y k +÷ B2 = èøç 22 2÷ æöhA11 B÷ h. gtç kk++,, x y k +÷ èøç 22 2÷ æö ç hA22 B÷ A3 = h. ftç kk++,, x y k +÷ B3 = èøç 22 2÷ æöhA22 B÷ h. gtç kk++,, x y k +÷ èøç 22 2÷ A4 = h.f(tk + h , xk + A3 , yk + B3) B4 = h.g(tk + h , xk + A3 , yk + B3) x » x + A1234 +++ 2A 2A A y » y + k + 1 k 6 k + 1 k B1234 +++ 2B 2B B 6 Như vậy lần lượt ta tìm được x1 , y1 , x2 ,y2 , . . . xn » x(t*) , yn » y(t*) Ví dụ 5.2 ⎪⎧x′=+xy2 Giải hệ ⎨ thoả mãn x(0) = 6 , y(0) = 4 ⎩⎪yxy′=+32 tìm xấp xỉ của x(0,1) và y(0,1) với h = 0,02. a- Bằng phương pháp Euler. b- Bằng phương pháp Runge-Kutta bậc bốn. Giải n = tt*- 0 = 0,1- 0 = 5 h 0,02 a- Phương pháp Euler x1 = x(0,02) = x0 + hf(t0 , x0 , y0) = 6 + 0,02(6 + 2.4) = 6,28 y1 = y(0,02) = y0 + hg(t0 , x0 , y0) = 4 + 0,02(3.6 + 2.4) = 4,52 x2 = x(0,04) = x1 + hf(t1 , x1 , y1) = 6,28 + 0,02(6,28 + 2.4,52) = 6,5864 y2 = y(0,04) = y1 + hg(t1 , x1 , y1) = 4,52 + 0,02(3.6,28 + 2.4,52) = 5,0776 x3 = x(0,06) = x2 + hf(t2 , x2 , y2) = 6,5864 + 0,02(6,5864 + 2.5,0776) = 6,931232 y3 = y(0,06) = y2 + hg(t2 , x2 , y2) = 5,0776 + 0,02(3.6,5864 + 2.5,0776) = 5,675888 x4 = x(0,08) = x3 + hf(t3 , x3 , y3) = 6,921232 + 0,02(6,921232 + 2.5,675888) » 7,286692 y4 = y(0,08) = y3 + hg(t3 , x3 , y3) = 5,675888 + 0,02(3.6,921232 + 2.5,675888) » 6,318197 x5 = x(0,10) = x4 + hf(t4 , x4 , y4) = 7,286692 + 0,02(7,286692 + 2.6,318197) 83
» 7,685115 y5 = y(0,10) = y4 + hg(t4 , x4 , y4) = 6,318197 + 0,02(3.7,286692 + 2.6,318197) » 7,008127 Vậy x(0,1) » 7,685115 , y(0,1) » 7,008127 b- Phương pháp Runge-Kutta bậc bốn Trước tiên ta tìm x1 = x(0,02) , y1 = y(0,02): A1 = hf(t0 , x0 , y0) = 0,02(6 + 2.4) = 0,28 B1 = hg(t0 , x0 , y0) = 0,02(3.6 + 2.4) = 0,52 æöé æöù ç hA11 B÷ ê 0,28ç 0,52 ÷ú A2 = h ftç 00++,, x y 0 +÷ = 0,02ê 6+++ 2ç 4 ÷ú = 0,2932 èø22 2 ë 22èøû æöé æöæöù ç hA11 B÷ ê çç0,28÷÷ 0,52 ú B2 = h gtç 00++,, x y 0 +÷ = 0,02ê 3çç 6+++÷÷ 2 4 ú = 0,5388 èø222 ë èøèø22û æöhA B é 0,2932æö 0,5388 ù ç 22÷ 0,02ê 6 2ç 4 ÷ú A3 = h ftç 00++,, x y0 +÷ = ê +++ç ÷ú = 0,293708 èø22 2 ë 22èøû æöhA B é æöæö0,2932 0,5388 ù ç 22÷ 0,02ê 3çç 6 ÷÷2 4 ú B3 = h gtç 00++,, x y 0 +÷ = ê çç+++÷÷ú = 0,541968 èø22 2 ë èøèø22û A4 = hf(t0 + h , x0 + A3 , y0 + B3) = 0,02[(6 + 0,293708) + 2.(4 + 0,541968)] = 0,307553 B4 = hg(t0 + h , x0 + A3 , y0 + B3) = 0,02[3.(6 + 0,293708) + 2.(4 + 0,541968)] = 0,559205 x = x + A1234 +++ 2A 2A A » 6,293546 1 0 6 y = y + B1234 +++ 2B 2B B » 4,539325 1 0 6 Tương tự ta tìm được: x2 = x(0,04) » 6,615622 , y2 = y(0,04) » 5,119486 x3 = x(0,06) » 6,968525 , y3 = y(0,06) » 5,743965 x4 = x(0,08) » 7,547432 , y4 = y(0,08) » 6,416533 x5 = x(0,10) » 7,776973 , y5 = y(0,10) » 7,141272 Vậy x(0,1) » 7,776973, y(0,1) » 7,141272 Program EULER-HEPTVP; Var K,N:Integer; T,X,Y:Array[0 15] of real; Tc:Real; Function F(T,X,Y:Real):Real; Begin F:= X+2*Y; End; Function G(T,X,Y:Real):Real; Begin G:=3*X+2*Y: End; 84
⎧⎪x′=+xy2 { Giải hệ phương trình vi phân ⎨ } ⎩⎪yxy′=32+ BEGIN Writeln(' Dieu kien ban dau'); Write(' T[0]=');Readln(T[0]); Write(' X[‘,T[0]:4:2,’]=');Readln(X[0]); Write(' Y[‘,T[0]:4:2,’]=');Readln(Y[0]); Write(' Can tinh gia tri xap xi cua X,Y tai Tc =');Readln(Tc); Write(' So diem chia N = ');Readln(N); H:= (Tc–T[0])/N; For K:=1 to N do T[K]:= T[K–1] + H; Writeln(’Gia tri gan dung tai cac diem chia’); For K:= 0 to N–1 do Begin X[K+1] := X[K] + F(T[K],X[K],Y[K]); Y[K+1] := Y[K] + G(T[K],X[K],Y[K]); Writeln('X[‘,T[K+1]:4:2,’]=', X[K+1],'Y[‘,T[K+1]:4:2,’]=', Y[K+1]); End; Readln; END. 5.3 Xấp xỉ nghiệm phương trình vi phân cấp 2 Tính gần đúng nghiệm x (t*) tại t* cho trước của phương trình x''(t) = f(t , x(t) , x'(t)) (5.3.1) thoả mãn điều kiện x(t0) = x0 , x'(t0) = x'0 (5.3.2) Ta có thể đưa (5.3.1) về hệ phương trình như sau ïì x'(t) = y(t) íï ï îï y'(t) = f(t,x(t), y(t)) (5.3.3) x(t0) = x0 , y(t0) = x'0 (5.3.4) Tìm xấp xỉ x(t*) (y(t*)). Ví dụ 5.3 Giải phương trình x''(t) + 4x'(t) + 5x(t) = 0 thoả mãn điều kiện x(0) = 3 , x'(0) = – 5 tìm giá trị gần đúng của x(0,4) với h = 0,1 Giải ïì x'(t) = y(t) íï * Ta có hệ tương ứng ï với điều kiện x(0) = 3 , y(0) = – 5. îï y'(t) =- 5x(t) - 4y(t) * n = tt*- 0 = 0,4- 0 = 4 h 0,1 Chẳng hạn bằng phương pháp Euler ta tìm được x(0,1) ≈ x1 = x0 + hy0 = 3 + 0,1(–5) = 2,5 85
y(0,1) ≈ y1 = y0 + hf(t0 , x0 , y0) = – 5 + 0,1[–5.3 – 4(– 5)] = – 4,5 x(0,2) ≈ x2 = x1 + hy1 = 2,5 + 0,1(– 4,5) = 2,05 y(0,2) ≈ y2 = y1 + hf(t1 , x1 , y1) = – 4,5 + 0,1[–5.2,5 – 4(– 4,5)] = – 3,95 x(0,3) ≈ x3 = x2 + hy2 = 2,05 + 0,1(– 3,95) = 1,655 y(0,3) ≈ y3 = y2 + hf(t2 , x2 , y2) = – 3,95 + 0,1[–5.2,05 – 4(– 3,95)] = – 3,395 x(0,4) ≈ x4 = x3 + hy3 = 1,655 + 0,1(–3,395) = 1,3155 Vậy x(0,4) » 1,3155 5.4 Xấp xỉ nghiệm phương trình đạo hàm riêng Trong thực tế ta thường gặp các bài toán sau đây * Các bài toán về trường điện từ, động lực học chất lỏng, trạng thái cân bằng nhiệt, khuyếch tán . . . dẫn ¶¶22uu đến phương trình Laplace + = 0. ¶¶x22y * Các bài toán về truyền nhiệt, thẩm thấu chất lỏng, chất khí dẫn đến phương trình 2 ¶¶uu2 = a . ¶ t ¶ x2 * Các bài toán về dao động dây, dao động điện . . . dẫn đến phương trình 22 ¶¶uu2 = a . ¶¶tx22 Chúng ta sẽ sử dụng phương pháp sai phân để giải gần đúng các phương trình trên đây. a- Giải gần đúng phương trình ¶¶22uu + = 0 ¶¶x22y (5.4.1) trong miền chữ nhật a ≤ x ≤ b , c ≤ y ≤ d (5.4.2) thoả mãn điều kiện biên u xacy=££, d= v1(y) (5.4.3) u xbcyd=££, = v2(y) (5.4.4) u axbyc££, == v3(y) (5.4.5) u axbyd££, = = v4(y) (5.4.6) v1(c) = v3(a) , v1(d) = v4(a) (5.4.7) v2(c) = v3(b) , v2(d) = v4(b) (5.4.8) Ta chia đều [a ; b] bởi các điểm x = a + ih , h = ba- , i = 0 , 1 , . . . , n. i n 86