FORTRAN数值求超越方程的根

方程在子程序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来求解更方便，毕竟上面的方法也是需要先作图观察得到猜测区间和初始值的

posted @ 2017-09-13 22:55  huangzc  阅读(1024)  评论(0)  编辑  收藏  举报