应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » 【求助】求教fortran 90高手！！
141/2返回列表12下一页
查看: 807  |  回复: 13
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前主题已经存档。

jianchaoyv

金虫 (小有名气)

[交流] 【求助】求教fortran 90高手！！

一个fortran 90程序如下：
program D3R1
!driver for routine TRAPZD
real,parameter::nmax=15,pio2=1.5707963
external Func,Fint
real::A=0.0,B=pio2
write(*,'(1x,a)') 'integral of Func with 2^(n-1) points'
write(*,'(1x,a,f10.6)')'actual value of integral is',&
               Fint(B)-Fint(A)
write(*,'(1x,t7,a,t16,a)')'n','Approx.integral'
do i=1,nmax
   call TRAPZD(Func,A,B,s,i)
   write(*,'(1x,i6,f20.6)')i,s
end do
end program D3R1
Function Func(x)
    Func=(x**2)*(x**2-2.0)*sin(x)
end function Func
function Fint(x)
!integral of Func
Fint=4.0*x*((x**2)-7.0)*sin(x)-((x**4)-14.0*(x**2)+28.0)*cos(x)
end function Fint



subroutine TRAPZD(Func,A,B,s,n)
integer,parameter::k1=selected_int_kind(9)
integer(kind=k1)::tnm
integer,parameter::k2=selected_real_kind(8,13)
real(kind=k2)::del,sum,x
if(n==1)then
  s=0.5*(b-a)*(Func(a)+Func(b))
  else
  tnm=2**(n-1)
  del=(b-a)/tnm
  x=A
  sum=0.0
  do j=2,tnm
    x=x+del
    sum=sum+Func(x)
  end do
   s=0.5*(Func(A)+Func(B)+2.0*sum)*del
end if
end subroutine TRAPZD
运行后结果：
integral of Func with 2^(n-1) points
actual value of integral is -0.479158
      n        Approx.integral
      1            0.905772
      2            6.166170
      3 -Infinity
      4 -Infinity
      5 -Infinity
      6 NaN
      7 NaN
      8 NaN
      9 NaN
     10 NaN

run-time error M6201: MATH
- sin: DOMAIN error
Image              PC        Routine            Line        Source
pp.exe             00407039  Unknown               Unknown  Unknown
pp.exe             00406E6B  Unknown               Unknown  Unknown
pp.exe             00406FF1  Unknown               Unknown  Unknown
pp.exe             00409268  Unknown               Unknown  Unknown
pp.exe             00429F90  Unknown               Unknown  Unknown
pp.exe             00426CFC  Unknown               Unknown  Unknown
pp.exe             004012A6  Unknown               Unknown  Unknown
pp.exe             00401451  Unknown               Unknown  Unknown
pp.exe             004011E5  Unknown               Unknown  Unknown
pp.exe             00433E29  Unknown               Unknown  Unknown
pp.exe             00426234  Unknown               Unknown  Unknown
kernel32.dll       7C817067  Unknown               Unknown  Unknown

