Giáo trình Matlab cơ bản - Chương 6: Đạo hàm và tích phân số

cobi. Các công thức tính đa thức Jacobi 
là: 
  ( , )0P (x) 1
α β =  
  [ ]( , )1P (x) 0.5 2( 1) ( 2)(x 1)α β = α + + α + β + −  
( , ) ( , )
( , ) 2n 3n n 1 4n n 2
n
1n
(b b x)P (x) b P (x)P (x)
b
α β α β
α β − −+ −=  
Với:  b1n = 2i(i + α + β)(2i + α + β ‐2) 
  b2n = (2i + α + β ‐1)(α2 ‐ β2) 
  b3n = (2i + α + β ‐ 2)(2i + α + β ‐ 1)(2i + α + β) 
  b4n = 2(i + α ‐ 1)(i + β ‐ 1)(2i + α + β) 
Các trọng số wi được xác định bằng: 
  [ ]
2i 1
i 22
i n i
(i 1) (i 1) 2 n!w
(i 1) (1 x ) V (x )
+α+β+Γ + α + Γ + β += Γ + α + β + ′−  
với: 
n
( , )
n n n
2 n!V P (x)
( 1)
α β= −  
Ta xây dựng hàm gaussjacobi() để tìm xi và wi: 
 350
function [x, w] = gaussjacobi(n, alfa, beta) 
%tinh cac trong so va hoanh do trong tich phan Gauss‐Jacobi 
p = [0.5*(alfa + beta + 2)   0.5*(alfa ‐ beta)]; 
a  = 1; 
b = p; 
for i = 2:n+1 
    b1 = 2*i*(i + alfa + beta)*(2*i + alfa + beta ‐2); 
    b2 = (2*i + alfa + beta ‐1)*(alfa^2 ‐ beta^2)/b1; 
    b3 = ((2*i + alfa + beta ‐2)*(2*i + alfa + beta ‐1 )*(2*i + alfa + beta))/b1; 
    b4= (2*(i + alfa ‐1)*(i + beta ‐ 1)*(2*i + alfa + beta))/b1; 
    s = [b3 b2]; 
    if i == n+1 
        pn1 = conv(s, b) ‐ [0 0 b4*a]; 
        break; 
    else 
        p = conv(s, b) ‐ [0 0 b4*a]; 
    end 
    a  = b; 
    b = p; 
end 
x = roots(p); 
w = zeros(n, 1); 
dv = polyder(p); 
if mod(n, 2) == 1 
    sign = ‐1; 
else 
    sign = 1; 
end 
dv = dv*(2^n)*factorial(n)/sign; 
pn1 = ‐pn1*(2^(n+1))*factorial(n+1)/sign;     
for i = 1:n 
    num = (2*n + alfa + beta +...  
    2)*gamma(n+alfa+1)*gamma(n+beta+1)*(2^(2*n+alfa+beta+1))*factorial(n); 
    den = (n + alfa + beta + 1)*gamma(n+alfa+beta+1)*polyval(dv,... 
    x(i))*polyval(pn1, x(i)); 
    w(i) = num/den; 
 351
end 
Tiếp theo ta xây dựng hàm intgaussjacobi() để tính tích phân: 
function J = intgaussjacobi(f, n, alf, bta) 
[t, w] = gaussjacobi(n, alf, bta); 
fx = feval(f, t); 
J = wʹ*fx; 
Để tính tích phân ta dùng  chương trình ctgaussjacobi.m: 
clear al, clc 
f = inline(ʹexp(x).*sin(x)ʹ,ʹxʹ); 
n = 6;%n <= 40 
alfa = 1; 
beta = 0; 
J = intgaussjacobi(f, n, alfa, beta) 
§17. TÍCH  PHÂN RADAU 
  Cầu phương Radau dùng để tính tích phân: 
1
1
J f(x)dx
−
= ∫                    (1) 
Theo công thức cầu phương Radau ta có: 
1 n
1 i i
i 21
J f(x)dx w f( 1) w f(x )
=−
= = − + ∑∫             (2) 
Điểm (‐1) là một một nút trong số các nút của cầu phương. Các hoành độ còn 
lại là nghiệm của đa thức: 
  n 1 nP (x) P (x)
