FORTRAN数值求超越方程的根

方程在子程序equ(x,y)中表示，有二分法bisection和newton-raphson法newraf，

可求多个根，对二分法需要输入多个区间，对newton-raphson法需要输入多个初始值

主程序（调用程序）

PROGRAM main
    USE FindRoots
    IMPLICIT NONE
    INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
    INTEGER::n=3,i
    REAL(iwp),ALLOCATABLE::sec(:,:),rt(:),x(:)
    REAL(iwp)::tol=1.0e-8
    ALLOCATE(sec(2,n),x(n),rt(n))
    !sec(:,1)=(/1.0_iwp,2.0_iwp/)
    !sec(:,2)=(/2.0_iwp,3.0_iwp/)
    !sec(:,3)=(/3.0_iwp,4.0_iwp/)
    !sec(:,4)=(/7.0_iwp,8.0_iwp/)
    !x=(/1.5_iwp,2.5_iwp,3.5_iwp,7.5_iwp/)
    sec(:,1)=(/4.0_iwp,5.0_iwp/)
    sec(:,2)=(/7.0_iwp,8.0_iwp/)
    sec(:,3)=(/10.0_iwp,12.0_iwp/)
    x=(/4.5_iwp,7.0_iwp,10.5_iwp/)
    CALL bisection(n,sec,tol,rt)
    DO i=1,n
        WRITE(*,'(f15.8)')rt(i)
    END DO
    CALL NewRaf(n,x,tol,rt)
    PRINT*,
    DO i=1,n
        WRITE(*,'(f15.8)')rt(i)
    END DO
END PROGRAM

模块（子程序）

MODULE FindRoots
    IMPLICIT NONE
    CONTAINS
    SUBROUTINE bisection(n,sec,tol,rt)
        IMPLICIT NONE
        INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
        INTEGER,INTENT(IN)::n
        REAL(iwp),INTENT(IN)::sec(:,:)
        REAL(iwp),INTENT(IN)::tol
        REAL(iwp),INTENT(INOUT)::rt(:)
        REAL(iwp)::x1,x2,x3,y1,y2,y3
        INTEGER::i
        DO i=1,n
            x1=sec(1,i); x2=sec(2,i)
            CALL equ(x1,y1); call equ(x2,y2)
            IF(y1*y2<=0.0_iwp)THEN
                DO WHILE(x2-x1>tol)
                    x3=(x2+x1)/2.
                    CALL equ(x3,y3)
                    IF(y3*y2<=0.0_iwp)THEN
                        x1=x3; y1=y3
                    ELSE
                        x2=x3; y2=y3
                    END IF
                END DO
                rt(i)=(x1+x2)/2.
            ELSE
                WRITE(*,'(a,i1,a)')"ROOT MIGHT NOT IN THE ",i," SECTION"
                rt(i)=0
            END IF
        END DO
    END SUBROUTINE

    SUBROUTINE equ(x,y)
        IMPLICIT NONE
        INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
        REAL(iwp),INTENT(IN)::x
        REAL(iwp),INTENT(INOUT)::y
        REAL(iwp),PARAMETER::pi=4.*ATAN(1.0_iwp)
        !y=2*EXP(-x*pi)-(1+EXP(-2*x*pi))*COS(x*pi)
        y=2.0_iwp-2.0_iwp*COS(x)*COSH(x)+SIN(x)*SINH(x)
    END SUBROUTINE

    SUBROUTINE NewRaf(n,x,tol,rt)
        IMPLICIT NONE
        INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
        INTEGER,INTENT(IN)::n
        REAL(iwp),PARAMETER::dx=1.0e-3
        REAL(iwp),INTENT(IN)::x(:)
        REAL(iwp),INTENT(IN)::tol
        REAL(iwp),INTENT(INOUT)::rt(:)
        REAL(iwp)::x1,x2,y1,yd,dy1
        INTEGER::i
        DO i=1,n
            x1=x(i)
            x2=x1
            DO
                x1=x2
                CALL equ(x1,y1)
                CALL equ(x1+dx,yd)
                dy1=(yd-y1)/dx
                x2=x1-y1/dy1
                IF(ABS(x2-x1)<tol) EXIT
            END DO
            rt(i)=(x2+x1)/2.
        END DO
    END SUBROUTINE
END MODULE

当然了，事实上我觉得用mathematica通过画图和FindRoot来求解更方便，毕竟上面的方法也是需要先作图观察得到猜测区间和初始值的