Incrementally linked image--PC correlation disabled.
问题：
（1）怎样解决运行后出现的情况（如上）
（2）语句：integer,parameter::k1=selected_int_kind(9)中若将9改为13将会出现返回值k1=-1，怎样解决？
(3)将语句:integer,parameter::k2=selected_real_kind(8,13)中(8,13改为(6,9)后运行结果:
integral of Func with 2^(n-1) points
actual value of integral is -0.479158
      n        Approx.integral
      1            0.905772
      2           -0.020945
      3           -0.361461
      4           -0.449584
      5           -0.471756
      6           -0.477307
      7           -0.478697
      8           -0.479042
      9           -0.479126
     10           -0.479146
     11           -0.479151
     12           -0.479152
     13           -0.479153
     14           -0.479411
     15           -0.479676
Press any key to continue
为何在 14           -0.479411
          15           -0.479676 处变大那么多?不是n越大最后的值越接近真实值吗?

[ Last edited by jianchaoyv on 2009-4-10 at 16:57 ]
回复此楼

» 猜你喜欢
1楼 2009-04-10 15:43:36
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

老虎大王

木虫 (著名写手)
★ ★ ★
kuhailangyu(金币+3,VIP+0):替虫友致谢！ 4-10 20:13
integer,parameter::k1=selected_int_kind(9)
integer(kind=k1)::tnm
integer,parameter::k2=selected_real_kind(8,13)
real(kind=k2)::del,sum,x

这几句表明了tnm, del,sum, x的类型（后三个是双精度的，8字节的浮点数，相当于real*8）

你调用了Func(x)，但在Func函数中没有说明x的类型，默认则为4字节浮点数，相当于real*4。这样类型不匹配，可能是造成cos或Sin函数出现问题的原因。
你后来的修改，等于是把参数都改成单精度，这样损失了精度，计算不准确。

简单的方法是在Func(x)中加上一行说明：
real*8::x
或者为了保持统一，加上两行说明：
integer,parameter::k2=selected_real_kind(8,13)
real（kind=k2)::x  也行。

这样运行程序，就好多了。结果为：
integral of Func with 2^(n-1) points
actual value of integral is -0.479158
      n        Approx.integral
      1            0.000000
      2           -0.473831
      3           -0.587905
      4           -0.562805
      5           -0.528367
      6           -0.505613
      7           -0.492849
      8           -0.486120
      9           -0.482668
     10           -0.480921
     11           -0.480042
     12           -0.479601
     13           -0.479380
     14           -0.479269
     15           -0.479214

这还可能不是最好的，因为还有参数a,  b 的问题。我建议你最好把在主程序和子程序中所有的浮点数（包括a，b 等）都说明成双精度数，可能会更好。

注：我没有考查你的算法，不知道结果对不对，但是从编程来讲，数据类型要尽量统一是没错的。呵呵。

[ Last edited by 老虎大王 on 2009-4-10 at 20:14 ]
赞一下(10人) 回复此楼
2楼2009-04-10 20:12:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

老虎大王

木虫 (著名写手)
★ ★ ★
jianchaoyv(金币+3,VIP+0):问题没全部解决,但是非常感谢你这位绝对是高手了!!! 4-10 21:25
索性我给你弄了：
程序改为：
program D3R1
!driver for routine TRAPZD
integer,parameter::nmax=15                ！这里有改动
real*8,parameter::pio2=1.5707963       ！这里有改，请你自己再把PiO2写准确一点，可能会更好
external Func,Fint
real*8 ::A=0.0,B=pio2           ！这里有改
write(2,'(1x,a)') 'integral of Func with 2^(n-1) points'
write(2,'(1x,a,f10.6)')'actual value of integral is',&
               Fint(B)-Fint(A)
write(2,'(1x,t7,a,t16,a)')'n','Approx.integral'
do i=1,nmax
   call TRAPZD(Func,A,B,s,i)
   write(2,'(1x,i6,f20.6)')i,s
end do
end program D3R1


Function Func(x)
real*8::x                      ！这里改了
real*8::Func                        ！这里改了
    Func=(x**2)*(x**2-2.0)*sin(x)
end function Func


function Fint(x)
!integral of Func
real*8::x           ！这里改了
real*8::Fint       ！这里改了
Fint=4.0*x*((x**2)-7.0)*sin(x)-((x**4)-14.0*(x**2)+28.0)*cos(x)
end function Fint



subroutine TRAPZD(Func,A,B,s,n)
integer,parameter::k1=selected_int_kind(9)
integer(kind=k1)::tnm   
integer,parameter::k2=selected_real_kind(8,13)
real*8::del,sum,x ,a,b  !(kind=k2)                       ！这里虽然改了，其实是一样的

