Phương pháp phi đạo hàm
Phương pháp chia đôi (bisection method).
Tìm nghiệm nhanh hơn phương pháp khoanh vùng nghiệm.
Có độ chính xác cao hơn.
Chỉ tìm được một nghiệm nào đó khoảng (a, b).
Xét khoảng (a, b) mà trong đó phương trình phi tuyến có nghiệm, tức là
Chọn điểm c là điểm giữa của (a, b).
Nếu như f (c) cùng dấu với f (a) thì thay khoảng (a, b) bằng (c, b).
Nếu như f (c) cùng dấu với f (b) thì thay khoảng (a, b) bằng (a, c).
Thực hiện qua trình lặp trên một số bước nào đó, hoặc khoảng chia đôi bé
hơn một thừa số cho trước (sai số)
49 trang |
Chia sẻ: nguyenlam99 | Lượt xem: 1024 | Lượt tải: 0
Bạn đang xem trước 20 trang tài liệu Tin học ứng dụng - Chương 1: Các phương pháp tính số, để xem tài liệu hoàn chỉnh bạn click vào nút DOWNLOAD ở trên
Ngo Van Thanh, IOP
11/2011
Phần II. Tin học ứng dụng
Chương 1: Các phương pháp tính số (LT: 5, TH:5)
Đạo hàm số
Tích phân số
Giải phương trình và hệ phương trình tuyến tính
Tìm nghiệm các phương trình phi tuyến
1.1 Đạo hàm số
Định nghĩa:
Đạo hàm dạng 2 điểm:
Đạo hàm dạng 3 điểm (End-Point):
Hoặc là (Mid-Point)
Đạo hàm dạng 5 điểm (Mid-Point):
Hoặc là (End-point)
Đạo hàm bậc 2:
Đạo hàm dạng 3 điểm (Mid-Point):
1.2 Tích phân số
Phương pháp Newton và phương pháp Simpson
Phương pháp FIT trực tiếp:
Các hệ số an được xác định bằng phương pháp khử Gauss
Biểu thức Newton-Cotes:
Chia nhỏ miền lấy tích phân (a, b) thành n + 1 điểm.
Đa thức nội suy Lagrange bậc n của hàm f :
Quy tắc hình thang (Trapezium).
Đa thức bậc 1,
Quy tắc Simpson.
Đa thức bậc 2,
Phương pháp tổ hợp MidPoint.
Chia miền tích phân thành n + 2 khoảng,
Phương pháp tổ hợp hình thang.
Chia miền tích phân thành n khoảng,
Phương pháp tổ hợp Simpson.
Chia miền tích phân thành n khoảng,
Phương pháp Monte-Carlo. Tích phân nhiều lớp
Tích Phân nhiều lớp:
Lấy tích phân theo phương pháp tổ hợp Simpson:
Chia miền tích phân thành n và m khoảng
Lấy tiếp tích phân theo phương pháp tổ hợp Simpson:
Trong đó:
Cuối cùng ta có:
Sai số:
1.3 Giải phương trình và hệ phương trình tuyến
tính
Phương pháp tính trực tiếp
Hệ M phương trình tuyến tính với N ẩn.
Biểu diễn dạng ma trận:
Phương pháp khử Gaussian
Xét hệ N phương trình tuyến tính với N ẩn số, dạng mở rộng ma trận:
Giải hệ phương trình bằng phương pháp khử theo hàng. Xét ví dụ đơn giản
Nhân hàng 1 cho 3 rồi cộng với hàng thứ 2, thay kết quả cho hàng 2.
Nhân hàng 1 cho -2 rồi cộng với hàng thứ 3, thay kết quả cho hàng 3.
Hàng 1 vừa rồi gọi là phương trình pivot hay hàng pivot, phần tử a11 gọi là phần tử
pivot.
Áp dụng phương pháp tương tự bắt đầu từ phương trình pivot thứ hai và phần tử
pivot thứ 2.
Hệ phương trình tương đương bây giờ có dạng
Giải lần lượt từ phương trình thứ 3 đến phương trình 1.
Trường hợp tổng quát:
Chọn hàng thứ nhất làm hàng pivot, ta có:
Thực hiện tiếp cho hàng pivot thứ hai:
Cuối cùng ta thu được ma trận chéo, trong đó các phần tử nằm phía dưới đường
chéo đều bằng 0.
Pivoting and scaling.
Nếu các phần tử pivot (các phần tử trên đường chéo) rất bé thì sẽ có sai số
do các phép tính làm tròn các con số sau dấu chấm thập phân.
Nếu các phần tử pivot bằng không thì không sử dụng được phương pháp
khử Gaussian.
Để giải quyết các vấn đề trên.
Nếu phần tử pivot bằng 0, thực hiện tráo đổi hàng pivot cho hàng kế tiếp.
Sử dụng các phương pháp scaling
Phương pháp scaling theo cột (partial pivoting).
Phương pháp scaling theo hàng (scaled partial pivoting).
Kết hợp cả hai phương pháp trên (full pivoting).
Cuối cùng là tìm các nghiệm bằng phương pháp thay thế ngược (backward
substitution)
Phương pháp scaling theo cột (partial pivoting):
Xét phần tử pivot akk, ta có
Tìm giá trị mki lớn nhất, tương đương với việc tìm phần tử có giá trị tuyệt đối
lớn nhất trong cột i: |aki|
Tráo đổi hàng pivot cho hàng tương ứng với
Ví dụ:
Phần tử |-3| lớn nhất trong các phần tử của cột 1.
Tráo đổi hàng thứ 2 cho hàng 1.
Tiếp tục thực hiện phương pháp này cho đến khi thu được ma trận rút gọn về
dạng chéo.
Phương pháp scaling theo hàng (scaled partial pivoting):
Xác định hệ số tỷ lệ là phần tử lớn nhất trong từng hàng
Tính tỷ số các phần tử trong cột k
Tìm hàng thứ p mà tỷ số này có giá trị lớn nhất.
Tráo đổi hàng pivot cho hàng thứ p đó.
Xét ví dụ:
Ta có:
Tính:
Vì a31/s3 lớn nhất nên tráo đổi hàng 1 cho hàng 3.
Phương pháp phân ly LU.
Giả thiết rằng ma trận A có thể viết dưới dạng tích của hai ma trận tam giác:
Suy ra
Ta có hai phương trình ma trận:
: giải phương trình này để tìm y
: giải phương trình này để tìm x
Giải phương trình:
Giải phương trình:
Phân ly LU để tìm hai ma trận L và U.
Các phần tử của các ma trận thoả mãn biểu thức:
Với i < j :
Với i = j :
Với i > j :
Đây là hệ N 2 phương trình với (N 2 + N ) ẩn số.
Thuật toán Crout:
Đầu tiên đặt N phần tử trên đường chéo của ma trận L bằng 1.
Chạy vòng lặp j = 1, 2,, N.
Bước 1: chạy tiếp vòng lặp trong i = 1,2, , j. Áp dụng điều kiện cho i < j ;
i = j và . Tính
Bước 2: chạy vòng lặp i = j + 1, j + 2, , N . Áp dụng điều kiện cho i > j. Tính
Kết hợp 2 ma trận L và U:
Chương trình đơn giản:
SUBROUTINE ludcmp(a,indx,d)
IMPLICIT NONE
REAL, DIMENSION(:,:), INTENT(INOUT) :: a
INTEGER, DIMENSION(:), INTENT(OUT) :: indx
REAL, INTENT(OUT) :: d
REAL, DIMENSION(size(a,1)) :: vv
REAL, PARAMETER :: TINY=1.0e-20 !A small number.
INTEGER :: j,n,imax
n=assert_eq(size(a,1),size(a,2),size(indx),‟ludcmp‟)
d=1.0 !No row interchanges yet.
vv=maxval(abs(a),dim=2)
if (any(vv == 0.0)) &
write(6,*) ‟ There is a row of zeros.‟
vv=1.0/vv
do j=1,n
imax=(j-1)+imaxloc(vv(j:n)*abs(a(j:n,j)))
if (j /= imax)
call swap(a(imax,:),a(j,:))
d = -d
vv(imax) = vv(j)
end if
indx(j) = imax
if (a(j,j) == 0.0) a(j,j)=TINY
a(j+1:n,j) = a(j+1:n,j)/a(j,j)
*
do i=1,j-1
a(i,j)=a(i,j)-sum(a(i,1:i-1)*a(1:i-1,j))
enddo
*
do i=j,n !This is i = j and i = j+1: ::N
a(i,j)=a(i,j)-sum(a(i,1:j-1)*a(1:j-1,j))
enddo
enddo
END SUBROUTINE ludcmp
Chương trình F90 kiểu parallel tối ưu:
SUBROUTINE ludcmp(a,indx,d)
IMPLICIT NONE
REAL, DIMENSION(:,:), INTENT(INOUT) :: a
INTEGER, DIMENSION(:), INTENT(OUT) :: indx
REAL, INTENT(OUT) :: d
REAL, DIMENSION(size(a,1)) :: vv
REAL, PARAMETER :: TINY=1.0e-20 !A small number.
INTEGER :: j,n,imax
n=assert_eq(size(a,1),size(a,2),size(indx),‟ludcmp‟)
d=1.0 !No row interchanges yet.
vv=maxval(abs(a),dim=2)
if (any(vv == 0.0)) &
write(6,*) ‟ There is a row of zeros.‟
vv=1.0/vv
do j=1,n
imax=(j-1) + imaxloc(vv(j:n)*abs(a(j:n,j)))
!Find the pivot row (only j:n elements).
if (j /= imax) then
!Do we need to interchange rows?
call swap(a(imax,:),a(j,:))
! Yes, do so...
d = -d !...and change the parity of d.
vv(imax)= vv(j)
!Also interchange the scale factor.
end if
indx(j) = imax
if (a(j,j) == 0.0) a(j,j) = TINY
a(j+1:n,j)=a(j+1:n,j)/a(j,j)
!Divide by the pivot element.
a(j+1:n,j+1:n)=a(j+1:n,j+1:n) - &
outerprod(a(j+1:n,j),a(j,j+1:n))
!Reduce remaining submatrix.
end do
END SUBROUTINE ludcmp
Nghịch đảo ma trận.
Xét hệ phương trình tuyến tính:
Nhân hai vế cho ma trận ngịch đảo A-1:
Nếu xác định được A-1 thì hệ phương trình hoàn toàn có thể tìm nghiệm dưới
dạng:
1.4 Tìm nghiệm các phương trình phi tuyến
Phương pháp khoanh vùng
Khoanh vùng nghiệm bằng cách vẽ đồ thị các hàm
Phương trình có nghiệm trong khoảng (a, b) nếu
liên tục trên khoảng (a, b).
Các bước thực hiện:
Vẽ đồ thị của hàm số bằng các phần mềm
như Gnuplot, Mathematica, Matlab
Dựa trên đồ thị, xác định khoảng (a, b) mà nghiệm nằm trong khoảng đó.
Xác định nghiệm gần đúng của phương trình x0.
Khoanh vùng nghiệm từ trong ra ngoài
Xét hai điểm bất kỳ
Tính tích
Nếu tích trên thì kết thúc chương trình tính.
Ngược lại, nếu
Nếu thì thay giá trị
Ngược lại: thì thay
Với b là thừa số tùy chọn.
Thực hiện các phép tính trên theo một số vòng lặp xác định.
INTEGER, PARAMETER :: NTRY=50
REAL, PARAMETER :: FACTOR=1.6
INTEGER :: j
REAL :: f1,f2
if (x1 == x2) RETURN
f1=func(x1)
f2=func(x2)
succes=.true.
do j=1,NTRY
if ((f1 > 0.0 .and. f2 < 0.0) .or. &
(f1 0.0)) RETURN
if (abs(f1) < abs(f2)) then
x1=x1+FACTOR*(x1-x2)
f1=func(x1)
else
x2=x2+FACTOR*(x2-x1)
f2=func(x2)
end if
end do
succes=.false.
Khoanh vùng nghiệm từ ngoài vào trong (cho phép khoanh vùng nhiều
nghiệm)
Xét hai điểm cho trước
Chia thành n khoảng và trong đó có tối đa là nb nghiệm.
Tính tích
Nếu tích trên đưa ra khoảng nghiệm
Thực hiện các phép tính trên cho n khoảng.
nbb=0
x=x1
dx=(x2-x1)/n
fp=fx(x)
do i=1,n !Loop over all intervals
x=x+dx
fc=fx(x)
if(fc*fp.le.0.) then
nbb=nbb+1
xb1(nbb)=x-dx
xb2(nbb)=x
if(nbb.eq.nb)goto 1
endif
fp=fc
enddo
1 continue
nb=nbb
END SUBROUTINE zbrak.
Phương pháp phi đạo hàm
Phương pháp chia đôi (bisection method).
Tìm nghiệm nhanh hơn phương pháp khoanh vùng nghiệm.
Có độ chính xác cao hơn.
Chỉ tìm được một nghiệm nào đó khoảng (a, b).
Xét khoảng (a, b) mà trong đó phương trình phi tuyến có nghiệm, tức là
Chọn điểm c là điểm giữa của (a, b).
Nếu như f (c) cùng dấu với f (a) thì thay khoảng (a, b) bằng (c, b).
Nếu như f (c) cùng dấu với f (b) thì thay khoảng (a, b) bằng (a, c).
Thực hiện qua trình lặp trên một số bước nào đó, hoặc khoảng chia đôi bé
hơn một thừa số cho trước (sai số).
fmid=func(x2)
f=func(x1)
if (f*fmid >= 0.0) write(6,*) 'rtbis:root must be bracketed„
if (f < 0.0) then
rtbis=x1
dx=x2-x1
else
rtbis=x2
dx=x1-x2
end if
do j=1,MAXIT !Bisection loop.
dx=dx*0.5
xmid=rtbis+dx
fmid=func(xmid)
if (fmid <= 0.0) rtbis = xmid
if (abs(dx) < xacc .or. fmid == 0.0) RETURN
end do
write(6,*) 'rtbis:too many bisections'
Phương pháp cát tuyến (Secant method).
Áp dụng cho các hàm trơn (smooth) ở gần nghiệm.
Tốc độ hội tụ nhanh hơn phương pháp bisection.
Xét khoảng (a, b) mà trong đó phương trình phi tuyến có nghiệm.
Chọn hai điểm ban đầu p0 = a, p1 = b.
Phương trình đường thẳng đi qua (p0, f(p0)) và (p1, f(p1))
Giao điểm với trục hoành tại (p2, 0):
suy ra
Tổng quát:
fl=func(x1)
f =func(x2)
if (abs(fl) < abs(f)) then
rtsec=x1
xl=x2
call swap(fl,f)
else
xl=x1
rtsec=x2
end if
do j=1,MAXIT !Secant loop.
dx=(xl-rtsec)*f/(f-fl)
xl=rtsec
fl=f
rtsec=rtsec+dx
f=func(rtsec)
if (abs(dx) < xacc .or. f == 0.0) RETURN ! Convergence.
end do
Phương pháp đạo hàm (Methods with derivatives)
Phương pháp lặp Newton (Newton iterative method).
Khai triển chuỗi Taylor:
Xét: suy ra
Hệ số góc của phương trình f(x):
Suy ra
Nghiệm của phương trình là khi
Nếu hàm số f(x) không tính được đạo hàm bằng giải tích thì phải tính f’ (x)
theo phương pháp gần đúng tại mỗi vòng lặp.
Thuật toán:
Run_newton(f; f‟; x0;N; tol)
(1) đặt x = x0; n = 0
(2) while n <= N
(3) tính fx f(x)
(4) if (|fx| < tol)
(5) return x
(6) tính fpx f‟(x)
(7) if |fpx| < tol
(8) return x
(9) Gán x = x – fx/fpx
(10) Gán n = n + 1
(11) end while
Chương trình:
INTEGER, PARAMETER :: MAXIT=20
INTEGER :: j
REAL :: df,dx,f
rtnewt=0.5*(x1+x2) Initial guess.
do j=1,MAXIT
call funcd(rtnewt,f,df)
dx=f/df
rtnewt=rtnewt-dx
if ((x1-rtnewt)*(rtnewt-x2) < 0.0)&
write(6,*) 'rtnewt:values jumped out of brackets'
if (abs(dx) < xacc) RETURN ! Convergence.
end do
write(6,*) 'rtnewt exceeded maximum iterations'
END FUNCTION rtnewt
Một số trường hợp mà phương pháp lặp Newton sẽ tính sai.
Kết hợp bisection và Newton-Raphson. Phương pháp bisection được áp dụng
cho trường hợp phương pháp Newton hội tụ chậm hoặc nghiệm tìm được
vượt ra ngoài khoảng nghiệm.
Chương trình:
call funcd(x1,fl,df)
call funcd(x2,fh,df)
if ((fl > 0.0 .and. fh > 0.0) .or. &
(fl < 0.0 .and. fh < 0.0)) &
write(6,*)'root must be bracketed in rtsafe'
if (fl == 0.0) then
rtsafe=x1
RETURN
else if (fh == 0.0) then
rtsafe=x2
RETURN
else if (fl < 0.0)
xl=x1
xh=x2
else
xh=x1
xl=x2
end if
*
rtsafe=0.5*(x1+x2) !Initialize the guess for root,
dxold=abs(x2-x1) !the “stepsize before last,”
dx=dxold !and the last step.
call funcd(rtsafe,f,df)
Chương trình:
do i=1,ntrial
call usrfun(x,fvec,fjac)
! User subroutine supplies function values at
! x in fvec and Jacobian matrix in fjac.
if (sum(abs(fvec)) <= tolf) RETURN
!Check function convergence.
!Right-hand side of linear equations.
p=-fvec
call ludcmp(fjac,indx,d)
!Solve linear equations using LU decomposition
call lubksb(fjac,indx,p)
!Update solution.
x=x+p
if (sum(abs(p)) <= tolx) RETURN
!Check root convergence.
end do
Phương pháp lặp Newton cho hệ các phương trình không tuyến tính
Xét hệ hai phương trình không tuyến tính bất kỳ:
Giao điểm của các đường mức chính là nghiệm của hệ phương trình.
Xét hệ N phương trình không tuyến tính
Khai triển chuỗi Taylor cho hệ phương trình trên
Ma trận các đạo hàm riêng chính là ma trận Jacobian
Suy ra:
Bỏ qua số hạng bậc cao, đặt , ta có
Giải phương trình trên để tìm bằng phương pháp phân ly LU. Nghiệm gần
đúng của hệ có dạng:
do j=1,MAXIT !Loop over allowed iterations.
if (((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f) > 0.0 .or. &
abs(2.0*f) > abs(dxold*df) ) then
! Bisect if Newton out of range, or not decreasing fast enough.
dxold=dx
dx=0.5*(xh-xl)
rtsafe=xl+dx
if (xl == rtsafe) RETURN !Change in root is negligible.
else !Newton step acceptable. Take it.
dxold=dx
dx=f/df
temp=rtsafe
rtsafe=rtsafe-dx
if (temp == rtsafe) RETURN
end if
if (abs(dx) < xacc) RETURN !Convergence criterion.
call funcd(rtsafe,f,df)
if (f < 0.0) then !Maintain the bracket on the root.
xl=rtsafe
else
xh=rtsafe
end if
end do
write(6,*) 'rtsafe:exceeded maximum iterations'
END FUNCTION rtsafe
Các file đính kèm theo tài liệu này:
- bai_giang_p2_01_1118.pdf