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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » 【求助】求教fortran 90高手！！
51/1返回列表
查看: 817  |  回复: 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的回帖

老虎大王

木虫 (著名写手)
★ ★ ★
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的回帖
查看全部 14 个回答

老虎大王

木虫 (著名写手)
★ ★ ★
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的回帖

老虎大王

木虫 (著名写手)
★ ★ ★ ★ ★ ★ ★ ★
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的回帖
查看全部 14 个回答
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
信息提示
关闭
请填处理意见
关闭