if(n==1)then
  s=0.5*(b-a)*(Func(a)+Func(b))
  else
  tnm=2**(n-1)
  del=(b-a)/tnm
  x=A
  sum=0.0
  do j=2,tnm
    x=x+del
    sum=sum+Func(x)
  end do
   s=0.5*(Func(A)+Func(B)+2.0*sum)*del
end if
end subroutine TRAPZD


结果为：

integral of Func with 2^(n-1) points
actual value of integral is -0.479158
      n        Approx.integral
      1            0.905773
      2           -0.020945
      3           -0.361461
      4           -0.449584
      5           -0.471756
      6           -0.477308
      7           -0.478696
      8           -0.479043
      9           -0.479130
     10           -0.479152
     11           -0.479157
     12           -0.479158
     13           -0.479159
     14           -0.479159
     15           -0.479159
赞一下(2人) 回复此楼
3楼2009-04-10 20:24:50
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

老虎大王

木虫 (著名写手)
呵呵呵。
回复此楼
4楼2009-04-10 20:25:07
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

老虎大王

木虫 (著名写手)
★ ★ ★ ★ ★ ★ ★ ★
kuhailangyu(金币+1,VIP+0):追着感谢！ 4-10 21:59
jianchaoyv(金币+7,VIP+0):感谢帮忙,交个朋友吧,以后多多探讨! 4-11 09:03
?什么地方没解决？

你还可以在主程序中这样写：
real*8::a,b,pio2
...
a=0.d0
pio2=acos(a)
b=pio2
这样，Pi/2就有很高的精度。
若把输出改成输出9位小数，我算了，收敛到-0.479158789，很不错的啊。你写的actual value of integral is -0.479158，非常接近啊。
有什么问题说出来，咱再研究。不过数学上的算法我不管啊。呵呵。

[ Last edited by 老虎大王 on 2009-4-10 at 21:51 ]
赞一下(11人) 回复此楼
5楼2009-04-10 21:50:00
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

jianchaoyv

金虫 (小有名气)
昨天晚上我按你的提示进行调试,发现有2个警告说call TRAPZD(Func,A,B,s,i)
时类型有问题,今天我再运行就可以了,哈哈...十分感谢你的帮忙!请问你是什么专业的?我是搞物理计算的,还有专业问题想请教你....
赞一下 回复此楼
6楼2009-04-11 09:20:16
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

老虎大王

木虫 (著名写手)
呵呵。是的，我忘记了。TRAPZD中的s也要搞成双精度。
我是搞材料计算的。平时主要用Fortran编程。
赞一下 回复此楼
7楼2009-04-11 11:10:41
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

jianchaoyv

金虫 (小有名气)
再向你请教个物理问题
问题如下：
     有十个氢原子在300K温度下其初速度满足麦克斯韦速率分布，即正比与exp[-mv**2/2kT]，用fortran语言编程求其初速度。
请指教怎样编程
赞一下 回复此楼
8楼2009-04-11 15:49:38
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

老虎大王

木虫 (著名写手)
你有没有《computer simulation of liquids》这本书？本论坛好像可以下载的。可以阅读书后附录G3。书中有一个生成初始速度的程序。我列在下面，你可以参考。
另外你可以参考一些编程书，如徐士良的《fortran常用算法程序集》，里面介绍了如何生成正态分布的随机数，你按自己的要求改写一下就可以了。

********************************************************************************
** FICHE F.24.  INITIAL VELOCITY DISTRIBUTION                                 **
** This FORTRAN code is intended to illustrate points made in the text.       **
** To our knowledge it works correctly.  However it is the responsibility of  **
** the user to test it, if it is to be used in a research application.        **
********************************************************************************