1 x
− −
+                   (3) 
Trong đó P(x) là đa thức Legendre. Các trọng số tương ứng được tính theo: 
  [ ]ii 22 n 1 i
1 xw
n P (x )−
−=                 (4) 
và điểm cuối có: 
  1 2
2w
n
=  
Ta xây dựng hàm radau() để tính các hoành độ xi và trọng số wi: 
 352
function [x , w] = radau(n) 
%Tinh cac hoanh do va trong so trong cau phuong Radau 
tol = 1e‐8; 
%danh gia ban dau cac hoanh do la ca nut Chebyshev‐Gauss‐Radau 
x(1:n,1) = ‐ cos(2.0*pi*(0:n‐1)/(2*n‐1))ʹ; 
p = zeros ( n, n+1 ); 
xold(1:n,1) = 2.0; 
while (tol < max(abs(x(1:n,1)‐xold(1:n,1)))) 
    xold = x;    
    p(1,1:n+1) = (‐1.0).^(0:n); 
    p(2:n,1) = 1.0; 
    p(2:n,2) = x(2:n,1);        
    for j = 2:n 
      p(2:n,j+1) = ((2*j ‐ 1)*x(2:n,1).*p(2:n,j)+ ( ‐j+1)*p(2:n,j‐1))/j; 
    end      
    x(2:n,1) = xold(2:n,1) ‐ ((1.0 ‐ xold(2:n,1))/n)... 
      .*(p(2:n,n)+p(2:n,n+1))./(p(2:n,n) ‐ p(2:n,n+1)); 
end 
w = zeros(n,1); 
w(1) = 2/n^2; 
Tiếp  theo  ta xây dựng hàm  intradau(), để  tính  tích phân. Trong hàm  ta đổi 
cận lấy tích phân trong khoảng [‐1, 1] thành tích phân trong khoảng: 
function J = intradau(f, n, a, b) 
[t, w] = radau(n); 
x = ((b ‐ a)*t + a + b)/2; 
fx = feval(f, x); 
J = wʹ*fx*(b ‐ a)/2; 
Để tính tích phân ta dùng chương trình ctradau.m: 
clear al, clc 
f = inline(ʹx.*sin(x)ʹ); 
n = 6;%n <= 40 
a = 1;  
 353
b = 3; 
J = intradau(f, n, a, b) 
§18. TÍCH  PHÂN CHEBYSHEV ‐ RADAU 
  Cầu phương Chebyshev – Radau dùng tính tích phân: 
1
1
J f(x)dx
−
= ∫                    (1) 
Theo công thức cầu phương Chebyshev ‐ Radau ta có: 
  [ ]1 n i i i
i 21
J f(x)dx w f(x ) f( x )
=−
= = − −∑∫             (2) 
Các hoành độ xi và trọng số tương ứng wi cho trong bảng: 
xi  0.3549416  0.6433097  0.7783202  0.9481574 
wi  0.1223363  0.1223363  0.1223363  0.1223363 
Ta xây dựng hàm chebradau() để chứa các giá trị x và w: 
function [x, w] = chebradau 
x(1) = 0.3549416; 
x(2) = 0.6433097; 
x(3) = 0.7783202; 
x(4) = 0.9481574; 
w(1) = 0.1223363; 
w(2) = 0.1223363; 
w(3) = 0.1223363; 
w(4) = 0.1223363; 
và hàm intchebradau() để tính tích phân: 
function J = intchebradau(f, a, b) 
[t, w] = chebradau; 
fx1 = feval(f, t); 
fx2 = feval(f, ‐t); 
J = (w*fx1ʹ ‐ w*fx2ʹ); 
 354
Để tính tích phân của một hàm cụ thể ta dùng chương trình ctchebradau.m: 
clear al, clc 
f = inline(ʹexp(x).*sin(x)ʹ,ʹxʹ); 
J = intchebradau(f) 
§19. TÍCH  PHÂN GAUSS – RADAU 
  Công thức cầu phương Gauss ‐ Radau có dạng: 
