CHƯƠNG 6: ĐẠO HÀM VÀ TÍCH PHÂN
CHƯƠNG 6: ĐẠO HÀM VÀ TÍCH PHÂN SỐ
§1. TÍNH ĐẠO HÀM BẬC NHẤT BẰNG PHƯƠNG PHÁP 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): h2 h3 h4 f( x + h) = f( x) + hf ′( x) + f ′′( x) + f ′′′( x) + f ( 4 ) ( x) + ⋅ ⋅ ⋅ (1) 2 3! 4! h2 h3 h4 ′( x) + f ′′( x) − f ′′′( x) + f ( 4 ) ( x) − ⋅ ⋅ ⋅ f( x − h) = f( x) − hf (2) 2 3! 4!...
CHƯƠNG 6: ĐẠO HÀM VÀ TÍCH PHÂN SỐ
§1. TÍNH ĐẠO HÀM BẬC NHẤT BẰNG PHƯƠNG PHÁP 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):
h2 h3 h4
f( x + h) = f( x) + hf ′( x) + f ′′( x) + f ′′′( x) + f ( 4 ) ( x) + ⋅ ⋅ ⋅ (1)
2 3! 4!
h2 h3 h4
f( x − h) = f( x) − hf ′( x) + f ′′( x) − f ′′′( x) + f ( 4 ) ( x) − ⋅ ⋅ ⋅ (2)
2 3! 4!
Trừ (1) cho (2) ta có:
2h 3 2h 5 ( 5)
f( x + h ) − f( x − h) = 2 hf ′( x) + f ′′′( x) + f ( x) + ⋅ ⋅ ⋅ (3)
3! 5!
Như vậy rút ra:
f( x + h) − f( x − h) h 2 h4
f ′( x) = − f ′′′( x) − f ( 5 ) ( x) − ⋅ ⋅ ⋅ (4)
2h 3! 5!
hay ta có thể viết lại:
1
f ′( x) = [f( x + h) − f( x − h)] + a 2 h 2 + a 4 h 4 + a 6 h 6 + ⋅ ⋅ ⋅ (5)
2h
trong đó các hệ số ai phụ thuộc f và x.
Ta đặt:
1
ϕ( h) = [f( x + h) − f( x − h)] (6)
2h
Như vậy từ (5) và (6) ta có:
D(1,1) = ϕ( h) = f ′( x) − a 2 h 2 − a 4 h 4 − a 6 h 6 − ⋅ ⋅ ⋅ (7)
⎛ h ⎞ = f ′( x) − a h − a h − a h − ⋅ ⋅ ⋅
2 4 6
D( 2 ,1) = ϕ⎜ ⎟ (8)
⎝2⎠
2 4 6
4 16 64
và tổng quát với hi = h/2i‐1 ta có :
D(i ,1) = ϕ( h i ) = f ′( x) − a 2 h i2 − a 4 h i4 − a 6 h 6 − ⋅ ⋅ ⋅
i (9)
Ta tạo ra sai phân D(1,1) ‐ 4D(2,1) và có:
ϕ( h) − 4ϕ⎛ ⎞ = −3f ′( x) − a 4 h 4 − a 6 h 6 − ⋅ ⋅ ⋅
h 3 15
⎜ ⎟ (10)
⎝2⎠ 4 16
Chia hai vế của (10) cho ‐3 ta nhận được:
4 D( 2 ,1) − D(1,1) 1 5
D( 2 ,2) = = f ′( x) + a 4 h 4 + a 6 h 6 + ⋅ ⋅ ⋅ (11)
4 4 16
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:
311
4 6
D( 3 ,2) = f ′( x) + a 4 ⎛ ⎞ + a 6 ⎛ ⎞ + ⋅ ⋅ ⋅
1 h 5 h
⎜ ⎟ ⎜ ⎟ (12)
4 ⎝ 2 ⎠ 16 ⎝ 2 ⎠
và khử số hạng có h4 bằng cách tạo ra:
15
D( 2 ,3) − 16 D( 3,2) = −15f ′( x) + ⋅ ⋅ ⋅ + a 6 h 6 (13)
64
Chia hai vế của (13) cho ‐15 ta có:
16 D( 3,2) − D( 2 ,2) 1
D( 3 ,3) = = f ′( x) − a 6 h 6 − ⋅ ⋅ ⋅ (14)
15 64
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ó:
4 j−1 D(i , j − 1) − D(i − 1, j − 1)
D(i , j) =
4 j−1 − 1
và giá trị khởi đầu là:
1
D(i , j) = ϕ( h i ) = [f( x + h i ) − f( x − h i )]
2h i
với hi = h/2i‐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.
Ta xây dựng hàm diffromberg() để thực hiên thuật toán trên:
function df = diffromberg(f, x, h, maxiter, tol)
%Tinh dao ham bang phuong phap Romberg
D(1, 1) = (feval(f,x+h) ‐ feval(f, x‐h))/(2*h);
for i = 1:maxiter
h = h/2;
D(i + 1, 1) = (feval(f,x+h) ‐ feval(f, x‐h))/(2*h);
for j = 1:i
D(i + 1, j + 1) = (4^j*D(i + 1, j) ‐ D(i, j))/(4^j ‐ 1);
312
end
if (abs( D(i + 1, i + 1) ‐ D(i, i) ) 2 2 D(2) (x, h) − D(2) (x, 2h) −f(x + 2h) + 16f(x + h) − 30f(x) + 16f(x + h) − f(x − 2h)
c2 c2
=
2 −1
2
12h 2
h 4 ( 5)
= f ′′(x) − f (x) +L
90
Do vậy:
−f(x + 2h) + 16f(x + h) − 30f(x) + 16f(x + h) − f(x − 2h)
D(2) (x, h) =
c2
12h 2
(4)
h 4 ( 5)
= f ′′(x) − f (x) +L
90
Nếu đạo hàm cấp được tính theo (4) thì sai số chỉ còn cỡ h4. Từ (4) ta có:
c f +c f +c f +c f +c f
D(2) (x, h) = 2 2 1 1 0 0 2 −1 −1 −2 −2
c2 (5)
h
Trong đó:
f2 = f(x + 2h)
f1 = f(x + h)
f0 = f(x)
f‐1 = f(x ‐ h)
f‐2 = f(x ‐ 2h)
Viết rõ các khai triển Taylor của f2, f1, f0, f‐1, f‐2 ta có:
⎧ ⎡ (2h)2 ⎤ ⎡ h2 ⎤ ⎫
⎪ c 2 ⎢f0 + 2hf0′ + f0′′ + L⎥ + c1 ⎢ f0 + hf0′ + f0′′ + L⎥ ⎪
1 ⎪ ⎣ 2! ⎦ ⎣ 2! ⎦ ⎪
Dc2 (x, h) = 2 ⎨
(2)
⎬
h ⎪ ⎡ h2 ⎤ ⎡ ( −2h)2 ⎤⎪
+c f + c −1 ⎢f0 − hf0′ + f0′′ − L⎥ + c −2 ⎢f0 − 2hf0′ + f0′′ − L⎥
⎪ 00 ⎣ 2! ⎦ ⎣ 2! ⎦⎪
⎩ ⎭
⎧(c 2 + c1 + c 0 + c −1 + c −2)f0 + h(2c 2 + c1 + c −1− 2c −2 )f0′ ⎫
1 ⎪ ⎪
D(2) (x, h) = 2 ⎨ 2 ⎛ 2 2 1 1 22 ⎞ ⎬ (6)
h ⎪+ h ⎜ c 2 + c1 + c −1 + c −2 ⎟ f0′′ + L
c2
⎪
⎩ ⎝ 2 2 2 2 ⎠ ⎭
Ta phải giải hệ phương trình sau để tìm các hệ số ci.
⎡ 1 1 1 1 1 ⎤ ⎡c 2 ⎤ ⎡0 ⎤
⎢ 2 1 0 −1 −2 ⎥ ⎢ c 1 ⎥ ⎢ 0 ⎥
⎢ 2 ⎥⎢ ⎥ ⎢ ⎥
⎢ 2 2! 1 2! 0 1 2! 2 2 2! ⎥ ⎢c 0 ⎥ = ⎢1 ⎥ (7)
⎢ 2 ⎥⎢ ⎥ ⎢ ⎥
⎢ 2 3! 1 3! 0 −1 3! − 2 3!⎥ ⎢c −1 ⎥ ⎢0 ⎥
3
⎢ 2 4 4! 1 4! 0 1 4! − 2 4 4!⎥ ⎢c −2 ⎥ ⎢0 ⎥
⎣ ⎦⎣ ⎦ ⎣ ⎦
314
Kết quả ta có c2 = ‐1/12, c1 = 4/3, c0 = ‐5/2, c‐1 = 4/3 c‐2 = ‐1/12. Do vậy:
−f + 16f1 − 30f0 + 16f−1 − f−2
D(2) (x, h) = 2
c2
12h 2
Tương tự ta có đạo hàm bậc 4 của hàm:
f − 4f1 + 6f0 − 4f−1 + f−2
D(4) (x, h) = 2
c2
12h 4
Ta xây dựng hàm diffn() để tính đạo hàm tới bậc 5:
function df = diffn(f, n, x)
% Tinh dao ham cap n cua f tai x
if n>5
error(ʹHam chi tinh duoc dao ham den bac 5ʹ);
return;
end;
N = 5;
xo = x;
T(1) = feval(f,xo);
h = 0.005;
tmp = 1;
for i = 1:N
tmp = tmp*h;
c = difapx(i,[‐i i]); %he so cua dao ham
dix = c*feval(f,xo + [‐i:i]*h)ʹ;
T(i+1) = dix/tmp; %dao ham
end
df = T(n+1);
h = 0.005;
tmp = 1;
for i = 1:N
tmp = tmp*h;
c = difapx(i,[‐i i]); %he so cua dao ham
dix = c*feval(f,xo + [‐i:i]*h)ʹ; %/h^i; %dao ham
T(i+1) = dix/tmp; %he so cua chuoi Taylor
end
df = T(n+1);
315
Để tính đạo hàm của hàm ta dùng chương trình ctdiffn.m
clear all, clc
f = inline(ʹx.^2 + atan(x)ʹ,ʹxʹ);
df = diffn(f, 5, 0)
§3. TÍNH ĐẠO HÀM BẰNG PHƯƠNG PHÁP NỘI SUY
Giả sử ta có hàm cho dưới dạng bảng:
x x0 x1 x0 ... xn
y y0 y1 y0 ... yn
Để tìm đạo hàm của hàm tại một điểm nào đó ta sẽ nội suy hàm rồi sau đó
tính đạo hàm của hàm tại điểm đã cho. Ta xây dựng hàm diffinterp() để thực
hiện công việc trên.
function df = diffinterp(x, y, n, x0)
%Tinh dao ham cap 1 hai 2 bang phuogphap noi suy
px = lagrange(x, y); % Tim da thuc noi suy Lagrange (x, y)
[p, dp, ddp] = peval(px, x0);
fprintf(ʹTri so cua ham la: %f\nʹ,p)
if n ==1
df = dp;
else
df = ddp;
end
fprintf(ʹDao ham cap %d la: %f\nʹ,n, df);
Để tính đạo hàm ta dùng chương trình ctdiffinterp.m:
clear, clc
x0 = pi/4;
x = [2:6]*pi/16;
y = sin(x);
x = [1.5 1.9 2.1 2.6 3.2];
y = [1.0628 1.3961 1.5432 1.8423 2.0397];
316
n = 2;
df = diffinterp(x, y, n, x0);
§4. TÍCH PHÂN XÁC ĐỊNH
Mục đích của tính tích phân xác định, còn gọi là cầu phương, là đánh
giá định lượng biểu thức:
b
y
J = ∫ f( x)dx B
a
trong đó f(x) là hàm liên tục trong khoảng A
[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 b
a x
hoành, các đường thẳng x = a và x = b. Tích
phân này thường được tính gần đúng bằng
công thức:
n
J = ∑ A i f(xi )
i =1
trong đó Ai là trọng số, phụ thuộc phương pháp tính tích phân.
Tất cả các phương pháp tính tích phân được suy ra từ phương pháp nội
suy hàm dưới dấu tích phân. Do vậy kết quả sẽ chính xác nếu hàm có thể xấp
xỉ bằng đa thức. Các phương pháp tính tích phân xác định bằng phương pháp
số được chia thành 2 nhóm: các phương pháp Newton ‐ Cotes và các phương
pháp Gauss. Khi dùng các phương pháp Newton ‐ Cotes khoảng lấy tích phân
được chia đều như trong phương pháp hình thang hay phương pháp
Simpson. Khi dùng các phương pháp Gauss, cácc diểm chia được chọn để đạt
độ chính xác cao nhất. Do phương pháp này cần ít lần tính giá trị hàm dươci
dấu tích phân nên thích hợp khi hàm f(x) khó tính.
§5. CÁC CÔNG THỨC NEWTON ‐ COTES
1. Khái niệm chung: Ta khảo sát tích phân
b
J= ∫ f(x)dx (1)
a
Ta chia miền lấy tích phân [a, b] thành (n ‐ 1) đoạn bằng nhau có
chiều dài mỗi đoạn h = (b ‐ a)/(n ‐ 1) như hình vẽ và kí hiệu các điểm chia là
317
x1, x2,.., xn. Sau đó ta xấp xỉ hàm f(x) bằng đa
thức bậc (n ‐ 1) đi qua các nút. Đa thức nội
suy Lagrange của f(x) có dạng:
n
Pn −1 (x) = ∑ f(xi )Li (x)
i =1
Như vậy, xấp xỉ tích phân (1) là:
n x1 x2 x3 xn
∑ f(x )∫ L (x)dx = ∑ A f(x )
b b b n
J= ∫ f(x)dx = ∫ Pn‐1 (x)dx = i i i i (2)
a a i =1 a i =1
Trong đó:
b
A i = ∫ Li (x)dx i = 1,2,...,n (3)
a
Công thức (2) là công thức Newton ‐ Cotes. Với n = 2 ta có công thức hình
thang và với n = 3 ta có công thức Simpson.
2. Phương pháp hình thang: Khi n = 2 ta có:
x − x2 x−b
L1 (x) = =−
x1 − x 2 h
x − x1 x − a
L 2 (x) = = h
x 2 − x1 h
b x1 = a x2 = b
∫
1 1 h
A1 = − (x − b)dx = (b − a)2 =
h 2h 2
a
b
∫
1 1 h
A2 = (x − a)dx = (b − a)2 =
h 2h 2
a
Vậy:
J = ⎡f(a) + f(b) ⎤
h
⎢
⎣ ⎥2
⎦
Trong thực tế, phương pháp hình thang được áp dụng trên từng đoạn. Trên
mỗi đoạn [xi, xi+1] ta có:
J i = ⎡f(x i ) + f(xi+1 )⎤
h
⎢
⎣ ⎥2
⎦
n
và: J = ∑ J i = ⎡f(x1 ) + 2f(x 2 ) + 2f(x 3 ) + L + 2f(x n −1 ) + 2f(x n ) ⎤
h
(7)
i =1
⎢
⎣ ⎥2
⎦
Ta gọi H = b ‐ a. Nếu tích phân trên được tính chỉ bởi k hình thang thì:
J1 = ⎡f(a) + f(b) ⎤
H
k = 1: (8)
⎢
⎣ ⎥2
⎦
318
⎡ ⎤
⎛ a + H ⎞ + f(b) ⎥ H = 1 J + f ⎛ a + H ⎞ H
J 2 = ⎢f(a) + 2f ⎜
k = 2: ⎟ ⎜ ⎟
⎝ 2⎠ ⎝ 2⎠2
1
⎢ ⎥4 2
⎣ ⎦
⎡ ⎤H
J 2 = ⎢f(a) + 2f ⎛ a + ⎞ + 2f ⎛ a + ⎞ + 2f ⎛ a +
H H 3H ⎞
⎜ ⎟ ⎜ ⎟ ⎜ ⎟ + f(b) ⎥
⎢ ⎝ 4⎠ ⎝ 2⎠ ⎝ 4 ⎠ ⎥8
⎣ ⎦
k = 3:
1 ⎡ 3H ⎞ ⎤ H
= J 2 + ⎢f ⎛ a + ⎞ + f ⎛ a +
H
⎜ ⎟ ⎜ ⎟⎥
2 ⎢ ⎝ 4⎠ ⎝ 4 ⎠⎥ 4
⎣ ⎦
Tổng quát, với k > 1 ta có:
2 k −1
Jk =
1
2 2
H
J k−1 + k−1
∑
i =1
f ⎡a +
⎢
⎣
(2i − 1)H ⎤
2 k −1 ⎥
⎦
k = 2,3,... (9)
Công thức (8) là công thức hình thang lặp. Ta thấy rằng tổng chỉ chứa các nút
mới xuất hiện khi số hình thang tăng gấp đôi. Tính dãy J1, J2,... bằng (8) và (9)
cần cùng một số lần tính như khi dùng (7). Nhưng khi dùng (8) và (9) ta kiểm
tra được tính hội tụ và có thể dừng lặp khi đạt độ chính xác cho trước.
Ta xây dựng hàm trapezoid() để thực hiện thuật toán trên.
function J = trapezoid(f, a, b, maxiter, tol)
% Quy tac hinh thang lap.
% Cu phap: J = trapezoid(f, a, b, k)
fa = feval(f, a);
fb = feval(f, b);
J1 = (fa + fb)*(b ‐ a)/2;
for k = 2:maxiter
n = 2^(k ‐2 ); % so diem moi
h = (b ‐ a)/n ; % khoang chia moi
x = a + h/2.0; % toa do diem moi thu nhat
sum = 0.0;
for i = 1:n
fx = feval(f, x);
sum = sum + fx;
x = x + h;
end
319
J = (J1 + h*sum)/2;
if abs(J1 ‐ J) t(t − 1) 2 ⎤
b 2
⎡
∫ P (x)dx = h ∫ ⎢ y
a
2
⎣
0
0 + t∆y 0 +
2!
∆ y 0 ⎥ dt
⎦
t =2
⎡ t2 1 ⎛ t3 t2 ⎞ 2 ⎤
= h ⎢ y 0 t + ∆y 0 + ⎜ − ⎟ ∆ y 0 ⎥
⎣ 2 2⎝ 3 2 ⎠ ⎦ t =0 (12)
⎡ ⎤
= h ⎢ 2y 0 + 2∆y 0 + ⎛ − ⎞ ∆ 2 y 0 ⎥
1 8 4
⎜ ⎟
⎣ 2⎝ 3 2⎠ ⎦
a+b⎞
( y0 + 4y1 + y 2 ) = ⎡f(a) + 4f ⎛ ⎤
h h
= ⎜ ⎟ + f(b) ⎥
3 3⎢ ⎣ ⎝ 2 ⎠ ⎦
Thực tế ta chia đoạn [a, b] thành 2n phần và tính tích phân trên mỗi đoạn con.
Cộng các tích phân trên các đoạn con ta có:
b
h⎡ ⎤
∫ f(x)dx = 3 ⎢ y0 + 4 ( y1 + y3 + ⋅ ⋅ ⋅ + y2n−1 ) + 2 ( y2 + y4 + ⋅ ⋅ ⋅ + y2n−2 ) + y 2n ⎥ (13)
a
⎣ ⎦
Công thức (13) đòi hỏi n là số chẵn.
Ta xây dựng hàm simpson() để thực hiện thuật toán trên
function s = simpson(f, a, b, n)
%n so khoang chia
%neu f chua trong mot file dung ki hieu @ de goi
% s = simpson(@f, a, b, n).
%neu f la ham inline
% s = simpson(f, a, b, n).
if mod(n, 2) ~= 0
n = n + 1
end
h = (b ‐ a)/(2*n);
s1 = 0;
s2 = 0;
for k = 1:n
x = a + h*(2*k‐1);
s1 = s1+ f(x);
end
for k = 1:(n‐1)
x = a + h*2*k;
s2 = s2 + f(x);
321
end
s = h*(f(a) + f(b) + 4*s1 + 2*s2)/3;
clc
Để tính tích phân ta dùng chương trình ctsimpson.m:
clear all, clc
f = inline(ʹexp(x).*sin(x)ʹ,ʹxʹ);
a = 0;
b = 1;
n = 6;
s = simpson(f, a, b, n)
3. Phương pháp cầu phương thích nghi: Trong
tích phân bằng phương pháp Simpson, các
đoạn được chia đều và làm cho sai số không
giống nhau trên các đoạn: sai số lớn trên các
đoạn hàm biến đổi nhiều và sai số nhỏ trên các
đoạn hàm tương đối bằng phẳng. Ngược lại
phương pháp cầu phương thích nghi chia các đoạn không đều: ngắn trên các
đoạn hàm thay đổi nhiều và dài trên các đoạn thay đổi ít và sẽ có sai số nhỏ
khi số đoạn chia nhỏ.
Thuật toán cầu phương thích nghi bắt đầu bằng việc tính tích phân int
đối với toàn bộ đoạn [a, b] và tổng tích phân int12 = int1 + int2 trên 2 đoạn
bằng nhau. Dựa trên int và int12 ta tính sai số. Nếu chưa đạt độ chính xác, ta
chia đôi mỗi đoạn và lặp lại quá trình tính. Ta dùng hàm adaptivesimpson()
để thực hiện thuật toán này:
function int = adaptivesimpson(f, a, b, tol)
mid = (b + a)/2.0;
int = simpsonapprox (f, a, b);
int12 = simpsonapprox (f, a, mid) + simpsonapprox (f, mid, b);
if( abs(int ‐ int12) int = leftint + rightint;
end
function int = simpsonapprox (f, a, b)
h = (b ‐ a)/2.0;
int = h*( feval(f, a) + 4.0*feval(f, (a + h)) + feval(f, b) )/3.0;
Để tính tích phân ta dùng chương trình ctadaptive.m:
clc, clear all
f = inline(ʹsqrt(x).*cos(x)ʹ);
a = 0;
b = 1;
tol = 1e‐5;
J = adaptivesimpson(f, a, b, tol)
§6. TÍCH PHÂN ROMBERG
Tích phân Romberg kết hợp quy tắc tích phân hình thang với phương
pháp ngoại suy Richardson. Trước hết ta đưa vào khái niệm:
Ri,1 = Ji
b
Trong đó Ji là giá trị xấp xỉ của ∫ f(x)dx có được bằng cách tính theo quy tắc
a
lặp hình thang lần thứ i.
Tích phân Romberg bắt đầu từ R1,1 = J1 (một hình thang) và R2,1 = J2 (hai
hìn thang). Sau đó tính R2,2 bằng cách ngoại suy:
2 2 R 2 ,1 − R 1,1 4 1
R 2 ,2 = = R 2 ,1 − R 1,1 (1)
2 −1
2
3 3
Để tiện dùng ta lưu các kết quả vào mảng dạng:
⎡ R 1,1 ⎤
⎢R
⎣ 2 ,1 R 2 ,2 ⎥
⎦
Bước tiếp theo là tính R3,1 = J3 (bốn hình thang) và lặp lại ngoại suy Richadson
ta có:
2 2 R 3,1 − R 2 ,1 4 1
R 3,2 = = R 3,1 − R 2 ,1 (2)
2 −1
2
3 3
323
2 4 R 3,2 − R 2 ,2 16 1
và: R 3,3 = = R 3,2 − R 2 ,2 (3)
24 − 1 15 15
Các phần tử của R bây giờ gồm:
⎡ R 1,1 ⎤
⎢ ⎥
⎢R 2 ,1 R 2 ,2 ⎥
⎢ R 3,1 R 3,2 R 3,3 ⎥
⎣ ⎦
Công thức tổng quát dùng trong sơ đồ này là:
4 j−1 R i ,j−1 − R i−1,j−1
R i ,j = i > 1, j = 2, 3,... (4)
4 j−1 − 1
Ta xây dựng hàm romberg() để thực hiện thuật toán trên:
function J = romberg(f, a, b, maxiter, tol)
m = 1;
h = b‐a;
err = 1;
j = 0;
R = zeros(4, 4);
R(1,1) = h*(f(a) + f(b))/2;
while((err > tol) & (j Để tính tích phân ta dùng chương trình ctromberg.m:
clear all, clc
f = inline(ʹexp(x).*sin(x)ʹ,ʹxʹ);
a = 0;
b = 1;
maxiter = 20;
tol = 1e‐6;
J = romberg(f, a, b, maxiter, tol)
§7. TÍCH PHÂN BOOL
Ta khảo sát hàm y = f(x) trên đoạn [x0, x4], với:
x1 = x0 + h, x2 = x0 + 2h, x3 = x0 + 3h, x4 = x0 + 4h
Theo Bool, tích phân:
x4
2h m
J = ∫ f(x)dx = ∑ 7f(x0 ) + 32f(x1 ) + 12f(x2 ) + 32f(x3 ) + 7f(x4 )
45 k=1
x0
Xét tích phân:
b
J = ∫ f(x)dx
a
b−a
Ta chia đoạn [a, b] thành 4m đoạn con đều nhau có độ rộng h = bởi các
4m
điểm chia xk = x0 + hk = a + hk, k = 0, 1,..., 4m. Công thức Bool cho 4m đoạn
con là:
b
2h m
J = ∫ f(x)dx = ∑ 7f(x0 ) + 32f(x1 ) + 12f(x2 ) + 32f(x3 ) + 7f(x4 )
45 k =1
a
Ta xây dựng hàm intbool() để thực hiện thuật toán này
function tp = intbool(f, a, b, m)
%Tinh tich phan bang phuong phap Bool
a = 0;
b = 2;
m = 2;
h = (b ‐ a)/(4*m);
for k = 1:4*m
x(k) = a + k*h;
325
end
tp = 0;
j = 1;
for k = 1:m
tp = tp + (7*feval(f, a) + 32*feval(f, x(j)) +...
12*feval(f, x(j+1)) + 32*feval(f, x(j+2)) + 7*feval(f, x(j+3)));
a = x(4*k);
j = 4*k + 1;
end
tp = tp*h*2/45;
Để tính tích phân của một hàm ta dùng chương trình ctintbool.m:
clear all, clc
format long
f = inline(ʹx.*sin(x)ʹ);
a = 0;
b = 2;
m = 2;
J = intbool(f, a, b, m)
§8. CÔNG THỨC TÍCH PHÂN FILON
Giả sử cần tính tích phân:
b
J = ∫ f(x)cos(ωx)dx
a
Lúc đó ta có thể dùng công thức tích phân Filon:
xn
∫ f(x)cos(tx)dx
x0
{
= h α(th) [ f2n sin(tx 2n )‐f2n sin(tx 2n )] + β(th)C 2n + γ(th)C 2n −1 +
2 4
45
th S′2n }
Trong đó:
a = x0, b = xn, t = ω
n
C 2n = ∑ f2i cos(tx 2i ) − 0.5 ⎡f2n cos(tx 2n ) + f0cos(tx0 )⎤
⎣ ⎦
i =0
326
n
C 2n −1 = ∑ f2i−1cos(tx 2i‐1 )
i =0
n
S′2n −1 = ∑ f2i−1 sin(tx 2i−1 )
′′′
i =1
1 sin 2θ sin 2θ
α(θ) = + − 3
θ 2θ 2 θ
⎡ 1 + cos θ sin 2θ ⎤
2
β(θ) = 2 ⎢ −
⎣ θ2 θ3 ⎥⎦
sin θ cos θ
γ(θ) = 4 ⎛ 3 − 2 ⎞
⎜ ⎟
⎝ θ θ ⎠
Ta xây dựng hàm filon() để thực hiện các công thức trên:
function int = filon(f, a, b, t, m, key)
% ham filon tinh gan dung tich phan
b
% ∫ f (x)cos(tx)dx neu key = 1,
a
% hay
b
% ∫ f (x)sin(tx)dx neu key = 2,
a
% dung m diem theo quy tac Filon (m le).
if (any(size(a) ~= [1 1]))
error (ʹThong so a phai la so.ʹ) ;
end
if (any(size(b) ~= [1 1]))
error (Thong so b nhap vao phai la so.ʹ) ;
end
if (any(size(t) ~= [1 1]))
error (ʹThong so t phai la so.ʹ) ;
end
if (any(size(m) ~= [1 1]))
error (ʹThong so m phai la so.ʹ) ;
end
if (any([(fix(m) ~= m) (rem(m, 2) == 0)]))
error (ʹThong so m phai la so le.ʹ) ;
end
327
if (m = 0.1)
s = sin(th) ;
c = cos(th) ;
alfa = (1.0 + s*(c ‐ 2.0*s/th)/th)/th ;
beta = 2.0*(1.0 + c*c ‐ 2.0*s*c/th)/thh ;
gamma = 4.0*(s/th ‐ c)/thh ;
else
alfa = th*thh*(2.0/45.0 + thh*(‐2.0/315.0 + 2.0*thh/4725.0)) ;
beta = 2.0/3.0 + thh*(2.0/15.0 + thh*(4.0/105.0 + 2.0*thh/567.0)) ;
gamma = 4.0/3.0 + thh*(‐2.0/15.0 + thh*(1.0/210.0 ‐ thh/11340.0)) ;
end
args = [a b];
fbounds = feval(f, args) ;
s1 = sin(a*t) ;
s2 = sin(b*t) ;
c1 = cos(a*t) ;
c2 = cos(b*t) ;
if (key == 1)
sum = s2*fbounds(2) ‐ s1*fbounds(1) ;
sum0 = 0.5*(c1*fbounds(1) + c2*fbounds(2)) ;
if (n > 2)
args = (a + (2:2:n‐2)*h)ʹ ;
sum0 = sum0 + cos(t*args)ʹ*feval(f, args) ;
end
args = (a + (1:2:n‐1)*h)ʹ ;
sum1 = cos(t*args)ʹ*feval(f, args) ;
328
else
sum = c1*fbounds(1) ‐ c2*fbounds(2) ;
%sum = ‐(c1*fbounds(1) ‐ c2*fbounds(2)) ;
sum0 = 0.5*(s1*fbounds(1) + s2*fbounds(2)) ;
%if (n == 2)
if (n > 2)
args = (a + (2:2:n‐2)*h)ʹ ;
sum0 = sum0 + sin(t*args)ʹ*feval(f, args) ;
end
args = (a + (1:2:n‐1)*h)ʹ ;
sum1 = sin(t*args)ʹ*feval(f, args) ;
end
int = h*(alfa*sum + beta*sum0 + gamma*sum1) ;
Khi tính tích phân ta dùng chương trình ctintfilon.m:
clear all, clc
a = 0;
b = 2;
key = 2;
t = 3;
m = 51;
f = inline(ʹ(x.^3+1).*sin(x)ʹ);
J = filon(f, a, b, t, key)
§9. QUY TẮC HARDY
b
Để tính tích phân J = ∫ f(x)dx ta có thể dùng công thức Hardy:
a
x7
∫ f(x)dx = 0.01h ( 28f
x1
1 + 162f2 + 220f4 + 162f6 + 28f7 )
Để tăng độ chính xác ta dùng phương pháp chia đoạn [a, b] thành m đoạn và
trên mỗi đoạn ta dùng công thức Hardy. Ta xây dựng hàm inthardy() để thực
hiện công thức trên:
function tp = inthardy(f, a, b, m)
329
%Tinh tich phan bang phuong phap Hardy
h = (b ‐ a)/(6*m);
for k = 1:6*m
x(k) = a + k*h;
end
tp = 0;
j = 1;
for k = 1:m
tp = tp + (28*feval(f, a) + 162*feval(f, x(j)) +...
220*feval(f, x(j+2)) + 162*feval(f, x(j+4)) + 28*feval(f, x(j+5)));
a = x(6*k);
j = 6*k + 1;
end
tp = tp*h*0.01;
Để tính tích phân ta dùng chương trình ctinthardy.m:
clear all, clc
format long
f = inline(ʹexp(x).*sin(x)ʹ,ʹxʹ);
a = 0;
b = 2;
m = 20;
J = inthardy(f, a, b, m)
§10. QUY TẮC DURANT
b
Để tính tích phân J = ∫ f(x)dx ta có thể dùng công thức Durant:
a
xn
⎛ 2 + 11 f + f + L + f + 11 f + 2 f ⎞
∫ f(x)dx = h ⎜ 5 f
x1 ⎝
1
10
2 3 n−2
10
n −1 n⎟
5 ⎠
Ta xây dựng hàm intdurant() để thực hiện công thức trên:
function tp = intdurant(f, a, b, n)
%Tinh tich phan bang phuong phap Durant
h = (b ‐ a)/(n);
330