C    *******************************************************************
C    ** CENTRE OF MASS AND ANGULAR VELOCITIES FOR LINEAR MOLECULES    **
C    **                                                               **
C    ** PRINCIPAL VARIABLES:                                          **
C    **                                                               **
C    ** INTEGER N                   THE NUMBER OF MOLECULES           **
C    ** REAL    RX(N),RY(N),RZ(N)   POSITIONS                         **
C    ** REAL    VX(N),VY(N),VZ(N)   VELOCITIES                        **
C    ** REAL    EX(N),EY(N),EZ(N)   ORIENTATIONS                      **
C    ** REAL    OX(N),OY(N),OZ(N)   SPACE-FIXED ANGULAR VELOCITIES    **
C    ** REAL    TEMP                REDUCED TEMPERATURE               **
C    ** REAL    INERT               REDUCED MOMENT OF INERTIA         **
C    **                                                               **
C    ** SUPPLIED ROUTINES:                                            **
C    **                                                               **
C    ** SUBROUTINE COMVEL ( TEMP )                                    **
C    **    SETS THE CENTRE OF MASS VELOCITIES FOR A CONFIGURATION OF  **
C    **    LINEAR MOLECULES AT A GIVEN TEMPERATURE.                   **
C    ** SUBROUTINE ANGVEL ( TEMP, INERT )                             **
C    **    SETS THE ANGULAR VELOCITIES FOR A CONFIGURATION OF LINEAR  **
C    **    MOLECULES AT A GIVEN TEMPERATURE.                          **
C    ** REAL FUNCTION RANF ( DUMMY )                                  **
C    **    RETURNS A UNIFORM RANDOM VARIATE ON THE RANGE ZERO TO ONE  **
C    ** REAL FUNCTION GAUSS ( DUMMY )                                 **
C    **    RETURNS A UNIFORM RANDOM NORMAL VARIATE FROM A             **
C    **    DISTRIBUTION WITH ZERO MEAN AND UNIT VARIANCE.             **
C    **                                                               **
C    ** UNITS:                                                        **
C    **                                                               **
C    ** WE ASSUME UNIT MOLECULAR MASS AND EMPLOY LENNARD-JONES UNITS  **
C    **       PROPERTY                      UNITS                     **
C    **       RX, RY, RZ           (EPSILON/M)**(1.0/2.0)             **
C    **       OX, OY, OZ           (EPSILON/M*SIGMA**2)**(1.0/2.0)    **
C    **       INERT                 M*SIGMA**2                        **
C    *******************************************************************

        SUBROUTINE COMVEL ( TEMP )

        COMMON / BLOCK1 / RX, RY, RZ, EX, EY, EZ,
     :                    VX, VY, VZ, OX, OY, OZ

C    *******************************************************************
C    ** TRANSLATIONAL VELOCITIES FROM MAXWELL-BOLTZMANN DISTRIBUTION  **
C    **                                                               **
C    ** THE DISTRIBUTION IS DETERMINED BY TEMPERATURE AND (UNIT) MASS.**
C    ** THIS ROUTINE IS GENERAL, AND CAN BE USED FOR ATOMS, LINEAR    **
C    ** MOLECULES, AND NON-LINEAR MOLECULES.                          **
C    **                                                               **
C    ** ROUTINE REFERENCED:                                           **
C    **                                                               **
C    ** REAL FUNCTION GAUSS ( DUMMY )                                 **
C    **    RETURNS A UNIFORM RANDOM NORMAL VARIATE FROM A             **
C    **    DISTRIBUTION WITH ZERO MEAN AND UNIT VARIANCE.             **
C    *******************************************************************

        INTEGER     N
        PARAMETER ( N = 108 )

        REAL        RX(N), RY(N), RZ(N), EX(N), EY(N), EZ(N)
        REAL        VX(N), VY(N), VZ(N), OX(N), OY(N), OZ(N)
        REAL        TEMP

        REAL        RTEMP, SUMX, SUMY, SUMZ
        REAL        GAUSS, DUMMY
        INTEGER     I

