Giáo trình Turbo C nâng cao và C++ - Chương 12: Tính gần đúng đạo hàm và tích phân xác định

204 Ch−ơng 12 : Tính gần đúng đạo hàm và tích phân xác định Đ1. Đạo hàm Romberg Đạo hàm theo ph−ơng pháp Romberg là một ph−ơng pháp ngoại suy để xác định đạo hàm với một độ chính xác cao . Ta xét khai triển Taylor của hàm f(x) tại (x+h) và (x-h) : ⋅⋅⋅++′′′+′′+′+=+ )x(f !4 h )x(f !3 h )x(f 2 h )x(fh)x(f)hx(f )4( 432 (1) ⋅⋅⋅−+′′′−′′+′−=− )x(f !4 h )x(f !3 h )x(f 2 h )x(fh)x(f)hx(f )4( 432 (2) Trừ (1) cho (2) ta có : ⋅⋅⋅++′′′+′=−−+ )x(f !5 h2 )x(f !3 h2 )x(fh2)hx(f)hx(f )5( 53 (3) Nh− vậy rút ra : ⋅⋅⋅−−′′′−−−+=′ )x(f !5 h )x(f !3 h h2 )hx(f)hx(f )x(f )5( 42 (4) hay ta có thể viết lại : [ ] ⋅⋅⋅++++−−+=′ 664422 hahaha)hx(f)hx(fh21)x(f (5) trong đó các hệ số ai phụ thuộc f và x . Ta đặt : )]hx(f)hx(f[ h2 1)h( −−+ϕ = (6) Nh− vậy từ (5) và (6) ta có : ⋅⋅⋅−−−−′=ϕ= 664422 hahaha)x(f)h()1,1(D (7) ⋅⋅⋅−−−−′=⎟⎠ ⎞⎜⎝ ⎛ϕ= 64 h a 16 h a 4 h a)x(f 2 h )1,2(D 6 6 4 4 2 2 (8) và tổng quát với hi = h/2 i-1 ta có : ⋅⋅⋅−−−−′=ϕ= 6i64i42i2i hahaha)x(f)h()1,i(D (9) Ta tạo ra sai phân D(1,1) - 4D(2,1) và có : ⋅⋅⋅−−−′−=⎟⎠ ⎞⎜⎝ ⎛ϕ−ϕ 6644 ha16 15 ha 4 3 )x(f3 2 h 4)h( (10) Chia hai vế của (10) cho -3 ta nhận đ−ợc : ⋅⋅⋅+++′=−= 6644 ha16 5 ha 4 1 )x(f 4 )1,1(D)1,2(D4 )2,2(D (11) Trong khi D(1,1) và D(2,1) sai khác f′(x) phụ thuộc vào h2 thì D(2,2) sai khác f′(x) phụ thuộc vào h4 . Bây giờ ta lại chia đôi b−ớc h và nhận đ−ợc : D f x a h a h(2, ) ( ) ( / ) ( / ) ...2 1 4 2 5 16 24 4 6 6= + + +′ (12) và khử số hạng có h4 bằng cách tạo ra : D D f x a h(2, ) ( , ) ( ) ( ) ...2 16 32 15 15 64 6 6− − ′= + + + (13) Chia hai vế của (13) cho -15 ta có : D D D f x a h(3, ) (3, ) (2, ) ( ) . ...3 16 2 2 15 1 64 6 6= = − − − ′ (14) 205 Với lần tính này sai số của đạo hàm chỉ còn phụ thuộc vào h6 . Lại tiếp tục chia đôi b−ớc h và tính D(4,4) thì sai số phụ thuộc h8 . Sơ đồ tính đạo hàm theo ph−ơng pháp Romberg là : D(1,1) D(2,1) D(2,2) D(3,1) D(3,2) D(3,3) D(4,1) D(4,2) D(4,3) D(4,4) . . . . . . . . . . . . trong đó mỗi giá trị sau là giá trị ngoại suy của giá trị tr−ớc đó ở hàng trên . Với 2 ≤ j ≤ i ≤ n ta có : D j D j D jj j(i, ) (i, ) (i , ) = − − − − − − − 1 1 4 1 1 1 4 1 và giá trị khởi đầu là : D h h f x h f x hi i i i (i, ) ( ) [ ( ) ( )]1 12= = + − −ϕ với hi = h/2 i-1 . Chúng ta ngừng lại khi hiệu giữa hai lần ngoại suy đạt độ chính xác yêu cầu. Ví dụ : Tìm đạo hàm của hàm f(x) = x2 + arctan(x) tại x = 2 với b−ớc tính h = 0.5 . Trị chính xác của đạo hàm là 4.2 201843569.4)]75.1(f)25.2(f[ 25.02 1 )1,2(D 207496266.4)]5.1(f)5.2(f[ 5.02 1 )1,1(D =−ì= =−ì= 200458976.4)]875.1(f)125.2(f[ 125.02 1 )1,3(D =−ì= 200492284.4 14 )2,2(D)2,3(D4 )3,3(D 200458976.4 14 )1,2(D)1,3(D4 )2,3(D 19995935.4 14 )1,1(D)1,2(D4 )2,2(D 21 2 1 1 1 1 == == == − − − − − − Ch−ơng trình tính đạo hàm nh− d−ới đây . Dùng ch−ơng trình tính đạo hàm của hàm cho trong function với b−ớc h = 0.25 tại xo = 0 ta nhận đ−ợc giá trị đạo hàm là 1.000000001. Ch−ơng trình12-.1 //Daoham_Romberg; #include #include #include #define max 11 float h; void main() { float d[max]; int j,k,n; float x,p; float y(float),dy(float); 206 clrscr(); printf("Cho diem can tim dao ham x = "); scanf("%f",&x); printf("Tinh dao ham theo phuong phap Romberg\n"); printf("cua ham f(x) = th(x) tai x = %4.2f\n",x); n=10; h=0.2; d[0]=dy(x); for (k=2;k<=n;k++) { h=h/2; d[k]=dy(x); p=1.0; for (j=k-1;j>=1;j--) { p=4*p; d[j]=(p*d[j+1]-d[j])/(p-1); } } printf("y'= %10.5f\n",d[1]); getch(); } float y(float x) { float a=(exp(x)-exp(-x))/(exp(x)+exp(-x)); return(a); } float dy(float x) { float b=(y(x+h)-y(x-h))/(2*h); return(b); } Đ2. Khái niệm về tích phân số Mục đích của tính tích phân xác định là đánh giá định l−ợng biểu thức : J f x a b = ∫ ( )dx trong đó f(x) là hàm liên tục trong khoảng [a,b] và có thể biểu diễn bởi đ−ờng cong y= f(x). Nh− vậy tích phân xác định J là diện tích SABba , giới hạn bởi đ−ờng cong f(x) , trục hoành , các đ−ờng thẳng x = a và x = b . Nếu ta chia đoạn [a,b] thành n phần bởi các điểm xi thì J là gới hạn của tổng diện tích các hình chữ nhật f(xi).(xi+1 - xi) khi số điểm chia tiến tới ∝, nghĩa là : a a b A B y x 207 J f x x x n i i n i i = − →∞ = + ∑lim ( )( ) 0 1 Nếu các điểm chia xi cách đều , thì ( xi+1- xi ) = h . Khi đặt f(xo) = fo , f(x1) = f1 ,... ta có tổng : n i i n S h f= = ∑ 0 Khi n rất lớn , Sn tiến tới J . Tuy nhiên sai số làm tròn lại đ−ợc tích luỹ . Do vậy cần phải tìm ph−ơng pháp tính chính xác hơn . Do đó ng−ời ta ít khi dùng ph−ơng pháp hình chữ nhật nh− vừa nêu . Đ3. Ph−ơng pháp hình thang Trong ph−ơng pháp hình thang , thay vì chia diện tích SABba thành các hình chữ nhật , ta lại dùng hình thang . Ví dụ nếu chia thành 3 đoạn nh− hình vẽ thì : S3 = t1 + t2 + t3 trong đó ti là các diện tích nguyên tố . Mỗi diện tích này là một hình thang : ti = [f(xi) + f(xi-1)]/ (2h) = h(fi - fi-1) / 2 Nh− vậy : S3 = h[(fo+f1)+(f1+f2)+(f2+f3)] / 2 = h[fo+2f1+2f2+f3] / 2 Một cách tổng quát chúng ta có : )f2f2f2f(n ab S n1n1on ++⋅⋅⋅++= − − hay : }f2ff{n ab S n 1i ion n ∑++− == Một cách khác ta có thể viết : f x dx f x hf a kh f a k h a b a kh a k h k n k n ( ) ( )dx { ( ) / [ ( ) ] / } ( ) ∫ ∫∑ ∑= ≈ + + + + + + + = − = −1 1 1 0 1 2 1 2 hay : f x h f a f a h f a n h f b a b ( )dx { ( ) / ( ) [ ( ) ] ( ) / }= + + + ⋅ ⋅ ⋅ + + − +∫ 2 1 2 Ch−ơng trình tính tích phân theo ph−ơng pháp hình thang nh− sau : Ch−ơng trình 12-2 //tinh tich phan bang phuong phap hinh_thang; #include #include #include float f(float x) { float a=exp(-x)*sin(x); return(a); }; 208 void main() { int i,n; float a,b,x,y,h,s,tp; clrscr(); printf("Tinh tich phan theo phuong phap hinh thang\n"); printf("Cho can duoi a = "); scanf("%f",&a); printf("Cho can tren b = "); scanf("%f",&b); printf("Cho so buoc n = "); scanf("%d",&n); h=(b-a)/n; x=a; s=(f(a)+f(b))/2; for (i=1;i<=n;i++) { x=x+h; s=s+f(x); } tp=s*h; printf("Gia tri cua tich phan la : %10.6f\n",tp); getch(); } Dùng ch−ơng trình này tính tích phân của hàm cho trong function trong khoảng [0 , 1] với 20 điểm chia ta có J = 0.261084. Đ4. Công thức Simpson Khác với ph−ơng pháp hình thang , ta chia đoạn [a,b] thành 2n phần đều nhau bởi các điểm chia xi : a = xo < x1 < x2 < ....< x2n = b xi = a+ih ; h = (b - a)/ 2n với i =0 , . . , 2n Do yi = f(xi) nên ta có : ∫∫∫ ∫ − +++= x x fdx... x x fdx b a x x fdxdx)x(f n2 2n2 4 2 2 0 Để tính tích phân này ta thay hàm f(x) ở vế phải bằng đa thức nội suy Newton tiến bậc 2 : y t2 )1t(t yty)x(P 0 2 002 ∆−∆ ++= và với tích phân thứ nhất ta có : dx)x(Pdx)x(f x x x x 2 0 2 0 2∫∫ = Đổi biến x = x0+th thì dx = hdt , với x0 thì t =0 và với x2 thì t = 2 nên : 209 |]y) 2 t 3 t( 2 1y 2 tty[h dt)y 2 )1t(1 yty(hdx)x(P 2t 0t0 2 23 0 2 0 0 2 0 2 0 02 x x 2 0 = =∆∆ ∆−∆∫∫ −++= ++= ]yy4y[ 3 h ]y) 2 4 3 8( 2 1y2y2[h 210 0 2 00 ++= −++= ∆∆ Đối với các tích phân sau ta cũng có kết quả t−ơng tự : ]yy4y[ 3 hdx)x(f 2i21i2i2 x x 2i2 i2 ++ ++=∫ + Cộng các tích phân trên ta có : ]y)yyy(2)yyy(4y[ 3 hdx)x(f n22n2421n231o b a ++⋅⋅⋅++++⋅⋅⋅+++= −−∫ Ch−ơng trình dùng thuật toán Simpson nh− sau : Ch−ơng trình 12-3 //Phuong phap Simpson; #include #include #include float y(float x) { float a=4/(1+x*x); return(a); } void main() { int i,n; float a,b,e,x,h,x2,y2,x4,y4,tp; clrscr(); printf("Tinh tich phan theo phuong phap Simpson\n"); printf("Cho can duoi a = "); scanf("%f",&a); printf("Cho can tren b = "); scanf("%f",&b); printf("Cho so diem tinh n = "); scanf("%d",&n); h=(b-a)/n; x2=a+h; x4=a+h/2; y4=y(x4); y2=y(x2); for (i=1;i<=n-2;i++) { 210 x2+=h; x4+=h; y4+=y(x4); y2+=y(x2); } y2=2*y2; y4=4*(y4+y(x4+h)); tp=h*(y4+y2+y(a)+y(b))/6; printf("Gia tri cua tich phan la : %10.8f\n",tp); getch(); } Dùng ch−ơng trình này tính tích phân của hàm trong function trong đoạn [0,1] với 20 khoảng chia cho ta kết quả J = 3.14159265.