Bài giảng Ứng dụng các phương pháp tính số - Chương 2: Giải phương trình đại số phi tuyến - Ngô Văn Thanh

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.2 Giải phương trình đại số phi tuyến. 2.1. Phương pháp khoanh vùng (Bracketing interval) 2.1.1. Khoanh vùng nghiệm bằng cách vẽ đồ thị các hàm (Bracketing a root by plotting the graph of functions). 2.1.2. Khoanh vùng nghiệm từ trong ra ngoài (Bracketing a root by outward from an initial interval). 2.1.3. Khoanh vùng nghiệm từ ngoài vào trong (Bracketing a root by inward on an initial interval). 2.2. Phương pháp phi đạo hàm (Methods without derivatives) 2.2.1. Phương pháp chia đôi (bisection method). 2.2.2. Phương pháp lặp đơn giản và thay thế liên tiếp (Simple iterative method with successive substitution). 2.2.3. Phương pháp cát tuyến (Secant method) 2.3. Phương pháp đạo hàm (Methods with derivatives) 2.3.1. Phương pháp lặp Newton (Newton iterative method). 2.3.2. Phương pháp lặp Newton cho hệ các phương trình không tuyến tính (Newton iterative method for the system of non-linear equations) @2009, Ngô Văn Thanh - Viện Vật Lý 2.1. Phương pháp khoanh vùng 2.1.1. Khoanh vùng nghiệm bằng cách vẽ đồ thị các hàm.  Một số kiểu nghiệm của một phương trình phi tuyến: @2009, Ngô Văn Thanh - Viện Vật Lý  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. @2009, Ngô Văn Thanh - Viện Vật Lý 2.1.2. 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. @2009, Ngô Văn Thanh - Viện Vật Lý 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. @2009, Ngô Văn Thanh - Viện Vật Lý 2.1.3. 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. @2009, Ngô Văn Thanh - Viện Vật Lý 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. @2009, Ngô Văn Thanh - Viện Vật Lý 2.2. Phương pháp phi đạo hàm 2.2.1 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ố). @2009, Ngô Văn Thanh - Viện Vật Lý 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' @2009, Ngô Văn Thanh - Viện Vật Lý 2.2.3 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: @2009, Ngô Văn Thanh - Viện Vật Lý 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 @2009, Ngô Văn Thanh - Viện Vật Lý 2.3. Phương pháp đạo hàm (Methods with derivatives) 2.3.1. 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 @2009, Ngô Văn Thanh - Viện Vật Lý  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. @2009, Ngô Văn Thanh - Viện Vật Lý 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 @2009, Ngô Văn Thanh - Viện Vật Lý 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 @2009, Ngô Văn Thanh - Viện Vật Lý  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. @2009, Ngô Văn Thanh - Viện Vật Lý 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) @2009, Ngô Văn Thanh - Viện Vật Lý 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 @2009, Ngô Văn Thanh - Viện Vật Lý 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 @2009, Ngô Văn Thanh - Viện Vật Lý 2.3.2. 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. @2009, Ngô Văn Thanh - Viện Vật Lý  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: @2009, Ngô Văn Thanh - Viện Vật Lý REAL, DIMENSION(:), INTENT(INOUT) :: x INTEGER, DIMENSION(size(x)) :: indx REAL, DIMENSION(size(x)) :: fvec,p REAL, DIMENSION(size(x),size(x)) :: fjac 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 @2009, Ngô Văn Thanh - Viện Vật Lý *********************************************************** INTEGER,FUNCTION assert_eq(n1,n2,n3,string) INTEGER, INTENT(IN) :: n1,n2,n3 CHARACTER(LEN=*),INTENT(IN) :: string * IF (n1==n2 .and. n2==n3) assert_eq = n1 ELSE write(6,*) 'nrerror: an assert_eq failed', string * END FUNCTION assert_eq ************************************************************* FUNCTION outerprod(a,b) REAL,DIMENSION(:), INTENT(IN) :: a,b REAL,DIMENSION(size(a),size(b)) :: outerprod * outerprod = spread(a,dim=2,ncopies=size(b))* & spread(b,dim=1,ncopies=size(a)) * END FUNCTION outerprod ************************************************************* @2009, Ngô Văn Thanh - Viện Vật Lý  Hàm SPREAD(A, DIM, N_copy)  DIM=1 : thêm các hàng, DIM=2 : thêm các cột.  Ví dụ A = (/2, 7/) SPREAD(A, 1, 4) SPREAD(A, 2, 4) @2009, Ngô Văn Thanh - Viện Vật Lý ******************************************************* SUBROUTINE swap(a,b) REAL, INTENT(IN) :: a,b REAL :: temp * temp = a a = b b = temp * END SUBROUTINE swap ******************************************************* INTEGER,FUNCTION imaxloc(arr) REAL, INTENT(IN) :: arr INTEGER, DIMENSION(1) :: imax * imax=maxloc(arr(:)) imaxloc=imax(1) * END FUNCTION imaxloc @2009, Ngô Văn Thanh - Viện Vật Lý