C    *******************************************************************

        RTEMP = SQRT ( TEMP )

        DO 100 I = 1, N

           VX(I) = RTEMP * GAUSS ( DUMMY )
           VY(I) = RTEMP * GAUSS ( DUMMY )
           VZ(I) = RTEMP * GAUSS ( DUMMY )

100     CONTINUE

C    ** REMOVE NET MOMENTUM **

        SUMX = 0.0
        SUMY = 0.0
        SUMZ = 0.0

        DO 200 I = 1, N

           SUMX = SUMX + VX(I)
           SUMY = SUMY + VY(I)
           SUMZ = SUMZ + VZ(I)

200     CONTINUE

        SUMX = SUMX / REAL ( N )
        SUMY = SUMY / REAL ( N )
        SUMZ = SUMZ / REAL ( N )

        DO 300 I = 1, N

           VX(I) = VX(I) - SUMX
           VY(I) = VY(I) - SUMY
           VZ(I) = VZ(I) - SUMZ

300     CONTINUE

        RETURN
        END



        SUBROUTINE ANGVEL ( TEMP, INERT )

        COMMON / BLOCK1 / RX, RY, RZ, EX, EY, EZ,
     :                    VX, VY, VZ, OX, OY, OZ

C    *******************************************************************
C    ** ANGULAR VELOCITIES FROM THE MAXWELL-BOLTZMANN DISTRIBUTION.   **
C    **                                                               **
C    ** THE DISTRIBUTION IS DETERMINED BY TEMPERATURE AND INERTIA.    **
C    ** THIS ROUTINE IS SPECIFIC TO LINEAR MOLECULES.                 **
C    ** IT CHOOSES THE DIRECTION OF THE ANGULAR VELOCITY RANDOMLY BUT **
C    ** PERPENDICULAR TO THE MOLECULAR AXIS. THE SQUARE OF THE        **
C    ** MAGNITUDE OF THE ANGULAR VELOCITY IS CHOSEN FROM AN           **
C    ** EXPONENTIAL DISTRIBUTION. THERE IS NO ATTEMPT TO SET THE      **
C    ** TOTAL ANGULAR MOMENTUM TO ZERO.                               **
C    **                                                               **
C    ** ROUTINE REFERENCED:                                           **
C    **                                                               **
C    ** REAL FUNCTION RANF ( DUMMY )                                  **
C    **    RETURNS A UNIFORM RANDOM VARIATE ON THE RANGE ZERO TO ONE  **
C    *******************************************************************

        INTEGER     N
        PARAMETER ( N = 108 )

        REAL        RX(N), RY(N), RZ(N), EX(N), EY(N), EZ(N)
        REAL        VX(N), VY(N), VZ(N), OX(N), OY(N), OZ(N)
        REAL        TEMP, INERT

        REAL        NORM, DOT, OSQ, O, MEAN
        REAL        XISQ, XI1, XI2, XI
        REAL        RANF, DUMMY
        INTEGER     I

C       ****************************************************************

        MEAN = 2.0 * TEMP / INERT

C    ** SET DIRECTION OF THE ANGULAR VELOCITY **

        DO 100 I = 1, N

C       ** CHOOSE A RANDOM VECTOR IN SPACE **

           XISQ = 1.0

1000       IF ( XISQ .GE. 1.0 ) THEN

              XI1  = RANF ( DUMMY ) * 2.0 - 1.0
              XI2  = RANF ( DUMMY ) * 2.0 - 1.0
              XISQ = XI1 * XI1 + XI2 * XI2

              GO TO 1000

           ENDIF

           XI    = SQRT ( 1.0 - XISQ )
           OX(I) = 2.0 * XI1 * XI
           OY(I) = 2.0 * XI2 * XI
           OZ(I) = 1.0 - 2.0 * XISQ

C       ** CONSTRAIN THE VECTOR TO BE PERPENDICULAR TO THE MOLECULE **

           DOT   = OX(I) * EX(I) + OY(I) * EY(I) + OZ(I) * EZ(I)
           OX(I) = OX(I) - DOT * EX(I)
           OY(I) = OY(I) - DOT * EY(I)
           OZ(I) = OZ(I) - DOT * EZ(I)

