Bài giảng Ứng dụng các phương pháp tính số - Chương 4: Các phương pháp cho hàm rời rạc - Ngô Văn Thanh
4.1. Đa thức nội suy (Polynomial interpolation) 4.1.1. Nội suy Lagrangian. 4.1.2. Nội suy Newton. 4.1.3. Nội suy bậc 3 (Cubic spline interpolation). 4.2. Phương pháp bình phương tối thiểu. 4.2.1. Làm khớp dữ liệu cho các hàm tuyến tính. 4.2.2. Làm khớp dữ liệu cho các hàm phi tuyến.
Bạn đang xem trước 20 trang tài liệu Bài giảng Ứng dụng các phương pháp tính số - Chương 4: Các phương pháp cho hàm rời rạc - Ngô Văn Thanh, để xem tài liệu hoàn chỉnh bạn click vào nút DOWNLOAD ở trên
TS. Ngô Văn Thanh,
Viện Vật lý.
Cao học vật lý – chuyên ngành Vật lý lý thuyết.
Chương 4. Các phương pháp cho hàm rời rạc.
4.1. Đa thức nội suy (Polynomial interpolation)
4.1.1. Nội suy Lagrangian.
4.1.2. Nội suy Newton.
4.1.3. Nội suy bậc 3 (Cubic spline interpolation).
4.2. Phương pháp bình phương tối thiểu.
4.2.1. Làm khớp dữ liệu cho các hàm tuyến tính.
4.2.2. Làm khớp dữ liệu cho các hàm phi tuyến.
@2009, Ngô Văn Thanh - Viện Vật Lý
4.1. Đa thức nội suy (Polynomial interpolation).
4.1.1. Nội suy Lagrangian.
¾ Xét N điểm dữ liệu (xi, yi), i =1, N: Đa thức nội suy bậc N - 1 có dạng:
¾ Đa thức nội suy chỉ có thể tính được các giá trị của
@2009, Ngô Văn Thanh - Viện Vật Lý
@2009, Ngô Văn Thanh - Viện Vật Lý
¾ Thuật toán Neville :
P1, P2, , PN là đa thức bậc 0 qua mỗi một điểm
P12, P23, , P(N-1)N là đa thức bậc 1 qua từng cặp điểm
Đa thức bậc N – 1 : P123N
Ví dụ với N = 4: xây dựng dạng bảng các đa thức nội suy
x1:
x2:
x3:
x4:
Biểu thức liên hệ giữa hai đa thức bậc thấp hơn gần nhất:
@2009, Ngô Văn Thanh - Viện Vật Lý
¾ Thuật toán Neville (tiếp theo) :
Định nghĩa :
@2009, Ngô Văn Thanh - Viện Vật Lý
¾ Chương trình:
*********************************************************
INTEGER,FUNCTION iminloc(arr)
REAL, INTENT(IN) :: arr
INTEGER, DIMENSION(1) :: imin
imin=minloc(arr(:))
iminloc=imin(1)
END FUNCTION iminloc
*********************************************************
SUBROUTINE polint(xa,ya,x,y,dy)
IMPLICIT NONE
REAL, DIMENSION(:), INTENT(IN) :: xa,ya
REAL, INTENT(IN) :: x
REAL, INTENT(OUT) :: y,dy
INTEGER :: m,n,ns
REAL, DIMENSION(size(xa)) :: c,d,den,ho
n=size(xa)
@2009, Ngô Văn Thanh - Viện Vật Lý
c=ya !Initialize the tableau of c's and d's.
d=ya
ho=xa-x
ns=iminloc(abs(x-xa)) !Find index ns of closest table entry.
y=ya(ns) !This is the initial approximation to y.
ns=ns-1
do m=1,n-1 !For each column of the tableau,
den(1:n-m)=ho(1:n-m)-ho(1+m:n)
if(any(den(1:n-m)==0.0)) write(6,*)'calculation failure'
den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m)
d(1:n-m)=ho(1+m:n)*den(1:n-m)
c(1:n-m)=ho(1:n-m)*den(1:n-m)
if (2*ns < n-m) then
dy=c(ns+1)
else
dy=d(ns)
ns=ns-1
endif
y=y+dy
enddo
END SUBROUTINE polint
@2009, Ngô Văn Thanh - Viện Vật Lý
4.1.2. Nội suy Newton
¾ Biểu thức Newton
¾ Xác định hệ số c người ta sử dụng biểu thức:
¾ Dạng tính kiểu truy hồi
¾ Dạng tổng quát:
@2009, Ngô Văn Thanh - Viện Vật Lý
4.1.3. Cubic spline interpolation
¾ Xét N điểm dữ liệu (xi, yi), i =1, N.
¾ Xét điểm x nằm trong khoảng hai điểm
¾ Ta có:
¾ Nếu ta có các giá trị đạo hàm bậc hai thì ta có thể viết
¾ Với
¾ Các đạo hàm:
@2009, Ngô Văn Thanh - Viện Vật Lý
¾ Chương trình:
*********************************************************
SUBROUTINE spline(x,y,yp1,ypn,y2)
IMPLICIT NONE
REAL, DIMENSION(:), INTENT(IN) :: x,y
REAL, INTENT(IN) :: yp1,ypn
REAL, DIMENSION(:), INTENT(OUT) :: y2
INTEGER :: n
REAL, DIMENSION(size(x)) :: a,b,c,r
n=size(x)
c(1:n-1)=x(2:n)-x(1:n-1) !Set up the tridiagonal equations.
r(1:n-1)=6.0*((y(2:n)-y(1:n-1))/c(1:n-1))
r(2:n-1)=r(2:n-1)-r(1:n-2)
a(2:n-1)=c(1:n-2)
b(2:n-1)=2.0*(c(2:n-1)+a(2:n-1))
b(1)=1.0
b(n)=1.0
@2009, Ngô Văn Thanh - Viện Vật Lý
if (yp1 > 0.99e30) then
r(1)=0.0
c(1)=0.0
else
r(1)=(3.0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
c(1)=0.5
end if
if (ypn > 0.99e30) then
r(n)=0.0
a(n)=0.0
else
r(n)=(-3.0/(x(n)-x(n-1))) &
*((y(n)-y(n-1))/(x(n)-x(n-1))-ypn)
a(n)=0.5
end if
call tridag(a(2:n),b(1:n),c(1:n-1),r(1:n),y2(1:n))
END SUBROUTINE spline
@2009, Ngô Văn Thanh - Viện Vật Lý
SUBROUTINE tridag(a,b,c,r,u)
bet=b(1)
if(bet==0.) write(6,*) 'Error 1 in tridag'
allocate(gam(n))
u(1)=r(1)/bet
do j=2,n
gam(j)=c(j-1)/bet
bet=b(j)-a(j-1)*gam(j)
if(bet==0.) write(6,*) 'Error 2 in tridag'
u(j)=(r(j)-a(j-1)*u(j-1))/bet
enddo
do j=n-1,1,-1
u(j)=u(j)-gam(j+1)*u(j+1)
enddo
deallocate(gam)
@2009, Ngô Văn Thanh - Viện Vật Lý
FUNCTION splint(xa,ya,y2a,x)
IMPLICIT NONE
REAL, DIMENSION(:), INTENT(IN) :: xa,ya,y2a
REAL, INTENT(IN) :: x
REAL :: splint
INTEGER :: khi,klo,n
REAL :: a,b,h
n=size(xa)
klo=max(min(locate(xa,x),n-1),1)
khi=klo+1
h=xa(khi)-xa(klo)
if (h == 0.0) write(6,*) 'The xa must be distinct'
a=(xa(khi)-x)/h
b=(x-xa(klo))/h
splint=a*ya(klo)+b*ya(khi)
& +((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.0
END FUNCTION splint
@2009, Ngô Văn Thanh - Viện Vật Lý
4.2. Phương pháp bình phương tối thiểu.
¾ Xét N điểm dữ liệu (xi, yi), i =1, N với M tham số hiệu chỉnh a1, a2, aM.
¾ Hàm làm khớp có dạng:
¾ Định nghĩa hàm bình phương tối thiểu
¾ Trong đó độ lệch chuẩn
¾ Cực tiểu của hàm bình phương tối thiểu
@2009, Ngô Văn Thanh - Viện Vật Lý
4.2.1. Làm khớp dữ liệu cho các hàm tuyến tính.
Xét N điểm dữ liệu (xi, yi), i =1, N: Làm khớp dữ liệu cho một hàm tuyến tính
bậc nhất có dạng:
¾ Hàm bình phương tối thiểu
¾ Cực tiểu của hàm bình phương tối thiểu
¾ Định nghĩa các tổng
@2009, Ngô Văn Thanh - Viện Vật Lý
¾ Nghiệm của hệ là
¾ Sai số:
Định nghĩa
Viết lại nghiệm a, b và các sai số:
@2009, Ngô Văn Thanh - Viện Vật Lý
Xét N điểm dữ liệu (xi, yi), i =1, N: Làm khớp dữ liệu cho một hàm tuyến tính
có dạng:
Trong đó Xk(x) gọi là hàm cơ sở
¾ Hàm bình phương tối thiểu
¾ Định nghĩa ma trận design A
¾ Các hệ số làm khớp cũng là một ma trận a
@2009, Ngô Văn Thanh - Viện Vật Lý
¾ Lấy đạo hàm của hàm bình phương tối thiểu
¾ Viết lại dạng phương trình ma trận
Trong đó
¾ Phương trình ma trận có thể viết dưới dạng:
@2009, Ngô Văn Thanh - Viện Vật Lý
¾ Các hệ số a
¾ Độ lệch:
@2009, Ngô Văn Thanh - Viện Vật Lý
SUBROUTINE lfit(x,y,sig,a,maska,covar,chisq,funcs)
IMPLICIT NONE
REAL, DIMENSION(:), INTENT(IN) :: x,y,sig
REAL, DIMENSION(:), INTENT(INOUT) :: a
LOGICAL, DIMENSION(:), INTENT(IN) :: maska
REAL, DIMENSION(:,:), INTENT(INOUT) :: covar
REAL, INTENT(OUT) :: chisq
INTERFACE
SUBROUTINE funcs(x,arr)
END SUBROUTINE funcs
END INTERFACE
INTEGER :: i,j,k,l,ma,mfit,n
REAL :: sig2i,wt,ym
REAL, DIMENSION(size(maska)) :: afunc
REAL, DIMENSION(size(maska),1) :: beta
n=size(x)
ma=size(maska)
mfit=count(maska)
@2009, Ngô Văn Thanh - Viện Vật Lý
if (mfit == 0) write(6,*) 'no parameters to be fitted'
covar(1:mfit,1:mfit)=0.0 !Initialize the matrix.
beta(1:mfit,1)=0.0
do i=1,n
call funcs(x(i),afunc) !the normal equations.
ym = y(i)
if (mfit< ma) ym=ym-sum(a(1:ma) &
*afunc(1:ma),mask=.not. maska)
sig2i=1.0/sig(i)**2; j = 0
do l=1,ma
if (maska(l)) then
j=j+1
wt=afunc(l)*sig2i
k=count(maska(1:l))
covar(j,1:k)=covar(j,1:k)&
+wt*pack(afunc(1:l),maska(1:l))
beta(j,1)=beta(j,1)+ym*wt
end if
enddo
enddo
@2009, Ngô Văn Thanh - Viện Vật Lý
call diagmult(covar(1:mfit,1:mfit),0.5)
covar(1:mfit,1:mfit)= covar(1:mfit,1:mfit) &
+transpose(covar(1:mfit,1:mfit))
call gaussj(covar(1:mfit,1:mfit),beta(1:mfit,1:1
a(1:ma)=unpack(beta(1:ma,1),maska,a(1:ma))
chisq=0.0 !Evaluate \chi^2 of the fit.
do i=1,n
call funcs(x(i),afunc)
chisq=chisq+((y(i)
dot_product(a(1:ma),afunc(1:ma)))/sig(i))**2
end do
call covsrt(covar,maska)
END SUBROUTINE lfit
@2009, Ngô Văn Thanh - Viện Vật Lý
4.2.2. Làm khớp dữ liệu cho các hàm phi tuyến.
@2009, Ngô Văn Thanh - Viện Vật Lý
