方程在子程序equ(x,y)中表示，有二分法bisection和newton-raphson法newraf，
可求多个根，对二分法需要输入多个区间，对newton-raphson法需要输入多个初始值
主程序（调用程序）
1 PROGRAM main 2 USE FindRoots 3 IMPLICIT NONE 4 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15) 5 INTEGER::n=3,i 6 REAL(iwp),ALLOCATABLE::sec(:,:),rt(:),x(:) 7 REAL(iwp)::tol=1.0e-8 8 ALLOCATE(sec(2,n),x(n),rt(n)) 9 !sec(:,1)=(/1.0_iwp,2.0_iwp/) 10 !sec(:,2)=(/2.0_iwp,3.0_iwp/) 11 !sec(:,3)=(/3.0_iwp,4.0_iwp/) 12 !sec(:,4)=(/7.0_iwp,8.0_iwp/) 13 !x=(/1.5_iwp,2.5_iwp,3.5_iwp,7.5_iwp/) 14 sec(:,1)=(/4.0_iwp,5.0_iwp/) 15 sec(:,2)=(/7.0_iwp,8.0_iwp/) 16 sec(:,3)=(/10.0_iwp,12.0_iwp/) 17 x=(/4.5_iwp,7.0_iwp,10.5_iwp/) 18 CALL bisection(n,sec,tol,rt) 19 DO i=1,n 20 WRITE(*,'(f15.8)')rt(i) 21 END DO 22 CALL NewRaf(n,x,tol,rt) 23 PRINT*, 24 DO i=1,n 25 WRITE(*,'(f15.8)')rt(i) 26 END DO 27 END PROGRAM
模块（子程序）
1 MODULE FindRoots 2 IMPLICIT NONE 3 CONTAINS 4 SUBROUTINE bisection(n,sec,tol,rt) 5 IMPLICIT NONE 6 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15) 7 INTEGER,INTENT(IN)::n 8 REAL(iwp),INTENT(IN)::sec(:,:) 9 REAL(iwp),INTENT(IN)::tol 10 REAL(iwp),INTENT(INOUT)::rt(:) 11 REAL(iwp)::x1,x2,x3,y1,y2,y3 12 INTEGER::i 13 DO i=1,n 14 x1=sec(1,i); x2=sec(2,i) 15 CALL equ(x1,y1); call equ(x2,y2) 16 IF(y1*y2<=0.0_iwp)THEN 17 DO WHILE(x2-x1>tol) 18 x3=(x2+x1)/2. 19 CALL equ(x3,y3) 20 IF(y3*y2<=0.0_iwp)THEN 21 x1=x3; y1=y3 22 ELSE 23 x2=x3; y2=y3 24 END IF 25 END DO 26 rt(i)=(x1+x2)/2. 27 ELSE 28 WRITE(*,'(a,i1,a)')"ROOT MIGHT NOT IN THE ",i," SECTION" 29 rt(i)=0 30 END IF 31 END DO 32 END SUBROUTINE 33 34 SUBROUTINE equ(x,y) 35 IMPLICIT NONE 36 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15) 37 REAL(iwp),INTENT(IN)::x 38 REAL(iwp),INTENT(INOUT)::y 39 REAL(iwp),PARAMETER::pi=4.*ATAN(1.0_iwp) 40 !y=2*EXP(-x*pi)-(1+EXP(-2*x*pi))*COS(x*pi) 41 y=2.0_iwp-2.0_iwp*COS(x)*COSH(x)+SIN(x)*SINH(x) 42 END SUBROUTINE 43 44 SUBROUTINE NewRaf(n,x,tol,rt) 45 IMPLICIT NONE 46 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15) 47 INTEGER,INTENT(IN)::n 48 REAL(iwp),PARAMETER::dx=1.0e-3 49 REAL(iwp),INTENT(IN)::x(:) 50 REAL(iwp),INTENT(IN)::tol 51 REAL(iwp),INTENT(INOUT)::rt(:) 52 REAL(iwp)::x1,x2,y1,yd,dy1 53 INTEGER::i 54 DO i=1,n 55 x1=x(i) 56 x2=x1 57 DO 58 x1=x2 59 CALL equ(x1,y1) 60 CALL equ(x1+dx,yd) 61 dy1=(yd-y1)/dx 62 x2=x1-y1/dy1 63 IF(ABS(x2-x1)<tol) EXIT 64 END DO 65 rt(i)=(x2+x1)/2. 66 END DO 67 END SUBROUTINE 68 END MODULE
当然了，事实上我觉得用mathematica通过画图和FindRoot来求解更方便，毕竟上面的方法也是需要先作图观察得到猜测区间和初始值的