C       ** RENORMALIZE **

           OSQ   = OX(I) * OX(I) + OY(I) * OY(I) + OZ(I) * OZ(I)
           NORM  = SQRT ( OSQ )
           OX(I) = OX(I) / NORM
           OY(I) = OY(I) / NORM
           OZ(I) = OZ(I) / NORM

C       ** CHOOSE THE MAGNITUDE OF THE ANGULAR VELOCITY **

           OSQ   = - MEAN * LOG ( RANF ( DUMMY ) )
           O     = SQRT ( OSQ )
           OX(I) = O * OX(I)
           OY(I) = O * OY(I)
           OZ(I) = O * OZ(I)

100     CONTINUE

        RETURN
        END



        REAL FUNCTION GAUSS ( DUMMY )

C    *******************************************************************
C    ** RANDOM VARIATE FROM THE STANDARD NORMAL DISTRIBUTION.         **
C    **                                                               **
C    ** THE DISTRIBUTION IS GAUSSIAN WITH ZERO MEAN AND UNIT VARIANCE.**
C    **                                                               **
C    ** REFERENCE:                                                    **
C    **                                                               **
C    ** KNUTH D, THE ART OF COMPUTER PROGRAMMING, (2ND EDITION        **
C    **    ADDISON-WESLEY), 1978                                      **
C    **                                                               **
C    ** ROUTINE REFERENCED:                                           **
C    **                                                               **
C    ** REAL FUNCTION RANF ( DUMMY )                                  **
C    **    RETURNS A UNIFORM RANDOM VARIATE ON THE RANGE ZERO TO ONE  **
C    *******************************************************************

        REAL        A1, A3, A5, A7, A9
        PARAMETER ( A1 = 3.949846138, A3 = 0.252408784 )
        PARAMETER ( A5 = 0.076542912, A7 = 0.008355968 )
        PARAMETER ( A9 = 0.029899776                   )

        REAL        SUM, R, R2
        REAL        RANF, DUMMY
        INTEGER     I

C    *******************************************************************

        SUM = 0.0

        DO 10 I = 1, 12

           SUM = SUM + RANF ( DUMMY )

10      CONTINUE

        R  = ( SUM - 6.0 ) / 4.0
        R2 = R * R

        GAUSS = (((( A9 * R2 + A7 ) * R2 + A5 ) * R2 + A3 ) * R2 +A1 )
     :          * R

        RETURN
        END



        REAL FUNCTION RANF ( DUMMY )

C    *******************************************************************
C    ** RETURNS A UNIFORM RANDOM VARIATE IN THE RANGE 0 TO 1.         **
C    **                                                               **
C    **                 ***************                               **
C    **                 **  WARNING  **                               **
C    **                 ***************                               **
C    **                                                               **
C    ** GOOD RANDOM NUMBER GENERATORS ARE MACHINE SPECIFIC.           **
C    ** PLEASE USE THE ONE RECOMMENDED FOR YOUR MACHINE.              **
C    *******************************************************************

        INTEGER     L, C, M
        PARAMETER ( L = 1029, C = 221591, M = 1048576 )

        INTEGER     SEED
        REAL        DUMMY
        SAVE        SEED
        DATA        SEED / 0 /

C    *******************************************************************

        SEED = MOD ( SEED * L + C, M )
        RANF = REAL ( SEED ) / M

        RETURN
        END
赞一下 回复此楼
9楼2009-04-11 17:43:24
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

jianchaoyv

金虫 (小有名气)
谢谢！！！太感谢了！！！！！
回复此楼
10楼2009-04-11 18:39:54
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 jianchaoyv 的主题更新
141/2返回列表12下一页
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
信息提示
关闭
请填处理意见
关闭