1 n 1
1 2 i i
i 21
f(x)dx w f( 1) w f(1) w f(x )
−
=−
= − + + ∑∫  
Ngoài hai điểm nút x = 1, các điểm nút khác  là nghiệm của đa  thức Pn(x) + 
Pn+1(x) , với P(x) là đa thức Legendre. Các trọng số được xác định bằng: 
  = +1 2
2w
(n 1)
  [ ]+
−= +
i
i 2
n 1 i
1 xw
(n 1) P (x )
Ta  xây dựng  hàm  gaussradau()  để  tính  các  hoành  độ  và  trọng  số  của  cầu 
phương: 
function [x, w] = gaussradau(N) 
% tinh cac nut va trong so cua cau phuong Gauss ‐ Radau 
N1 = N + 1; 
% dung cac nut Chebyshev‐Gauss‐Radau lam xap xi dau tien 
x = ‐cos(2*pi*(0:N)/(2*N+1))ʹ; 
P = zeros(N1, N1+1); 
xold = 2; 
free = 2:N1; 
while max(abs(x‐xold)) > eps 
    xold = x;     
    P(1, :) = (‐1).^(0:N1);    
    P(free, 1) =1;     
    P(free, 2) = x(free); 
    for k = 2:N1 
        P(free, k+1) = ( (2*k ‐ 1)*x(free).*P(free, k) ‐ (k ‐ 1)*P(free, k‐1) )/k; 
    end     
 355
x(free)  =  xold(free)  ‐  ((1‐xold(free))/N1).*(P(free,  N1)  +  P(free, 
N1+1))..../(P(free, N1) ‐ P(free, N1+1)); 
end P = P(1:N1, 1:N1); 
w = zeros(N1, 1); 
w(1) =2/N1^2; 
w(free) = (1‐x(free))./(N1*P(free, N1)).^2; 
Ta dùng hàm intgaussradau() để tính tích phân của hàm f(x) trên đoạn [a, b]: 
function J = intgaussradau(f, n, a, b) 
[t, w] = gaussradau(n); 
x = ((b ‐ a)*t + a + b)/2; 
fx = feval(f, x); 
J = wʹ*fx*(b ‐ a)/2; 
Để tính tích phân của hàm ta dùng chương trình ctgaussradau.m: 
clear all, clc 
f = inline(ʹexp(x).*sin(x)ʹ,ʹxʹ); 
n = 6; 
a = 1;  
b = 3; 
J = intgaussradau(f, n, a, b) 
§20. ĐA THỨC NỘI SUY VÀ TÍCH  PHÂN SỐ 
  Khi hàm  được  cho dưới dạng bảng  số,  để  tính  tính phân  của hàm  ta 
thực hiện các bước sau: 
  ‐ Tìm đa thức nội suy, có thể là đa thức Lagrange, đa thức Newton... 
  ‐ Tìm tích phân của đa thức nội suy 
Sau đây chúng ta sẽ xây dựng chương trình ctinterp.m để tính tích phân dùng 
đa thức nội suy Lagrange: 
clear all, clc 
x = [1.0000    1.2000    1.4000    1.6000    1.8000    2.0000]; 
y = [1.9221    1.9756    1.6517    0.8501   ‐0.4984   ‐2.4199]; 
l = lagrange(x, y); 
 356
n = length(l); 
p = conv(l, [1 0]); 
n = length(p); 
for i = 1:n‐1 
    q(i) = p(i)/(n‐i); 
