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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » 【求助】求教fortran 90高手！！
北京石油化工学院2026年研究生招生接收调剂公告
51/1返回列表
查看: 877  |  回复: 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 个回答
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 22408 359分调剂 +4 Qshers 2026-03-27 5/250 2026-03-28 21:26 by zhq0425
[考研] 083000学硕274求调剂 +8 Li李鱼 2026-03-26 8/400 2026-03-28 20:33 by 加油向未来啊
[考研] 286求调剂 +12 PolarBear11 2026-03-26 12/600 2026-03-28 12:14 by zllcz
[考研] 311求调剂 +3 希望上岸阿小杨 2026-03-23 3/150 2026-03-28 07:57 by 热情沙漠
[考研] 315分求调剂 +7 26考研上岸版26 2026-03-26 7/350 2026-03-28 04:05 by fmesaito
[考研] 0856调剂 +5 求求让我有书读� 2026-03-26 6/300 2026-03-27 15:12 by caszguilin
[考研] 348求调剂 +4 小懒虫不懒了 2026-03-27 5/250 2026-03-27 12:47 by 果果妈咪
[论文投稿] Journal of Mechanical Science and Technology +3 Russ_ss 2026-03-25 5/250 2026-03-27 10:49 by 陆小果画大饼
[考研] 一志愿吉大071010，316分求调剂 +3 xgbiknn 2026-03-27 3/150 2026-03-27 10:36 by guoweigw
[考研] 0703化学338求调剂！ +6 Zuhui0306 2026-03-26 7/350 2026-03-27 10:35 by shangxh
[考研] 求调剂323材料与化工 +7 1124361 2026-03-24 7/350 2026-03-27 10:22 by wangjy2002
[考研] 材料调剂 +8 匹克i 2026-03-23 8/400 2026-03-27 08:11 by hypershenger
[考研] 329求调剂 +5 1() 2026-03-22 5/250 2026-03-26 20:40 by fmesaito
[考研] 281求调剂 +3 亚克西good 2026-03-26 5/250 2026-03-26 19:48 by 不吃魚的貓
[考研] 296求调剂 +4 汪！？！ 2026-03-25 7/350 2026-03-25 16:41 by 汪！？！
[考研] 【2026考研调剂】制药工程 284分 求相关专业调剂名额 +4 袁奂奂 2026-03-25 8/400 2026-03-25 14:32 by lbsjt
[考研] 284求调剂 +15 Zhao anqi 2026-03-22 15/750 2026-03-25 12:51 by wht0531
[考研] 282求调剂 +3 wcq131415 2026-03-24 3/150 2026-03-25 12:16 by userper
[考博] 26申博自荐 +3 whh869393 2026-03-24 3/150 2026-03-24 09:55 by 21018060
[考研] 336化工调剂 +4 王大坦1 2026-03-23 5/250 2026-03-23 18:32 by allen-yin
信息提示
关闭
请填处理意见
关闭