| 查看: 885 | 回复: 13 | |||
| 当前主题已经存档。 | |||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | |||
[交流]
【求助】求教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 ] |
回复此楼» 猜你喜欢
085601材料工程找调剂
已经有11人回复
-
总分293求调剂
已经有9人回复
-
334分 一志愿武理 材料求调剂
已经有5人回复
-
327求调剂
已经有6人回复
-
329求调剂
已经有10人回复
-
275求调剂
已经有15人回复
-
291求调剂
已经有6人回复
-
309求调剂
已经有7人回复
-
297求调剂
已经有11人回复
-
一志愿南航 335分 | 0856 | GPA 4.07 | 有科研经历
已经有6人回复
| 我有这本书《computer simulation of liquids》中后面有37个程序，我都没有，请问在哪能搞到？书中提到的网址:http://www.dl.ac.uk/CCp/ccp5/main.html怎么打不开啊？如果你方便的话能否把你拥有的《computer simulation of liquids》中提到的37个程序给我一些，我的QQ：526068505 |
11楼2009-04-11 20:48:35
老虎大王
木虫 (著名写手)
- 应助: 26 (小学生)
- 贵宾: 0.17
- 金币: 4774.1
- 散金: 8
- 红花: 42
- 帖子: 1361
- 在线: 215.2小时
- 虫号: 659094
- 注册: 2008-11-21
- 专业: 金属结构材料
★ ★ ★
kuhailangyu(金币+3,VIP+0):替虫友致谢！ 4-10 20:13
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 ] |
2楼2009-04-10 20:12:12
老虎大王
木虫 (著名写手)
- 应助: 26 (小学生)
- 贵宾: 0.17
- 金币: 4774.1
- 散金: 8
- 红花: 42
- 帖子: 1361
- 在线: 215.2小时
- 虫号: 659094
- 注册: 2008-11-21
- 专业: 金属结构材料
★ ★ ★
jianchaoyv(金币+3,VIP+0):问题没全部解决,但是非常感谢你这位绝对是高手了!!! 4-10 21:25
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 |
3楼2009-04-10 20:24:50
老虎大王
木虫 (著名写手)
- 应助: 26 (小学生)
- 贵宾: 0.17
- 金币: 4774.1
- 散金: 8
- 红花: 42
- 帖子: 1361
- 在线: 215.2小时
- 虫号: 659094
- 注册: 2008-11-21
- 专业: 金属结构材料
★ ★ ★ ★ ★ ★ ★ ★
kuhailangyu(金币+1,VIP+0):追着感谢！ 4-10 21:59
jianchaoyv(金币+7,VIP+0):感谢帮忙,交个朋友吧,以后多多探讨! 4-11 09:03
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 ] |
5楼2009-04-10 21:50:00