end 
q(n) = p(n); 
tp = polyval(q, x(length(x))) ‐ polyval(q, x(1)); 
§21. TÍCH PHÂN KÉP 
1. Khái  niệm  chung: Ta  khảo  sát  tích  phân  của  hàm  z  =  f(x,  y)  trên miền 
{ }= ≤ ≤ ≤ ≤R (x,y|a x b,c(x) y d(x)  như hình vẽ. Ta cần tính tích phân: 
⎧ ⎫⎪ ⎪= = ⎨ ⎬⎪ ⎪⎩ ⎭∫∫ ∫∫
b
d(x)
R c(x)
a
J f(x,y)dxdy f(x,y)dy dx  
Công thức gần đúng của tích phân là: 
  [ ]
= =
= ∑ ∑m ni j i i ,j
i 1 j 1
J a,b,c(x),d(x) w v f(x ,y )  
với các trọng số wi, vj tuỳ thuộc vào cáh tính tích phân hàm một biến.  
Ta xây dựng hàm int2simpson() để tính tích phân kép bằng công thức 
Simpson.  
function J = int2simpson(f, a, b, c, d, m, n) 
%  tich phan kep  cua ham  f(x,yen mien R =  {(x,y)|a <= x <= b,  c(x) <= y <= 
d(x)} 
a 
y 
d(x)
hx0,y1 
b 
hx1 x0 
x 
c(x)
xm x2 x1  hx2  hxm 
hx0,y2 
hx1,y2 
hx1,y1 
 357
% dung quy tac Simpson 
if ceil(m) ~= floor(m) %be rong co dinh cua cac doan tren x 
    hx = m;  
    m = ceil((b ‐ a)/hx); 
end 
if mod(m, 2) ~= 0 
    m = m + 1;  
end 
hx = (b ‐ a)/m;  
i = 1:m+1;  
x = a + (i ‐ 1)*hx; 
if isnumeric(c) 
    cx(i) = c; %neu c la hang so 
else  
    cx(i) = feval(c,x(i)); %khi c la ham cua c x 
end 
if isnumeric(d) 
    dx(i) = d; %neu d la hang so 
else  
    dx(i) = feval(d,x(i)); %khi d la ham cua x 
end 
if ceil(n) ~= floor(n) %be rong co dinh theo y 
    hy = n;  
    nx(i) = ceil((dx(i)‐ cx(i))/hy); 
    ind = find(mod(nx(i),2) ~= 0);  
    nx(ind) = nx(ind) + 1; 
else %so khoang co dinh 
    if mod(n, 2) ~= 0 
        n = n + 1;  
    end 
    nx(i) = n; 
end 
for i = 1:m + 1 
    sx(i) = simpsonfxy(f,x(i),cx(i),dx(i),nx(i)); 
end 
kodd = 2:2:m;  
 358
keven = 3:2:m ‐ 1; %the set of odd/even indices 
J = hx/3*(sx(1) + sx(m + 1) + 4*sum(sx(kodd)) + 2*sum(sx(keven))); 
function Jf = simpsonfxy(f, x, c, d, n) 
%tich phan mot bien cua f(x,y) voi Ry = {c <= y <= d} 
if nargin < 5 
    n = 100;  
end 
if abs(d ‐ c)< eps | n <= 0 
    Jf = 0;  
    return;  
end 
if mod(n, 2) ~= 0 
    n = n + 1;  
end 
h = (d ‐ c)/n;  
y = c + [0:n]*h;  
fxy = feval(f,x,y); 
fxy(find(fxy == inf)) = realmax;  
fxy(find(fxy == ‐inf)) = ‐realmax; 
kodd = 2:2:n;  
keven = 3:2:n ‐ 1; 
Jf = h/3*(fxy(1) + fxy(n + 1) + 4*sum(fxy(kodd)) + 2*sum(fxy(keven))); 
Để tính thể tích của hình cầu ta dùng chương trình ctint2simp.m: 
clear all, clc 
%Tinh the tich hinh cau 
x = [‐1:0.05:1];  
y = [0:0.05:1];  
[X,Y] = meshgrid(x,y); 
f = inline(ʹsqrt(max(1 ‐ x.*x ‐ y.*y,0))ʹ,ʹxʹ,ʹyʹ); 
Z = f(X,Y);  
mesh(x,y,Z); 
a = ‐1;  
b = 1;  
c = 0;  
 359
d = inline(ʹsqrt(max(1 ‐ x.*x,0))ʹ,ʹxʹ); 
Vs1 = int2simpson(f, a, b, c, d, 100, 100)% so diem cho truoc 
error1 = Vs1 ‐ pi/3 
Vs2 = int2simpson(f, a, b, c, d, 0.01, 0.01) %be rong cac doan cho truoc 
error2 = Vs2 ‐ pi/3