Chương 6: Giải gần đúng phương trình vi phân

Chương 6GIẢI GẦN ĐÚNGPHƯƠNG TRÌNH VI PHÂNI. GIẢI GẦN ĐÚNG PTVP CẤP 1 : Xét bài toán Cauchy : tìm nghiệm y=y(x) của phương trình vi phân với giá trị ban đầu y0	y’ = f(x, y), x  [a,b]	y(a) = y0Các phương pháp giải gần đúng : Công thức Euler Công thức Euler cải tiến Công thức Runge-Kutta1. Công thức Euler : Để tìm nghiệm gần đúng của bài toán Cauchy ta chia đoạn [a,b] thành n đoạn nhỏ bằng nhau với bước h = (b-a)/nxo= a, x1 = x0 +h, ... , xk = x0 + kh, ... , xn = bNghiệm gần đúng của bài toán là dãy {yk} gồm các giá trị gần đúng của hàm tại xk 	Ta có yk  y(xk) , k =0, nGiả sử bài toán có nghiệm duy nhất y(x) có đạo hàm đến cấp 2 liên tục trên [a,b]. Khai triển Taylor ta có y(xk+1) = y(xk) + (xk+1-xk) y’(xk) + (xk+1-xk)2 y’’(k)/2 với k  (xk, xk+1)Công thức Euler : 	yk+1 = yk + h f(xk, yk) , k = 0, n-1	với h = xk+1 - xkVí dụ : Dùng công thức Euler tìm nghiệm gần đúng của bài toán Cauchy 	y’ = y – x2 +1, 	0≤x≤1	y(0) = 0.5với n = 5 Tính sai số biết nghiệm chính xác là : 	y(x) = (x+1)2 – 0.5exgiảita có h = 0.2x0 = 0, x1 = 0.2, x2 = 0.4, x3 = 0.6, x4 = 0.8, x5 = 1Công thức Euler y0 = 0.5 yk+1 = yk + h f(xk, yk) = yk + 0.2 (yk - xk2 +1) kxkyky(xk)|y(xk) - yk |000.50.5010.20.80.82929860.029298620.41.1521.21408770.062087730.61.55041.64894060.098540640.81.988482.12722950.1387495512.4581762.64085910.1826831A = 0B = 0.5B = B + 0.2(B – A2 + 1) : A=A+0.2:(A+1)2-0.5eA:Ans-B* Nhận xét : công thức Euler đơn gian, nhưng sai số còn lớn nên ít được sử dụng 2. Công thức Euler cải tiến : 	yk+1 = yk + (k1+k2)/2	k = 0,1, ..., n-1	k1 = hf(xk, yk), 	k2 = hf(xk+h, yk + k1)	với h = xk+1 - xkVí dụ : Dùng công thức Euler cải tiến tìm nghiệm gần đúng của bài toán Cauchy 	y’ = y – x2 +1, 	0≤x≤1	y(0) = 0.5với n = 5 Tính sai số biết nghiệm chính xác là : 	y(x) = (x+1)2 – 0.5exgiảita có h = 0.2x0 = 0, x1 = 0.2, x2 = 0.4, x3 = 0.6, x4 = 0.8, x5 = 1Công thức Euler cải tiến	yo = 0.5	yk+1 = yk + (k1 +k2) /2	k1= 0.2(yk - xk2 +1) 	k2 = 0.2(yk + k1 – (xk+0.2)2 +1)kxkyky(xk)|y(xk) - yk |000.50.5010.20.8260.82929860.003320.41.206921.21408770.007230.61.63724241.64894060.011740.82.11023572.12722950.0170512.61768762.64085910.0232A = 0 (xk)B = 0.5 (yk)C = 0.2(B – A2 + 1) : D = 0.2(B + C - (A+0.2)2 + 1):B=B + (C+D)/2:A=A+0.2:(A+1)2-0.5eA:Ans-B3. Công thức Runge Kutta bậc 4 : Ví dụ : Xét bài toán Cauchy 	y’ = 2.7xy + cos (x+2.7y), 	1.2≤x	y(1.2) = 5.4Dùng công thức Runge-Kutta tính gần đúng y(1.5) với bước h = 0.3xo = 1.2, yo = 5.4y1 = y0 + (K1+ 2K2+ 2K3+ K4) /6Công thức Runge-Kutta bậc 4giảiK1= 0.3(2.7xoyo + cos(xo+2.7yo)) K2= 0.3(2.7(xo+0.3/2)(yo+K1/2) +cos(xo+0.3/2 +2.7(yo+K1/2))K3= 0.3(2.7(xo+0.3/2)(yo+K2/2) +cos(xo+0.3/2 +2.7(yo+K2/2)) K4= 0.3(2.7(xo+0.3)(yo+K3) +cos(xo+0.3 +2.7(yo+K3) Bấm máy ta được K1 = 4.949578057	K2 = 8.367054617K3 = 10.33000627	K4 = 19.41193853y(1.5) = 15.69260639  15.6926Ví dụ : Dùng công thức Runge-Kutta tìm nghiệm gần đúng của bài toán Cauchy 	y’ = y – x2 +1, 	0≤x≤1	y(0) = 0.5với n = 5 Tính sai số biết nghiệm chính xác là : 	y(x) = (x+1)2 – 0.5exgiảita có h = 0.2x0 = 0, x1 = 0.2, x2 = 0.4, x3 = 0.6, x4 = 0.8, x5 = 1A = 0 (xk)B = 0.5 (yk)C = 0.2(B – A2 + 1) : D = 0.2(B + C/2 - (A+0.1)2 + 1):E = 0.2(B + D/2 - (A+0.1)2 + 1):F = 0.2(B + E - (A+0.2)2 + 1):B =B + (C+2D+2E+F)/6:A =A+0.2:(A+1)2-0.5eA:Ans-Byk+1 = yk + (K1+ 2K2+ 2K3+ K4) /6Công thức Runge-Kutta bậc 4K2 = 0.2 [yk + 0.1(yk - xk2 +1) –(xk+0.1)2 +1 ] = 0.2(1.1 yk – 1.1xk2 – 0.2xk + 1.09)K1= 0.2(yk - xk2 +1) K3 = 0.2[ yk + 0.1(1.1yk – 1.1xk2 – 0.2xk + 1.09) – (xk+0.1)2 +1 ] = 0.2(1.11yk – 1.11xk2 – 0.22xk + 1.099)K4 = 0.2[ yk+0.2(1.11yk–1.11xk2–0.22xk+1.099) – (xk+0.2)2 +1 ] = 0.2(1.222yk–1.222xk2–0.444xk+1.1798)y0 = 0.5yk+1 = yk+0.2(6.642yk–6.642xk2–1.284xk+6.5578)/6kxkyky(xk)|y(xk) - yk |000.50.5010.20.82929330.82929860.000005320.41.21407621.21408770.000011530.61.64892201.64894060.000018640.82.12720272.12722950.0000269512.64082272.64085910.0000364II. GIẢI GẦN ĐÚNG HỆ PTVP : Xét hệ phương trình vi phân cấp 1y’1 = f1(x, y1, y2, ..., ym)y’2 = f2(x, y1, y2, ..., ym)	. . .y’m = fm(x, y1, y2, ..., ym)với a≤ x ≤ b và thỏa điều kiện ban đầu	y1(a) = 1, y2(a) = 2, .... , ym(a) = mNghiệm y = (y1, y2, , ym)Để tìm nghiệm gần đúng, ta chia đoạn [a,b] thành n đoạn nhỏ bằng nhau với bước h = (b-a)/n và các điểm chiaxo= a, x1 = x0 +h, ... , xk = x0 + kh, ... , xn = bCông thức Euler : 	yi k+1 = yi k + h fi(xk, y1 k,  , ym k)	i=1..m; k = 0.. n-1Nghiệm gần đúng là dãy { yk=(y1 k, y2 k, , ym k)}với yi k  yi(xk)Công thức Euler cải tiến : 	yi k+1 = yi k + (K1 i + K2 i) / 2	K1 i = h fi(xk, y1 k,  , ym k)	K2 i = h fi(xk+h, y1 k+K1 1,  , ym k+K1 m)	i=1,m; k = 0, n-1Công thức Runge-Kutta bậc 4 : 	yi k+1 = yi k + (K1 i+2K2 i+2K3 i+K4 i) / 6	K1 i = h fi(xk, y1 k,  , ym k)	K2 i = h fi(xk+h/2, y1 k+K11/2,  , ym k+K1 m/2)	K3 i = h fi(xk+h/2, y1 k+K21/2,  , ym k+K2 m/2)	K4 i = h fi(xk+h, y1 k+K31,  , ym k+K3 m)	i=1,m; k = 0, n-1Ví dụ : Sử dụng công thức Euler giải gần đúng hệ pt vi phân	y’1 = 3y1 + 2y2 – (2x2 +1)e2x	y’2 = 4y1 + y2 + (x2 +2x –4) e2xvới 0 ≤x≤0.5điều kiện ban đầu y1(0)=y2(0)=1bước h = 0.1So sánh với nghiệm chính xác	y1(x) = 1/3e5x –1/3e-x+e2x	y2(x) = 1/3e5x +2/3e-x+x2e2xCông thức Eulery1 0 = 1y1 k+1 = y1 k + h (3y1k + 2y2 k – (2xk2 +1)e2xk)y2 0 = 1y2 k+1 = y2 k + h (4y1k + y2 k + (xk2 +2xk –4) e2xk)xky1ky1(xk)y2ky2(xk)011110.11.41.46941.11.16500.21.91542.12501.30711.51160.32.59033.06911.67292.15180.43.48704.46512.27323.26600.54.69406.57693.21875.1448A=0 (x)B=1 (y1k)C=1 (y2k)D=B + 0.1 (3B + 2C – (2A2 +1)e2A):C=C + 0.1 (4B + C + (A2 +2A –4) e2A):B=D:A=A+0.1A=0e5A/3–e-A/3+e2A:e5A/3+2/3e-A/3+A2e2A:A=A+0.1III. GIẢI GẦN ĐÚNG PTVP CẤP CAO: Xét phương trình vi phân bậc m 	y(m)(x) = f(x, y, y’, ... , y(m-1)), a≤x≤bvới điều kiện ban đầu 	y(a) = 1, y’(a) = 2, .... , y(m-1)(a) = mĐặt	y1 = y, y2 = y’, y3 = y”, ... , ym = y(m-1)Ta chuyển phương trình vi phân bậc m về hệ m phương trình vi phân cấp 1 với điều kiện ban đầu 	y1(a) = 1, y2(a) = 2, .... , ym(a) = m,y’1 = y2y’2 = y3. . .y’m-1 = ymy’m = f(x, y, y’, ... , y(m-1))Ví dụ : Sử dụng công thức Euler giải gần đúng pt vi phân cấp 2	y “ – 2 y’ + 2y = sinx e2x , 0≤x≤0.5 điều kiện ban đầu 	y(0) = -0.4, y’(0) = -0.6với bước h = 0.1 So sánh với nghiệm chính xác biết nghiệm CX 	y1(x) = 0.2e2x (sinx – 2cosx)	y2(x) = 0.2e2x(4sinx - 3cosx)=y’đặt y1 = y, y2 = y’ chuyển pt về hệ	y’1 = y2	y’2 = sinx e2x– 2 y1 + 2y2điều kiện y1(0) = -0.4, y2(0) = -0.6Công thức Euler	y1 0 = -0.4	y1 k+1 = y1 k + 0.1 y2k 	y2 0 = -0.6	y2 k+1 = y2 k + 0.1 (sinxke2xk - 2y1k +2y2 k)xky1 ky1(xk)y2 ky2(xk)=y’(xk)0-0.4-0.4-0.6-0.60.1-0.46-0.4617-0.64-0.63160.2-0.524-0.5256-0.6638-0.64010.3-0.5904-0.5886-0.6621-0.61370.4-0.6566-0.6466-0.6226-0.53660.5-0.7189-0.6936-0.5292-0.3887A=0B=-0.4C=-0.6D=B+0.1CC=C+0.1(sinAe2A – 2B + 2C)B=DA=A+0.1